Strain Localization Modes within Single Crystals Using Finite Deformation Strain Gradient Crystal Plasticity

The present paper aims at providing a comprehensive investigation of the abilities and limitations of strain gradient crystal plasticity (SGCP) theories in capturing different kinds of localization modes in single crystals. To this end, the small deformation Gurtin-type SGCP model recently proposed by the authors, based on non-quadratic defect energy and the uncoupled dissipation assumption, is extended to finite deformation. The extended model is then applied to simulate several single crystal localization problems with different slip system configurations. These configurations are chosen in such a way as to obtain idealized slip and kink bands as well as general localization bands, i.e., with no particular orientation with respect to the initial crystallographic directions. The obtained results show the good abilities of the applied model in regularizing various kinds of localization bands, except for idealized slip bands. Finally, the model is applied to reproduce the complex localization behavior of single crystals undergoing single slip, where competition between kink and slip bands can take place. Both higher-order energetic and dissipative effects are considered in this investigation. For both effects, mesh-independent results are obtained, proving the good capabilities of SGCP theories in regularizing complex localization behaviors. The results associated with higher-order energetic effects are in close agreement with those obtained using a micromorphic crystal plasticity approach. Higher-order dissipative effects led to different results with dominant slip banding.


Introduction
Strain localization is an important plastic instability process occurring prior to fracture. It is usually observed in the form of narrow bands of intense plastic shear strain in deformed bodies undergoing severe inhomogeneous deformation. This process has received strong scientific interest since the pioneering work of Considère [1], and numerous studies have been published on the subject [2][3][4][5][6][7][8][9]. Several types of localization modes have been revealed by these studies, depending on the microstructure of materials. Considering a single crystal with multiple activated slip systems, it has been shown that strain localization may occur along general shear bands, which can have an arbitrary orientation with respect to the crystallographic directions [3,10,11].
When only one slip system is activated, two types of shear bands, known as slip and kink bands, may occur at incipient plasticity according to the seminal work of Asaro and Rice [4] bifurcation analysis. Slip bands are sharp structures formed by intense dislocation slip on a few neighboring crystallographic planes. They are parallel to the active slip plane and correspond to the main signature of plastic slip in single crystals. Kink bands, which were first observed by Orowan [12], are localization bands orthogonal to the slip direction. They are accompanied by a strong lattice rotation in the bands compared to the surrounding regions. These bands are bounded by high lattice curvature boundaries formed by geometrically necessary dislocations (GNDs) walls.
These bands present a much larger width than slip bands. They are known to occur when strain incompatibility arises and there are not enough slip systems available to accommodate this incompatibility. According to the experimental evidence of Di Gioacchino and Quinta da Fonseca [13], such bands are characterized by the presence of a dense distribution of short slip bands along the slip plane normal direction. An interesting discussion about the occurrence of these bands can be found in [9].
Conventional crystal plasticity (CCP) theories have widely been applied in the literature to study strain localization within single crystals [5,[14][15][16][17]. Although these theories are able to capture several kinds of localization modes, including slip and kink bands, they present mesh-dependence difficulties. Localization bands obtained using such theories present mesh-dependent widths and plastic slip peak values [6,18]. In addition, CCP theories identically predict slip and kink bands, which appear in this framework as equivalent bifurcation modes.
This result was demonstrated by Asaro and Rice [4] using bifurcation analysis for single crystals undergoing single slip. A confirmation of such results for poly-crystals was recently provided by Marano et al. [8] using Fast Fourier Transform (FFT) poly-crystal simulations. In these simulations, systematically comparable distributions of slip and kink bands were obtained. This seems not to be in accordance with strain localization in real materials where overwhelming dominant slip bands are commonly observed. In the context of CCP theories, selection between slip and kink bands is only due to structural effects.
From a physical point of view, the equivalence between these bands is questioned as they rely on quite different formation processes. Compared to slip bands, more dislocation sources are generally required to form kink bands. Consequently, CCP theories are not suitable to study localization phenomena in single crystals. These theories include no internal length scale(s) allowing for stabilizing localization, which theoretically will occur in a set of zero measures.
A solution to overcome the limitations of the aforementioned theories consists in applying nonlocal plasticity approaches. Including internal length scale(s), these approaches provide a natural framework to capture nonlocal effects. Their promising features have motivated the development of various nonlocal models for both single and poly-crystal structures over the last two decades [6,8,[19][20][21][22][23][24][25][26][27][28]. Although most of these models have been developed to investigate the size effect at small scales, a few works applying them to localization problems exist in the literature.
In the context of single crystals, Forest [6] proposed a small deformation Cosserat crystal plasticity model to investigate the formation of slip and kink bands. Taking into account lattice rotation, the proposed model allows for breaking the equivalence between the two bands. Furthermore, it enables the regularization of kink bands, which present a finite width. However, this is not the case for slip bands, which are not necessarily accompanied by lattice curvature. Recently, a reduced micromorphic crystal plasticity model relying on the gradient of a cumulative scalar slip variable was proposed by Ling et al. [29].
This model presents the unique feature of regularizing both slip and kink bands. However, the same internal length scale is involved in the regularization of these bands, which could be problematic considering their different physical features and formation processes. The application of the proposed model to study strain localization within single crystals undergoing single slip revealed a complex localization behavior with competition between kink and slip bands. Experimental evidence of this behavior can be found in [17].
Another class of nonlocal approaches, which presents auspicious features to capture localization phenomena in single crystals, is the class of strain gradient crystal plasticity (SGCP) theories. This class has been the subject of a large number of recent works mostly focusing on size effects [11,24,[30][31][32][33][34]. However, only a few works applying SGCP theories to study localization phenomena can be found in the literature. Using a finite deformation SGCP model, Kuroda [11] investigated the effects of the finite element type on the formation of localization bands within single crystals. Scherer et al. [35] proposed a reduced SGCP model based on the work of Ling et al. [29] to study the formation of localization bands in voided irradiated materials.
This model relies on the gradient of a cumulative scalar slip variable, and this wa implemented using an approach combining micromorophic and Lagrange multiplier methods. Recently, Marano et al. [9] applied a small deformation version of the model of Gurtin [21], based on the full GND density tensor, to study the influence of the strain gradient plasticity on the modeling of localization bands within poly-crystal solids. In almost all existing studies of localization phenomena in single crystals, only higher-order energetic effects have been considered. Higher-order dissipative effects on these phenomena, particularly the competition between slip and kink bands, have not yet been explored.
The present paper aims at tackling these tasks by providing a comprehensive investigation of the abilities and limitations of SGCP theories in capturing different localization modes while considering both higher-order energetic and dissipative effects. To this end, the small deformation Gurtin-type SGCP model presented in [34], based on non-quadratic defect energy and uncoupled dissipation assumption, is extended to a finite deformation framework following the procedure detained in [21,23].
The extended model is then applied to simulate various localization problems implying single crystals with different slip system configurations. A focus is given to the ability of this model to capture the complex localization behavior characterizing single crystals undergoing single slip. The effects of higher-order energetic and dissipative stresses on this behavior are discussed.
Following this introduction, the paper is organized as follows. Section 2 details the theoretical formulation of the proposed Gurtin-type finite deformation SGCP model. A two-dimensional (2D) version of this model is derived and implemented in Section 3. Section 4 validates this implementation using simple shear tests of a constrained layer, for which small and finite deformation results exist in the literature. In Section 5, the 2D version of the proposed SGCP model is applied to key localization problems in single crystals. Finally, the paper ends with our concluding remarks in Section 6.

Finite Deformation Strain Gradient Crystal Plasticity (SGCP) Model
In this section, the small deformation Gurtin-type SGCP model proposed by Jebahi et al. [34], based on generalized defect energy and uncoupled dissipation, is extended to the finite deformation framework. Following [21,23], the finite deformation model is formulated with respect to the reference configuration to simplify its numerical implementation.

Kinematics and Single Crystal Hypothesis
As in [21,23], within the framework of finite deformation, the elastic-plastic kinematics are described by the Kröner-Lee decomposition of the deformation gradient F [36,37]: where F e represents the elastic distortion including stretch and rotation of the lattice and F p represents the plastic distortion due to plastic deformation caused by the flow of defects. The plastic incompressibility is also assumed with det(F p ) = 1. Considering (1), the velocity gradient tensor can be written: with the elastic and plastic distortion rate tensors L e and L p defined by: The single crystal hypothesis requires L p to be governed by slip rates on the individual slip systems. These systems are defined in the intermediate space, which is also called lattice or structural space. In this space, the slip direction and the slip plane normal direction associated with a slip system α, respectively noted s α and m α , are constant lattice vectors: s α · m α = 0 and |s α | = |m α | = 1 (4) These vectors as well as the glide direction vector (l α = m α ∧ s α ) are expressed in the reference configuration as [38]: Using single crystal hypothesis, the plastic distortion rate tensor L p can be written within the intermediate space as: where S α = s α ⊗ m α denotes the Schmid tensor associated with the slip system α,γ α represents the plastic slip rate on this system, and q is the number of slip systems. For the sake of consistency with the small deformation SGCP model [34], the same symbol (γ α ) is used to designate the plastic slip rate on the slip system α. However, within the finite deformation framework,γ α cannot be interpreted as the material time derivative of a physically-meaningful field variable γ α . Indeed, slips (as fields on individual slip systems) are not generally well-defined quantities from a physical point of view [23].

Principle of Virtual Power and Governing Equations
Considering the generalized power density of internal forces [39,40], the internal and external virtual power expenditures within an arbitrary subregion V 0 of boundary S 0 can be expressed in the reference configuration as (neglecting body forces): where P = JσF T−1 is the first Piola-Kirchhoff stress tensor (power-conjugate toḞ), J = det(F) is the volumetric Jacobian, σ is the Cauchy stress tensor, κ α is a scalar microscopic stress (power-conjugate toγ α ), ξ α is a vector microscopic stress (power-conjugate to grad X (γ α )), t is a macroscopic surface traction, χ α is a microscopic surface traction, and grad X (•) is the Lagrangian gradient operator (i.e., w.r.t. the Lagrangian coordinates). This operator can be related to the lattice gradient operator grad # (•) and the Eulerian gradient operator grad x (•) by: Introducing the first and second Piola elastic stress tensors, which are, respectively, defined by: P e = JσF eT−1 and S e = J F e−1 · σ · F eT−1 (10) and defining the resolved shear (Schmid) stress τ α as: with C e = F eT F e is the right Cauchy-Green elastic deformation tensor. The term P :Ḟ in the expression of P int (7) can be rewritten as: which leads to: with the scalar microscopic stress π α = κ α + τ α . Application of the generalized principle of virtual power [41,42] results in the following macroscopic and microscopic balance equations and boundary conditions: with n 0 the outward unit normal to S 0 and div X (•) the Lagrangian divergence operator.

Free Energy Imbalance
The free energy imbalance (second thermodynamic law) requires that the temporal increase of free energy in an arbitrary subregion V 0 be less than or equal to the external power expended on V 0 . Mathematically, this can be formulated as: where Ψ denotes the free energy measured in the reference configuration. Defining the elastic Green-Lagrange strain tensor E e as: it can easily be shown that P e :Ḟ e = S e :Ė e . Since V 0 is arbitrary, the local free energy imbalance can be obtained:

Energetic Constitutive Laws
The free energy in (16) consists of an elastic strain energy Ψ e and a defect energy Ψ ρ . Building on the work of Gurtin [23], Ψ ρ is assumed to be a function of dislocation densities: where ρ α and ρ α are, respectively, the edge and screw dislocation densities associated with the slip system α, and are, respectively, defined by: These densities are in units of length −1 and may be positive or negative. Bearing in mind Equations (5) and (9), it can easily be demonstrated that: The elastic strain energy Ψ e is assumed to be a quadratic form of E e : where C is a symmetric and positive-definite elasticity tensor. In rate form, this quantity can be written as: As in [34], a generalized power-law defect energy with adjustable order-controlling index n is adopted in this work: where X 0 is a constant representing the energetic slip resistance, and l en is an energetic length scale. To ensure the convexity of Ψ ρ , the defect energy index n must be greater than or equal to 1 (n ≥ 1). Its rate form is given by: Substituting (19) and (21) into (16), the local free energy imbalance can be rewritten as: In the following, it is assumed that S e is fully energetic, π α is fully dissipative, and ξ α is divided into energetic and dissipative parts: ∂Ψ e ∂E e , π α = π α dis , and ξ α = ξ α en + ξ α dis It should be noted that higher-order dissipation in the context of strain gradient (crystal) plasticity theories still remains an open question. Although ignored in some works, several authors consider this kind of dissipation in their models [34,[43][44][45][46]. As reported by Fleck and Willis [45], direct experimental measurement shows that the core energy of dislocations stored during plastic deformation is much smaller than the plastic work dissipated in dislocation motion.
Consequently, both statistically stored and geometrically necessary dislocations would contribute more to plastic dissipation than to a change in free energy. The source of the dissipation process can be attributed to the motion of the mobile dislocations. The motion of dislocations yields energy dissipation, which, in turn, results in resistance to the dislocation motion and increases in the yield strength due to size effects [44].
The energetic constitutive laws can then be obtained: It can be verified from (22) 2 that the energetic microscopic stress ξ α en is tangent to the αth slip plane. Using (22), the dissipation inequality becomes: ξ α dis characterizes the dissipative microscopic forces associated with the evolution of dislocations in the αth slip plane. As the motion of such dislocations is tangent to this plane, ξ α dis is also required to be tangential. Therefore, without losing generality, the gradient term in (23) can be replaced by the corresponding tangential gradient defined by:

Dissipative Constitutive Laws
Referring to [34], uncoupled effective plastic strain measures are used to describe the first-and higher-order dissipative effects. These measures, which are referred to hereafter as first-and higher-order measures, are, respectively, defined by: where l dis is dissipative internal length scale. To ensure positive dissipation, the scalar and vector dissipative stresses associated with slip system α are taken as: whereγ α 0 > 0 is a constant strain rate representative of the flow rates of interest, m > 0 is a constant characterizing the rate-sensitivity of the considered material, and S α i (i ∈ {π, ξ}) are stress-dimensional internal state variables (S α π > 0 and S α ξ ≥ 0). These variables are called, hereafter, dissipative slip resistances. Their evolutions are assumed to be governed by: are positive hardening moduli.

Flow Rule
Based on the expressions of the microscopic stresses, the microscopic balance equation associated with slip system α can be reformulated as follows: Substituting (22), (25) and (26) into (28), the flow rule associated with slip system α can be obtained: Term I in the above equation, being energetic represents a backstress, which leads to kinematic hardening (Bauschinger) effects. Terms II and III are dissipative and lead to dissipative hardening effects (a reversal of the flow direction, i.e.,γ α → −γ α , simply changes the sign of these terms).

Microscopic Boundary Conditions
Within the framework of higher-order gradient theories, the presence of higher-order stresses can result in an expenditure of microscopic power on the boundary of the studied domain (S 0 ). For each slip system α, the microscopic power expended, per unit area, on this boundary is given by: The present work focuses on simple microscopic boundary conditions (BCs) that result in a null power expenditure on S 0 : Specifically, two limit cases of microscopic BCs are considered in this work: with S h and S f complementary subsurfaces of S 0 . The micro-hard BCs (32) 1 correspond to impenetrable surfaces with no slip, whereas the micro-free BCs (32) 2 represent free surfaces without slip resistance.

Numerical Implementation
In this section, a two-dimensional version of the finite deformation SGCP model is implemented under plane strain conditions with planar slip systems, using a total Lagrangian finite element formulation. An extension to the generalized three-dimensional case can easily be derived.

Two-Dimensional Version of the Proposed Model
Assuming plane strain conditions with planar slip systems, it can easily be demonstrated that screw dislocation densities vanish [34,43,47]: Therefore, the energetic microscopic stresses become: The dissipative microscopic stresses are given by (25) and (26), with the slip resistances S α π and S α ξ assumed, hereafter, to evolve linearly with the effective plastic strain measures: Note that i ∈ {π, ξ}, S π0 > 0 and S ξ0 ≥ 0. Using the above expressions of the microscopic stresses, a simplified form of the flow rule (29) can be obtained: In summary, key equations characterizing the simplified two-dimensional SGCP model are recalled below: These equations will be solved numerically within a total Lagrangian finite element framework. Several methods are used in the literature to solve these equations. Among them, one can cite the approach proposed by Scherer et al. [48] combining the advantages of micromorphic and Lagrange multiplier methods. In this approach, a scalar micromorphic variable is linked to a cumulative plastic slip measure by a Lagrange multiplier to ensure uncoupling between constitutive nonlinearily and spatial nonlocality.
Other methods based on an independent thermodynamic variable for each slip system can be found, allowing for more flexible control of gliding on the slip systems. Two common examples of these methods can be cited: the ρ−method [11,[49][50][51][52] and the γ−method [34,43,53]. The former considers dislocation densities or their rates as primary variables (additional degrees of freedom), whereas plastic slips or plastic slip rates are considered as primary variables in the latter.
According to Klusemann et al. [50], the two methods generally lead to similar results. However, γ−methods provide more flexibility in constraining the dislocation flow on the boundaries. Furthermore, it does not require a local loop (at the integration points) to calculate the plastic slip ratesγ α , which are given by solving the global system of equations. This method is chosen in the present work. In addition to the displacement field u, the scalar field variables defined on each slip system α by: are taken as primary variables, in consistency with the implemented small deformation version of the proposed SGCP model [34]. It is worth recalling that these variables, which can be interpreted in small deformation as plastic slips on the slip systems, have no real physical interpretation in finite deformation (only their increments are physically relevant).
As will be seen in the validation section, this choice of primary variables has no impact on the numerical results, as the integration of the constitutive laws is only based on the increments of these variables.

Integration of Constitutive Laws
Numerical integration of constitutive laws consists in updating the values of stresses and internal variables, given the increments of kinematic variables. Knowing the increments of the primary field variables (∆u and ∆γ α ), the displacement, the deformation gradient, and the plastic distortion rate quantities can be calculated at a given point as: with I as the second order unity tensor. The plastic distortion tensor can be determined based on Equation (3) 2 . Using the implicit approach, this equation can be written in discretized form as: which leads to: To ensure plastic incompressibility (det F p = 1), a unimodular operator is applied to the term between brackets: where, for a second order tensor A, the term uni(A) is defined by: Considering (1), the elastic part of deformation gradient is calculated by: Using the above equation, the elastic Green-Lagrange strain tensor E e n+1 can be obtained: Knowing the updated kinematic quantities, the stresses and internal variables can be updated. Algorithm 1 summarizes the steps of integration of constitutive laws, where the subscripts "n" and "n + 1" are omitted for clarity.

Weak Forms and Finite Element Discretization
Assume the virtual fields δu and δγ α , which are kinematically admissible to zero on the portions of the boundary of the studied domain on which macroscopic and microscopic Dirichlet boundary conditions are respectively imposed. The weak forms of the macroscopic and microscopic balance equations (first lines of (12) and (13)) can be obtained as follows: where S t 0 and S χ α 0 are the portions of the domain boundary on which macroscopic and microscopic traction forces, noted, respectively, as t and χ α , are imposed. These weak forms are solved numerically using the total Lagrangian finite element method. To this end, a User-ELement (UEL) subroutine is implemented within the commercial finite element package ABAQUS/Standard. In this UEL, displacement and plastic scalar fields (u and γ α ) are considered as primary variables (degrees of freedom).
Several finite element types can be used to approximate these fields. In this work, an element type employing eight-noded (quadratic) interpolation of u, four-noded (linear) interpolation of γ α and 2 × 2 Gauss integration is adopted (Figure 1). Based on Kuroda [11], this element type is well suited to study localization problems. Within each element, dis-placements u i (i = 1, 2) and plastic scalar variables γ α (α = 1, 2, · · · q) can be approximated based on the associated nodal values as follows: where U i k and Γ α k are nodal values of displacements and plastic scalar variables, and N k u and N k γ are the associated interpolation functions.
Algorithm 1 Integration of constitutive equations. Using the above field approximations, the macroscopic and microscopic weak forms can be written within a representative finite element, in a simplified matrix form, as follows: where N u and B u are the interpolation and gradient matrices associated with the displacement field, N γ and B γ are the interpolation and gradient matrices associated with the plastic scalar fields, P, ξ, π are vector representations of the macroscopic stress and microscopic stresses on all slip systems, τ is vector of resolved shear stresses on all slip systems, t and χ are vectors of macroscopic and microscopic traction forces.
Considering the principle of virtual power (G e u = 0 and G e γ = 0, for any arbitrary δU e and δΓ e ), the associated residuals can be written: Linearization of these equations with respect to small variations of U e and Γ e leads to an elementary system of equations, which can be written in matrix form as: with The global system of linear equations can be obtained by assembling all the elementary systems of equations associated with the overall finite elements. This system is solved by means of a Newton-Raphson iterative solution scheme. At each iteration, updated values of the overall increments of the displacement and plastic scalar fields (∆U and ∆Γ, respectively) are obtained and used to numerically solve the constitutive equations at the Gauss points (Algorithm 1).

Validation of the Model Implementation
To ensure that the finite deformation SGCP model is correctly implemented, it is applied to simulate a simple shear problem of a constrained layer under plane strain conditions. The considered layer is first subjected to small shear loading, allowing for comparison with the small deformation version of the proposed model [34]. This makes it possible to draw conclusions about the effectiveness of the model implementation to correctly reproduce different energetic and dissipative size effects, for which large deformation results are not necessarily available in the literature.
Then, large shear loading is applied for a particular case of quadratic defect energy (n = 2) with no higher-order dissipation (l dis = 0). For this case, the present model is equivalent to the model of Kuroda and Tvergaard [49], for which large deformation shear results are available. These results are used to assess the ability of the model implementation to correctly capture finite deformation effects. Figure 2 presents the geometrical model used in this section. A strip of infinite width in the x 1 −direction and height H in the x 2 −direction is subjected to simple shear loading under plane strain conditions. Two slip systems symmetrically oriented with respect to the x 1 −direction (θ 1 = 60°and θ 2 = −60°) are considered. The macroscopic boundary conditions are given by: where Γ is the prescribed shear strain. Micro-hard boundary conditions are imposed on the top and bottom edges to hinder dislocation flow at them (impenetrable to dislocations): To model the infinite width of the strip, periodic boundary conditions are imposed at its left and right edges: With these conditions, the strip width W does not affect the results and can arbitrarily be chosen. To save computation time, W is chosen to have only one finite element along the x 1 −direction. A numerical sample consisting of a column of 100 finite elements in the x 2 −direction is then adopted. For the case of small shear loading, a complete loading cycle is imposed (i.e., 0 → Γ → −Γ → Γ), with the prescribed shear strain restricted to Γ = 0.02. The material parameters are chosen as in [34]. These parameters are presented in Table 1. Figure 3 compares the associated results obtained using the present finite deformation model with those obtained using the small deformation one [34], for different energetic and dissipative parameters. Solid line curves correspond to the finite deformation results and those with markers correspond to the small deformation ones.
For comparison purposes, results obtained using conventional small and finite deformation crystal plasticity with no hardening (H π = 0) are also presented in Figure 3 (curves named "Perfect CP"). It can be seen that the stress-strain responses obtained using the finite deformation formulation show very good agreement with those obtained using the small deformation one. Note that the stress-strain curves obtained using n = 1.7 present inflection points, reflecting nonlinear kinematic hardening type III (KIII) of Asaro [54]. Further discussions about this hardening type can be found in [34,55,56]. Applying small shear loading, the nonlinear geometrical effects are insignificant. In this case, the small and finite deformation formulations must lead to similar results, which is the case according to Figure 3. This validates the capabilities of the implemented finite deformation model in correctly reproducing the energetic and dissipative effects, at least in the range of the prescribed shear loading.  To examine the effectiveness of the present model implementation with regard to nonlinear geometrical effects, the simple shear example of [49], which implies large shear loading, is reproduced in this work. In this example, the same geometrical sample as in Figure 2 is used, with imposed shear loading Γ = 0.24. The defect energy index and the dissipative length scale are set as n = 2 and l dis = 0, respectively. With these choices, the present finite deformation model is equivalent to that of [49], in which only energetic back stress is considered.
The other material parameters, which are taken from [49], are presented in Table 2. Figure 4 compares the results of the present model with those obtained by [49] for two values of energetic length scale l en . Very good agreement between the results is obtained, regardless of the value of l en . This shows the ability of the present model implementation to correctly reproduce finite deformation effects. As can be seen from Figure 4a, in the range of small shear strains, until approximately Γ = 0.025, nonlinear geometrical effects are insignificant and the shear stress σ 12 increases almost linearly with Γ after yielding. Beyond approximately Γ = 0.025, nonlinear stressstrain dependence takes place and becomes more striking at larger strains, due to nonlinear geometrical effects. Figure 4b presents the distributions of the total plastic variables γ α associated with the considered slip systems along x 2 −direction, as obtained using the present and the [49] models.
As shown, the nonlinear geometrical effects lead to asymmetric evolution of these variables, despite the initial symmetry of the slip systems with respect to the loading direction. A distinct deviation between the γ 1 and γ 2 profiles, which coincide with those obtained by [49], is observed. Such a deviation was also reported by Bittencourt [32] using a similar simple shear problem under monotonic large shear loading. This can be explained by an apparent hardening increase in one system and an apparent hardening decrease in the other system due to lattice rotation effects [32,49].

Application to Strain Localization in Single Crystals
Several works studying strain localization in single crystals undergoing single and multiple slips can be found in the literature [3,4,6,[8][9][10]. Using multiple slip systems, it has been shown that strain localization may occur along narrow shear bands, which can be arbitrarily oriented with respect to the initial crystallographic directions [3,10,11]. In the case of a single slip system, two types of shear bands, known as slip and kink bands, can take place [4]. Slip bands correspond to slip localization along a crystallographic slip plane, whereas kink bands are orthogonal to the slip direction ( Figure 5). Due to lack of internal length scale(s), conventional crystal plasticity theories are not able to correctly capture localization bands. As an alternative, Forest [6] has applied a small deformation Cosserat single crystal plasticity theory to investigate the formation of slip and kink bands within single crystals. As shown by the author, this theory allows for regularizing kink bands that will have a finite thickness. However, it does not regularize slip bands, which are not necessarily associated with the lattice curvature [6,57].
Although not well investigated in the literature, SGCP theories, which are classically based on GND densities, are also expected to present difficulties in regularizing idealized slip bands. Indeed, a slip gradient along the normal to the slip planes does not result in the storage of GNDs. Ling et al. [29] proposed a finite deformation micromorphic crystal plasticity model that is capable of regularizing both kink and slip bands. Application of this model to strain localization in single crystals presenting one slip system has brought to light a complex localization behavior with competition between kink and slip bands. Recently, Marano et al. [9] applied a small deformation version of the SGCP model of Gurtin [21], which ignores higher-order dissipation to study localization phenomena within poly-crystal solids. Interesting discussion on the capabilities of the applied model in capturing the competition between slip and kink bands was provided in this study, using fully higher-order energetic effects.
In almost all existing studies of localization phenomena in single crystals, only higherorder energetic effects were considered. Higher-order dissipative effects on these phenomena, particularly the competition between slip and kink bands, have not yet been explored. The objective of the present section is to provide an in-depth investigation of the capabilities and limitations of the implemented finite deformation SGCP model in predicting different kinds of localization modes within single crystals, considering both higher-order energetic and dissipative effects. Particular focus will be given to the ability of such a model to reproduce the complex localization behavior encountered considering a single slip system, with competition between slip and kink bands.

Regularization of Localization Bands by SGCP Theories
Before applying the implemented SGCP model to study complex localization phenomena, the present subsection aims at reviewing its ability to capture different kinds of localization bands, including idealized slip and kink bands as well as general shear bands (i.e.,with no particular orientation with respect to the initial crystallographic directions).

Simple Shear with Single Slip System: Regularization of Idealized Slip/Kink Bands
Simple shear simulations using single crystals with single slip are performed in the present subsection. The geometrical model used in these simulations is similar to that presented in Figure 2, except for the slip systems, which are replaced by only one system. To obtain idealized slip and kink bands, this system is oriented such that the slip plane is either parallel (horizontal slip case) or normal (vertical slip case) to the shear direction.
The same microscopic and macroscopic boundary conditions as in Figure 2 are assumed. The material parameters are those presented in Table 1, except for H π , which is taken as H π = −500 MPa to model softening. To trigger strain localization, a material imperfection is considered in a small central region of height H/11. In this region, the initial first-order slip resistance S π0 is linearly decreased to 0.9 S π0 at the center.
Considering the case of horizontal slip, Figure 6 presents the associated results in terms of the stress-strain curves and distributions of plastic scalar variable γ along the x 2 −direction for different mesh refinements. Although these results are obtained considering only higher-order energetic effects (l en /H = 0.1 and l dis = 0), it has been verified that activation of higher-order dissipative effects (l dis = 0) do not affect them. The SGCP model failed to regularize the obtained idealized slip band, which shows mesh-dependence. As can be seen in Figure 6b, the width of the obtained band is set by the mesh spacing and decreases with increasing the number of finite elements.
For the considered slip configuration, the present SGCP model behaves as a conventional crystal plasticity model. Indeed, slip gradients in the direction normal to the slip plane do not correspond to GNDs [6,23,43]. Consequently, no size effects are predicted by GND-based SGCP theories for the development of such normal slip gradients. This is also the case for Cosserat crystal plasticity theories in which size effects are only triggered by lattice curvature [6].  Figure 7 presents the shear results associated with the case of vertical slip in terms of the stress-strain responses and distributions of plastic scalar variable γ along the */++x 2 −direction for different mesh refinements. Using a sufficiently fine mesh, the stress-strain response becomes nearly mesh-independent (Figure 7a). Note that the slope inversion and the plateau-shaped form observed on the stress-strain curves beyond Γ = 0.008 is due to the competition between softening caused by negative H π and hardening caused by higher-order energetic effects (l en = 0). As shown in Figure 7b, thick kink bands that spread over several finite elements are obtained. For a relatively fine mesh, the thickness of these bands is nearly mesh-insensitive. For Cosserat crystal plasticity models [6], the proposed SGCP model enables meshindependent simulation of kink bands, which imply slip gradients in the slip direction. On the contrary, it fails to regularize idealized slip bands, which are not accompanied by the development of GNDs. However, the idealized slip conditions used to obtain the latter bands are not commonly encountered in real material deformation. The following subsection discusses the ability of the present SGCP model to correctly capture strain localization within a more practical problem, for which shear bands do not necessarily have particular orientations with respect to the initial crystallographic directions.

Uniaxial Tension with Two Slip Systems: Regularization of General Localization Bands
A uniaxial tension problem with two slip systems symmetrically oriented with respect to the tension direction is considered in this subsection. Figure 8 presents the associated geometrical model, which is a plate of dimensions 2 L × 2 W = 20 × 10 mm having two slip systems inclined at 60°from both sides of x 1 −direction (θ 1 = −θ 2 = 60°). For this problem, according to bifurcation analysis [3,10], it is expected that localization will occur along two crossed shear bands whose orientations are close to those of the slip planes. Due to symmetry, only one quarter (upper right part) of the plate is analyzed. The macroscopic boundary conditions applied on this portion of the plate are set as: Furthermore, all the edges of this plate portion are assumed to be micro-free (transparent to dislocations). As in the previous simple shear problem, H π is assumed to be negative to model softening (H π = −500 MPa). The other material constants are given in Table 3, following several previous works designated for single crystal tension problems [11,15]. To trigger localization, a material imperfection is considered within a square region of 0.25 mm side located at the bottom left-hand corner of the considered plate portion (Figure 8). In this region, the initial first-order slip resistance S π0 is replaced by 0.9 S π0 . To quantify strain localization, the accumulated plastic slip variable defined by γ acc = ∑ α |γ α | is adopted.
Mat. imperf.   Figure 9 presents the associated results in terms of nominal stress-strain curves and accumulated plastic slip variable distributions along the x 2 −direction for different mesh refinements. The nominal stress is calculated as the reaction force F on the top edge over the initial width W 0 of the studied plate portion (the thickness is assumed to be 1) and the nominal strain as the ratio between the elongation ∆L and the initial length L 0 . The values of γ acc in Figure 9b are determined at nodes belonging to the center line of the studied portion of the plate (dashed blue line in Figure 8). Figure 10 presents the deformed meshes and contours of the accumulated plastic slip variable γ acc for different mesh refinements. Concerning the nominal stress-strain response, mesh convergence is obtained as soon as a sufficiently fine mesh is applied (Figure 9a). In terms of strain localization, it seems that finer meshes lead to a larger γ acc peak value in the obtained shear band. However, this band presents a finite width, which converges toward an intrinsic value covering several finite elements.
Based on Figure 10, such a band is oriented at an angle of approximately 54°with respect to the x 1 −direction, which is slightly smaller than the angle of the first slip direction (θ 1 = 60°). This is in agreement with the bifurcation and numerical results of Kuroda [11]. Using a similar configuration of activated slip systems with θ 1 = 54.7°, these authors demonstrated that strain localization occurs along a line inclined at approximately 49°from the x 1 −direction.
To better analyze the obtained results, the lattice rotation is approximated based on the polar decomposition of the elastic deformation gradient F e = R e U e , where R e and U e are, respectively, the elastic rotation and elastic stretch tensors. Neglecting the elastic distortion described by U e , R e can be interpreted as the lattice rotation, and the corresponding angle θ is calculated using [8]: The angle sign is determined considering the sign of the off-diagonal terms of R e . By doing so, we find that strain localization is accompanied by a lattice rotation of approximately −5.5°. Therefore, the localization band obtained in the considered portion of the plate is not an idealized slip band, as it is associated with lattice curvature, which explains its regularization by the present SGCP model. However, the obtained band could be considered as a slip band since it is nearly parallel to the rotated first slip direction within the band (60°− 5.5°= 54.5°).
Although GND-based SGCP theories fail to regularize idealized slip bands, the present investigation shows that they are capable of regularizing general localization bands. The next subsection aims at investigating the ability of such theories to correctly capture complex localization behaviors with competition between the kink and slip bands.

Investigation of Localization Modes in Single Crystals Undergoing Single Slip
The present finite deformation SGCP model is applied, hereafter, to reproduce the uniaxial tension problem of a single crystal undergoing single slip, which has been studied by Ling et al. [29] using a micromorphic crystal plasticity approach. This problem is characterized by complex localization behavior with competition between the kink and slip bands.
The geometrical model is similar to that used by Ling et al. [29]. It consists of a single crystal plate of dimensions W and L in, respectively, the x 1 − and x 2 − directions with an aspect ratio of L/W = 6 and a single slip system oriented to an angle of θ = −33.7°with respect to x 1 −direction. This plate is subjected to the following macroscopic boundary conditions: Furthermore, all the plate edges are assumed to be micro-free (transparent to dislocations). The material parameters are the same as those used in the previous uniaxial tension problem (Table 3), except for H π , which is set to a small positive value (H π = 1 MPa). For this problem, it is not needed to introduce material softening or imperfection to trigger strain localization. Indeed, the latter can occur due to geometrical softening induced by cross-section reduction and lattice rotation.
Although the material parameters used in the present and [29] investigations are not the same, comparison between the associated results is possible from a qualitative point of view, which is sufficient to meet the objectives of this subsection. As in [29], the plastic scalar variable defined by (37) on the considered slip system is adopted to investigate the strain localization and formation of slip and kink bands. Since the main difference between these bands is that kink bands are accompanied by a lot of lattice rotation, the lattice rotation angle (53) is used to detect the presence of the latter bands.

Regularization of the Mesh-Dependence of the Stress-Strain Response
To check the mesh-regularization capability of the present SGCP model regarding the stress-strain response, different mesh refinements are considered: 8 × 48, 16 × 96, and 24 × 144 meshes. Figure 11 presents the associated results obtained using different constitutive configurations. In this figure, the CP, SGCP-EN, and SGCP-DIS curves are, respectively, associated with the conventional CP (l en = l dis = 0), SGCP including only higher-order energetic effects (l en /W = 0.25 and l dis /W = 0), and SGCP including only higher-order dissipative effects ( l en /W = 0 and l dis /W = 0.25). For conventional crystal plasticity, the results are highly mesh-dependent and display spurious oscillations. Similar results were obtained by [29].
Activation of higher-order effects, regardless of their energetic or dissipative nature, allows for overcoming these mesh difficulties. No spurious oscillations were observed for the results associated with SGCP-EN and SGCP-DIS cases. Furthermore, mesh convergence is obtained for both cases as soon as a sufficiently refined mesh is used. The results associated with 16 × 96 and 24 × 144 meshes are in very close agreement for both these cases. The apparent increase of yield stress obtained for the SGCP-DIS case is due to the strengthening effects of the dissipative length scale l dis , as widely recognized in the literature [34,43]. In the following subsections, the considered plate is discretized using a 24 × 144 mesh.

Localization Modes Using Conventional Crystal Plasticity
For comparison purposes, localization modes in the considered plate are first revisited using conventional CP (l en = l dis = 0). Figure 12a shows the evolution of the plastic slip variable γ (37) within this plate at different nominal strains. To investigate the formation of kink bands, evolution of the lattice rotation angle within the plate is presented in Figure 12b. At the beginning of straining (∆L/L 0 ≈ 0.005), both slip and kink bands start to emerge (crossed bands near the top and bottom edges of the plate). This is in agreement with the bifurcation analysis of Asaro and Rice [4].
According to this analysis, slip and kink bands can simultaneously appear at incipient plasticity when using conventional CP. By increasing ∆L/L 0 , the kinking phenomenon becomes dominant, leading to the formation of two pronounced kink bands. These bands are accompanied by considerable lattice rotation, which results in a local increase of the Schmid factor and then a local decrease of the resolved shear stress (RSS) for the considered slip system. This geometrical softening mechanism acting in the kink regions promotes the development of kink bands.
As straining increases, strain localizes in the kink band located near the bottom edge of the plate. In addition, the direction of localization starts to shift gradually, which is due to formation of a slip band crossing the initial kink region. Due to hardening effects, the RSS in the kink region increases and can reach a level for which the glide along the "unrotated" slip direction becomes easier, hence, resulting in the appearance of the slip band.
It can be observed from Figure 12b that formation of this band does not give rise to lattice rotation. Finally, the formed slip band intensifies with straining, leading to final necking of the plate.

Localization Modes Using Purely Higher-Order Energetic Effects
The above complex localization behavior is studied, hereafter, using the present SGCP model with only higher-order energetic effects. Figure 13 presents the contours of the plastic slip variable γ and the lattice rotation angle θ within the studied plate at different nominal strains, with the energetic length scale l en = 0.25 W. Compared to conventional crystal plasticity results, two broader kink bands appear at the beginning of straining (∆L/L 0 = 0.006 in Figure 13a). These bands rapidly aggregate to form a larger kink region with concentrated lattice rotation at the center of the studied plate. In the meanwhile, a central slip band appears and starts to interact with this region as shown at straining ∆L/L 0 = 0.015.
The kink-slip transition region shows a much smoother form compared to that obtained using conventional CP (Figure 12a). As straining increases, the formed slip band becomes more and more intensified, leading to final necking of the plate.
(a) Evolution of the plastic slip variable γ.  To investigate the influence of the higher-order energetic effects on the present localization behavior, different values of l en were tested. Figure 14 presents the associated results in terms of plastic slip contours in the studied plate at straining ∆L/L = 0.006, which corresponds to the emergence of kink bands. The width of the obtained bands as well as the distance between them are functions of l en . By increasing this parameter, larger and closer kink bands are obtained. For sufficiently large l en , a single kink band located at the center of the considered plate is formed. Figure 15 presents the nominal tensile stress-strain curves and the plastic slip variable distributions along the plate length direction for different energetic length scales. Although no change in the initial yield is observed, the stress level after yielding increases with increasing l en due to strain gradient hardening effects (Figure 15a). In addition, larger l en leads to a larger final localization band with lower peak value of plastic slip γ in this band (Figure15b). This is in accordance with several investigations based on nonlocal approaches [26,29,53,58,59]. The energetic length scale behaves as a "localization limiter" and governs the width of localization bands.

Localization Modes Using Purely Higher-Order Dissipative Effects
The present subsection aims at investigating the higher-order dissipative effects on the localization behavior in the considered plate, while neglecting the higher-order energetic effects (l en = 0). One distinctive feature of the present SGCP model is that the higherorder dissipative effects are uncoupled from the first-order ones. The advantages of this uncoupled constitutive assumption have already been discussed in [34]. It is worth recalling here that such an assumption makes it possible to flexibly choose the initial value of the higher-order dissipative slip resistance S ξ0 . This parameter can also be set to zero in order to obtain progressive activation of the higher-order dissipative effects, as with the energetic ones involved in the previous subsection. In this case, a nonzero higher-order hardening modulus H ξ is necessary to activate such effects according to (34). Using S ξ0 = 0 and H ξ = 60.84 MPa, the obtained results (not presented in the paper) are qualitatively quite similar to those obtained in the above subsection.
Applying nonzero S ξ0 , finite higher-order dissipative stresses take place and instantaneously activate with the emergence of strain gradients. The influence of these stresses is investigated, hereafter, assuming S ξ0 = S π0 = 60.84 MPa, as in classical SGCP theories with coupled dissipation, while keeping H ξ = 60.84 MPa. Figure 16 presents the associated results in terms of the plastic slip variable and lattice rotation distributions in the plate at different nominal tensile strains.
These results are quite different from those obtained using conventional CP and SGCP, including progressive activation of higher-order effects. Instead of two kink bands, a large central region parallel to the slip direction appears at the beginning of straining (∆L/L 0 = 0.0035) with a nearly plateau-shaped distribution of γ. The width of this region progressively decreases with straining.
A less-pronounced kink region appears at the plate center and interacts with the initial slip region, as can be shown by analyzing the lattice rotation distributions in the plate (Figure 16b). The width of both these regions continues to reduce with increasing ∆L/L 0 until the final necking of the plate. The final localization band is oriented to an angle close to that of the slip direction. Compared to those obtained in the previous subsections, this band is much more diffuse. Using higher-order dissipative effects with nonzero S ξ0 , slip localization seems to outweigh kinking, although the latter phenomenon is not completely removed.   Effects of l dis on the localization response of the considered plate are given in Figure 17, which is obtained assuming S ξ0 = S π0 = 60.84 MPa and H ξ = 60.84 MPa. This parameter leads to similar effects as for l en regarding the increase of hardening after yielding and the growth of the width of the final localization band. In addition to these effects, l dis presents strengthening effects, which are illustrated by an increase of the initial yield.

Conclusions
This contribution focused on the capabilities of strain gradient crystal plasticity (SGCP) in modeling localization phenomena in single crystals, considering both higherorder energetic and dissipative effects. To this end, the Gurtin-type SGCP model pre-sented in Jebahi et al. [34] was extended to finite deformation following the procedure of Gurtin [21,23].
The capabilities of this model in regularizing idealized slip and kink bands were first revisited using a simple shear problem with a single slip system. The obtained results show the effectiveness of the model in capturing kink bands. However, this is not the case for idealized slip bands. Indeed, slip gradients in the direction normal to the slip plane do not correspond to GNDs.
Consequently, the present model behaves as a conventional crystal plasticity model with regard to idealized slip bands. To investigate the effectiveness of the present model in dealing with more general localization bands, it was applied to simulate a uniaxial tension problem with double slip systems symmetrically oriented with respect to the tension direction. For this problem, regularized localization bands slightly rotated with respect to the slip direction were obtained, showing the ability of the applied model to regularize general slip bands.
The good capabilities of the applied model in dealing with the latter bands motivated the application of this model to study the complex localization behavior with competition between slip and kink bands, which is encountered within single crystals undergoing single slip. A uniaxial tension problem of a single crystal plate having a single slip system was simulated using the present SGCP model. Different constitutive configurations, including purely higher-order energetic effects and purely higher-order dissipative effects, were considered. The associated results show that mesh-dependence difficulties, encountered using conventional crystal plasticity, can be removed by including higher-order effects, either energetic or dissipative.
Considering purely higher-order energetic effects, the results of the competition between slip and kink bands are qualitatively quite similar to those obtained by [29] using the micromorphic approach. At the beginning of straining, two thick kink bands appear and rapidly aggregate to form a larger kink region with concentrated lattice rotation at the center of the studied plate. The width and the distance between the formed kink bands depend on the energetic length scale l en . For sufficiently large l en , a single kink band at the center of the plate is formed at the beginning of straining. As straining increases, a central slip band appears and starts to interact with the kink region. With further straining, the formed slip band becomes increasingly intensified, leading to final necking of the plate.
Concerning higher-order dissipative effects, distinct results were obtained depending on the value of the initial higher-order dissipative slip resistance S ξ0 . The latter parameter can easily and independently be controlled in the present SGCP model due to the implied uncoupled dissipation assumption. Using zero S ξ0 (representing the progressive activation of the higher-order dissipative effects, as it is the case for the higher-order energetic ones), the obtained results are qualitatively quite similar to those obtained using purely higherorder energetic effects.
However, the use of a nonzero S ξ0 (as in classical Gurtin-type SGCP theories including higher-order dissipation) leads to quite different results. This choice leads to a localization behavior that is overwhelmingly dominated by slip banding, although analysis of lattice rotation in the studied plate shows the coexistence of kinking.

Conflicts of Interest:
The authors declare no conflict of interest.