Modeling of surface effects in crystalline materials within the framework of gradient crystal plasticity

A finite-deformation gradient crystal plasticity theory is developed, which takes into account the interaction between dislocations and surfaces. The model captures both energetic and dissipative effects for surfaces penetrable by dislocations. By taking advantage of the principle of virtual power, the surface microscopic boundary equations are obtained naturally. Surface equations govern surface yielding and hardening. A thin film under shear deformation serves as a benchmark problem for validation of the proposed model. It is found that both energetic and dissipative surface effects significantly affect the plastic behavior. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
In the last two to three decades, many effort s have been devoted to studying the mechanical properties of smallscaled and/or fine-grained crystalline materials widely used in small-scaled engineering such as microelectronics and microelectromechanical systems to ensure their performance and reliability in practical applications. It has been found in various experiments (see e.g., Fleck et al., 1994;Gu et al., 2012;McElhaney et al., 1998;Swadener et al., 2002 ) that crystalline materials at small scales usually possess plastic behaviors quite different from those of their bulk counterparts, among which is the size-dependence of the yield stress and strain hardening behavior. Since conventional plasticity theories are essentially size-independent, these size effects are beyond their predictability. As is well known, the macroscopic plastic deformation in crystals mainly results from the microscopic dislocation gliding on the individual crystallographic slip planes, thus those unusual plastic behaviors are closely relevant to the distinctive dislocation activities in small-scaled crystals. According to Ashby (1970) , dislocations in crystals can be classified into two types, namely statistically stored dislocations (SSDs) and geometrically necessary dislocations (GNDs). The latter is stored in crystals to accommodate lattice curvature. As discussed by Nix and Gao (1998) , GNDs may be responsible for size effects due to strain gradients since they are more likely to accumulate when inhomogeneous plastic deformation takes place. Furthermore, with the decrease of grain size and/or geometry size of crystals, the volume to grain boundary (GB)/surface ratio increases dramatically, which may prompt dislocations to reach and then interact with the GB/surface. For instance, as demonstrated by atomistic simulations (see e.g., Spearot et al., 2005;Van Swygenhoven et al., 2002;Yuasa et al., 2010 ), grain boundaries act as sources and sinks, emitting and absorbing dislocations. In single crystals, dislocations are observed to escape towards and nucleate at the surface, which dramatically changes the strain hardening behavior ( Oh et al., 2009;Shan et al., 2008;Wang et al., 2012 ). In order to construct predictive theories for the plastic behavior in small-scaled crystalline materials, the dislocation-relevant mechanisms both in the bulk and at the GB/surface need to be properly considered on the continuum level.
So far, since the pioneering work by Aifantis (1984) , numerous strain gradient plasticity theories have been proposed by incorporating the plastic strain gradient or the GND density, and have successfully predicted experimentally observed sizeeffects (see e.g., Chen and Wang, 20 0 0;Fleck et al., 1994 ). As discussed in Kuroda and Tvergaard (2008) , the existing strain gradient theories can be classified into two types namely the lower-order theories (see e.g., Acharya and Bassani, 20 0 0; Chen and Wang, 20 0 0; Evers et al.,20 02 ) and the higher-order theories (see e.g., Bargmann et al., 2011;Evers, 2004;Fleck and Hutchinson, 1993;1997;Groma et al., 2003;Gurtin, 20 0 0;20 02;Yefimov et al., 2004 ) based on how the effects of GNDs are considered. Of particular interest are higher-order theories, in which in addition to the traditional force balance equation and the traction boundary condition, the microscopic governing equations and the corresponding microscopic boundary conditions are involved. 1 The higher-order theories provide the possibility to model dislocation activities at a GB/surface through properly choosing microscopic boundary conditions. For simplicity, the microhard and microfree conditions representing two extreme cases are usually imposed, and the underlying assumption is that the GB/surface is completely impenetrable resp. freely penetrable to dislocations. However, these two extreme GB/surface models are not sufficient to capture dislocation activities at the GB/surface. To date, a few intermediate GB/surface models have been proposed. Cermelli and Gurtin (2002) developed a dissipative grain boundary model where a visco-plastic relation for grain boundary microtraction is assumed. By considering the interface energy contributed by the interaction between dislocations and grain boundaries, other researchers proposed energetic grain boundary models (see e.g., Aifantis et al., 2006;Fredriksson and Gudmundson, 2005;Gudmundson, 2004 ). These models are purely phenomenological. Gurtin (2008) constructed a grain boundary model where the effects of grain misorientation and grain boundary orientation are considered through the exactly derived grain boundary Burgers tensor. In van Beers et al. (2013) , a grain boundary model similar to that of Gurtin (2008) was put forward, in which a vectorial measure of the density of grain boundary defect was introduced. As illustrated by Gottschalk et al. (2016) , these two models are essentially identical if considering planar problems. van Beers et al. (2015) extended the previous model ( van Beers et al., 2013 ) by further accounting for the redistribution of the grain boundary defects. Ekh et al. (2011) suggested a microflexible boundary condition by which the crystallographic misorientation between adjacent grains is taken into account. The thermal effects on the model are further considered by Bargmann and Ekh (2013) . In the grain boundary model proposed by Wulfinghoff et al. (2013) , a criterion for interface yielding is introduced to account for the resistance of grain boundaries against slip occurrence. Kuroda (2017) proposed a strategy for modeling various microscopic interface boundary conditions by introducing a scalar quantity to control the slipping rate at the interface, based on which an orientation-dependent grain boundary model is proposed and then implemented numerically to study the plastic behavior of bicrystal micropillars.
Some models study dislocation-surface interactions in single crystals. Huang and Svendsen (2010) proposed an energetic surface model by considering the change of surface energy contributed by surface steps due to dislocation absorption by surfaces. Similarly, Hurtado and Ortiz (2012) derived an explicit expression for the density of surface energy by considering the relative orientation between slip systems and surfaces. In both models, the density of surface energy depends linearly on the plastic slips at the surface. By further considering the interaction energy between surface steps, Peng and Huang (2015) constructed a physically based energetic surface model within the framework of work-conjugated higher-order gradient crystal plasticity. In this model, a surface yielding condition is derived naturally, and interactions between slip systems are considered. Recently,  suggested a flexible surface boundary model by establishing a direct relation between the surface slip rate and the microforce vector projected onto the boundary. Using a surface yield criterion, a non-idealized surface behavior was addressed including dissipative effects of size-dependent surface yielding.
As reviewed above, modeling GB/surface effects has recently been the forefront in the field of crystal plasticity at small scales. Yet, in the existing GB/surface models, less attention has been paid on the underlying physical mechanisms. To incorporate the underlying physical mechanisms relevant to interactions between dislocations and surfaces as illustrated by in situ experiments into the continuum model, we revisit modeling dislocation absorption by surfaces. On the one hand, after dislocation absorption by surfaces, surface steps form (see e.g., Shan et al., 2008;Wang et al., 2012;Zhong et al., 2017 ), which results in the change of surface energy. Thus, such surface-driven change in the free energy must be properly considered. On the other hand, as is observed in Oh et al. (2009) , dislocations stop for a while near a surface contaminated by an oxide layer before they are absorbed. This is strong evidence for the existence of an increased resistance against dislocation absorption at surfaces affected by oxide layers, coatings or initial defects. Furthermore, according to the experiment by Zhong et al. (2017) , dislocation absorption by surfaces may be viewed as a dissipative process. All those indicate the dissipative nature of dislocation absorption by surfaces. In addition, models confined to small deformation are inadequate for situations with severe plastic deformation. Consequently, an advanced surface model with energetic and dissipative surface effects is proposed in the framework of finite deformation gradient crystal plasticity. Energetic surface effects are modeled by extending the physically motivated surface model proposed in Peng and Huang (2015) . To consider dissipative surface effects, a rate-dependent constitutive relation for the dissipative surface microforce is proposed within the new surface model. The theory is numerically implemented by using an in-house finite element code based on a dual-mixed finite element solution strategy to study a benchmark problem.

Basic kinematics of finite deformation crystal plasticity
In finite deformation crystal plasticity, the basic assumption is that the deformation gradient tensor F is multiplicatively decomposed into an elastic part F e and a plastic part F p as follows (1) The plastic part F p describes the plastic deformation from the reference configuration to the stress-free intermediate configuration resulting from dislocation gliding on the crystallographic planes. The elastic part F e by which the intermediate configuration is deformed to the current configuration comprises the lattice distortion and rotation. It is assumed that the plastic deformation does not change the volume such that det F p = 1 with det denoting the determinant of a tensor. In the following, if necessary, the subscripts i and c are used to identify quantities in the intermediate and the current configuration, respectively. Based on relation (1) , the velocity gradient tensor L is decomposed as with the elastic velocity gradient L e and the plastic velocity gradient L p defined as and where [ ∇a ] i j = a i, j denotes the gradient of a vector a , the superposed dot denotes the material-time derivative, the superscript −1 denotes the inverse of a tensor, and u is the displacement vector. Let us consider a crystal in which each slip system α is identified by a slip direction s ( α) and a slip plane normal m ( α) attached to the lattice space in the intermediate configuration.
The isoclinic assumption is adopted such that s ( α) and m ( α) remain unaltered from the reference configuration to the intermediate configuration.
The relevant strain measures, namely the Green-Lagrange strain tensor, and the right Cauchy-Green stretch tensor and their elastic parts are respectively defined as and where the superscript T denotes the transpose of a tensor and I is the identity tensor.
The evolution of edge and screw GND densities ρ ge (α) i and ρ gs (α) and ˙ ρ gs (α) where b is the magnitude of the Burgers vector, and and ρ gs (α) 0i being the initial GND densities.

Balance equations and boundary conditions
In the following, the macroscopic and the microscopic balance equations are derived via the principle of virtual power. Following Gurtin (2002) , in the intermediate configuration, the internal power in the bulk is assumed to be expended by the second Piola-Kirchhof stress tensor S e power-conjugate to the rate of the Green-Lagrange strain tensor ˙ E e , the microforce π (α) i power-conjugate to the slip rate˙ γ (α) and the microstress ξ (α) i power-conjugate to the slip rate gradient ∇ i ˙ γ (α) for each slip system α. In the present work, an additional contribution to the internal power expended by the microforce η (α) i power-conjugate to the slip rate ˙ γ (α) at the surface S i is introduced to consider the effect of dislocation absorption by surfaces. Thus, the total internal power P int is expressed in the intermediate configuration as where the second Piola-Kirchholf stress tensor S e is defined as S e = F e −1 · Jσ · F e −T with σ being the Cauchy stress tensor.
Accordingly, the external power P ext in the absence of body forces is written as where T i is the prescribed traction and (α) i is the microtraction defined in the intermediate configuration.
According to the principle of virtual power, the variation of the internal power with respect to the velocity ˙ u and the slip rate ˙ γ (α) equals that of the external power, which gives From Eqs. (2) , (3), (5) and (7) , the variation of ˙ E e is expressed as Then, by substituting Eq. (13) into Eq. (12) and using the divergence theorem, the principle of virtual power is rewritten as where N i is the surface outward normal vector, Div denotes the divergence operator, and the Mandel stress tensor M is expressed as First, considering the special case without plastic slip ( δ˙ γ (α) = 0 ), the requirement that Eq. (14) should be satisfied for arbitrary δ˙ u results in the macroscopic force balance equation and the corresponding standard traction condition in the intermediate configuration, i.e., and at the part of surface S T i with prescribed traction. Then, given the validity of Eq. (14) for arbitrary δ˙ γ (α) in the case without macroscopic displacement ( δ˙ u = 0 ), the microscopic balance equation and the associated microscopic boundary condition for each slip system α are obtained as follows being the Schmid stress, and

Constitutive relations: theory with surface effects
The elastic behavior of the crystal is assumed to be compressible neo-Hookean such that the strain energy density ψ e i has the following form with J = det ( 2 E e + I ) , and λ and μ being the Lamé constants. Then, the hyperelastic constitutive relation for S e is derived as Generally, the microstress ξ (α) i may comprise an energetic part and a dissipative part ( Bargmann et al., 2014;Gurtin and Anand, 2005 ). Here, for simplicity, ξ (α) i is assumed to be purely energetic. According to Kuroda and Tvergaard (2008) and Ertürk et al. (2009) , the energetic microstress ξ (α) i in the work-conjugate theories can be essentially related to the back in the non-work-conjugate theories (see e.g., Bayley et al., 2006 ) as follows In Ertürk et al. (2009) , based on relation (22) and the expression of the back stress τ b (α) i derived by Bayley et al. (2006) through considering the elastic interactions between GNDs within the same slip system and those from different slip systems, the constitutive relation for ξ (α) i is given. Here, we directly adopt the relation therein, thus ξ (α) i is expressed as where ν is Poisson's ratio and l denotes the radius of the domain within which the interaction between GNDs is considered. Alternatively, the interaction between different slip systems can be accounted for in the relation for GND densities in Eqs. (8) and (9) as done in Bargmann et al. (2011) and  , but this way is not pursued here.
The scalar microforce π (α) i being dissipative in nature is assumed to follow a visco-plastic power-law, i.e., where sgn() denotes the sign function, m b is the rate-sensitivity exponent in the bulk, and ˙ γ b 0 is the reference slip rate. The slip resistance R b (α) i resulting from the random trapping processes between SSDs has the following linear hardening form, where γ (β ) The microforce η (α) i at the surface is divided into an energetic part η en (α) i and a dissipative part η dis (α) i . According to Peng and Huang (2015) , the power expended by the energetic microforce η en (α) i at the surface equals the change of the surface energy resulting from surface steps formed after dislocation absorption, where ˙ ψ S i is the rate of the surface energy density resulting from surface steps formed after dislocation absorption. In Peng and Huang (2015) , the expression for the density of surface energy in the small deformation case was strictly derived by considering both the surface energy change due to the increase of surface area and the interaction energy between surface steps. Following the idea therein and considering the effect of finite deformations, the expression for ˙ ψ S i can be readily derived, and without going into the details, the result is directly given here, i.e., The first and second terms represent contributions from the increase of surface area and the interaction between surface steps, respectively. (α) 1 denotes the surface free energy per unit surface step area and (αβ ) 2 is the surface hardening modulus measuring the interaction strength between surface steps from slip systems α and β. In Eq. (29) , the surface outward normal vector N i in the intermediate configuration depends on deformation, and, hence, the explicit expression for the density of surface energy is not available, which is different from the small deformation case where the surface outward normal vector remains unaltered. Then, from Eq. (28) , η en (α) i is expressed as The dissipative surface microforce η dis (α) i is introduced to consider the dissipative nature of dislocation absorption by surfaces as indicated by in situ experiments ( Oh et al., 2009;Zhong et al., 2017 ). To construct a constitutive relation for η dis (α) i physically, detailed information on the process for a dislocation absorbed by a surface is necessary, which nevertheless is not yet available in the literature. As an estimate, the resistance force against dislocation absorption is proportional to the rate of dislocation absorption (measured by the plastic slip rate at the surface). Moreover, the resistance force increases with the accumulated density of surface steps. These can be accounted for in the following rate-dependent power-law relation ˙ γ s 0 is the reference slip rate at the surface and m s is the surface rate-sensitivity exponent. The surface slip resistance R s (α) i is expressed as where γ s (β ) acc = ˙ γ (β ) d t is the accumulated slip at the surface, R s (α) 0 denotes the initial surface slip resistance, and H s( αβ) is the dissipative surface hardening modulus. For imperfect surfaces with surface coatings, oxide layers or initial defects, the surface parameters depend on the initial surface state. It is worth mentioning that in some grain boundary models ( van Beers et al., 2013 ), a similar form for dissipative grain boundary microforce is adopted, but the slip resistance is simply assumed to be constant, by which the dissipative hardening is ignored.
Neglecting the microtraction (α) i , the microscopic boundary condition Eq. (19) is rewritten as which constitutes the governing equation for dislocation absorption by surfaces. Both energetic and dissipative surface effects are accounted for.

Finite element implementation
For the rate-dependent initial boundary value problem from the present rate-dependent model to be solved for arbitrary geometries and boundary conditions, a dual-mixed finite element formulation ( Ekh et al., 2007 ) is adopted. The force balance Eq. (16) and the evolution Eqs. (8) and (9) for the GND densities are chosen as governing equations, and, hence, the displacement u and the GND densities ρ ge( α) and ρ gs( α) are considered as primary variables. The plastic slip γ ( α) is evaluated locally at the integration points. The treatment for the balance of momentum (16) is straightforward, but due to the surface effects, the dicretization and implementation of the GND evolution equations deserve some attention. Since the evolution equations for the two kinds of dislocations are similar, it is only explained for edge GND densities for the sake of brevity. To obtain the weak forms of the governing equations, they are multiplied by weighting functions δu and δ˙ ρ ge (α) , which vanish at the part of the surface with natural boundary conditions, and then integrated over the volume.
Next, applying the divergence theorem, we obtain where the surface traction boundary condition (17) has been substituted into Eq. (34) . The plastic slip rate ˙ γ (α) at the surface in Eq. (35) is evaluated locally via the microscopic surface boundary condition (33) , and the initial GND densities are neglected (but are straightforward to incorporate). For the spatial discretization, the volume of the body is divided into finite elements in each of which, following the standard Galerkin approach, the unknown fields of the displacement and GND densities and the weighting functions are approximated by their nodal values multiplied by shape functions. For time integration within a time increment [ t n , t n +1 ] , the implicit backward Euler scheme is adopted such that ˙ γ (α) where γ (α) = γ (α) n +1 − γ (α) n , ρ ge (α) = ρ ge (α) n +1 − ρ ge (α) n , t = t n +1 − t n and the subscript n identifies the quantities at the n th time step. Then, the weak forms (34) and (35) reduce to a highly nonlinear, strongly coupled system of equations at each time step which are solved by the Newton-Raphson method.
For the incorporation of the surface model into the finite element formulation, the evaluation of slip rate at the surface involved in the surface integration in Eq. (35) is required. To this end, the increment of plastic slip γ ( α) at the surface is determined by the microscopic boundary condition (33) as where plastic slips required to evaluate the terms at the right hand side are extrapolated from the bulk. Particularly, if the surface dissipative microforce is ignored, the increment of plastic slip γ ( α) at the surface is evaluated by which constitutes the microscopic boundary conditions for a purely energetic surface.

Numerical example
A plane strain benchmark problem of plane constrained shear of a single crystalline thin film is studied. As illustrated in Fig. 1 , the film with thickness h in the x 2 -direction is assumed to be infinitely long in the x 1 -direction, and without loss of generality, two active slip systems with each of them defined by a slip direction s (α) = [ cos θ α , sin θ α ] and a slip plane normal m (α) = [ −sin θ α , cos θ α ] are considered. For the present plane strain problem, only edge GNDs are involved. The thin film suffers a homogeneous external shear rate ˙ γ ext with its lower surface being constrained such that the macroscopic boundary conditions read where u 1 and u 2 denote the displacement increments after time increment t . Microscopically, the single crystal thin film is free-standing, and therefore the governing Eq. (33) for dislocation absorption by surfaces acts as the microscopic surface boundary condition. In addition, to model the infinite extension of the film in the x 1 -direction, a representative part of the thin film with length L is considered and periodic boundary conditions are imposed for both the displacement and the GND densities where the periodic length L can be arbitrarily chosen without affecting the solution. For the spatial discretization, 4-node bilinear quadrilateral elements with 2 × 2 full integration are chosen for both the displacements and the GND densities. To numerically evaluate the surface integral, two additional integration points are placed on the surface edges of elements located at the surface. The

Validation of the FEM solution
First, to validate the finite element implementation summarized in Section 3 , the FEM solution is compared to the analytical solution of ( Peng and Huang, 2017 ). The analytical solution in Peng and Huang (2017) only applies to the small deformation case without dissipative effects. However, if the applied load is small and dissipative microforces both in the bulk and the surface are ignored, the solution of the present problem should approach the analytical one. Thus, in the following comparison, the dissipative forces π (α) i and η dis (α) i are temporarily ignored. To avoid repetition, the results from the analytical solution are directly given here without presenting the detailed formulation which can be found in Peng and Huang (2017) . For both the numerical and analytical cases, the associated parameters other than those given above are taken as h = 1 μm , θ 1 = π / 12 , θ 2 = 5 π / 12 , 1 = 10 N / m and 2 = 100 N / m . In Fig. 2 a, the shear stress σ 12 , which is homogeneous along the x 2 -direction, is plotted as a function of the externally applied shear γ ext . The numerical results correspond very well to the analytical results for γ ext < 0.005. If γ ext increases, a slight difference is noticeable which is attributed to the finite deformation framework applied here. The stress-strain curves are divided into three stages by two yielding points which denote the onset of dislocation absorption by surfaces for the two slip systems. As dissipative resistances are ignored at this stage, plastic deformation occurs immediately after loading, and, hence, there is no bulk yielding point on the stress-strain curve. In Fig. 2 b and c, the distributions of plastic slips and GND densities along the x 2 -direction are shown, where three values of γ ext namely 0.0 01, 0.0 02 and 0.005 representing the three stages respectively are taken. For both plastic slips and GND densities, the FEM curves coincide almost perfectly with the analytical curves, indicating the excellent match between the FEM and analytical solutions.

Influence of energetic and dissipative surface effects
In the present surface model, both energetic and dissipative surface effects are involved. Here, to illustrate their influence separately, two specific cases, i.e., the case with only energetic surface effects and the one with only dissipative surface effects are considered, respectively. In both cases, θ 1 = π / 3 and θ 2 = 2 π / 3 are taken.
First, we consider the case with only energetic surface effects. For different values of 1 and 2 , the stress-strain curve, the evolution of the total average GND density ρ ge = h 0 ρ ge (1) + ρ ge (2) d x 2 /h, the evolution of slip rate ˙ γ s at the surface, and the distribution of plastic slip γ for γ ext = 0 . 01 are plotted in Fig. 3 a-d, respectively. Since the results for the two slip systems are similar due to symmetry, only those for slip system 1 are displayed. In Fig. 3 a, the stress-strain curves for the energetic surface model are generally divided into three stages by two points respectively denoting macroscopic bulk yielding and the onset of dislocation absorption by surfaces (or alternately surface yielding). Similar stress-strain curves with both bulk and surface/interface yielding points are also predicted by some other models with surface/interface effects (see e.g., Aifantis et al., 2006;van Beers et al., 2013;Fredriksson and Gudmundson, 2005;. After bulk yielding, the surface initially remains impenetrable due to the lack of driving force for dislocation absorption, and the hardening rate is the same as that predicted by microhard boundary condition. If the externally applied shear reaches a critical value, the driving force is sufficient to overcome the energetic surface resistance for dislocation absorption, surface yielding occurs, resulting in a considerable decrease in the hardening rate. From Fig. 3 b, after surface yielding due to the absorption of dislocations by the surface, the accumulation rate of the total average GND density decreases in accordance with the decrease of hardening rate. Comparing the results for different values of 1 and 2 reveals that the critical point for surface yielding only depends on the surface energy parameter 1 . A larger 1 results in surface yielding at a larger γ ext . The hardening rate and the accumulation rate of GND density after surface yielding are governed by the energetic hardening modulus 2 in a way that a larger 2 yields a larger hardening rate. As shown in Fig. 3 c, after surface yielding, the slip rate at the surface (measuring the rate of dislocation absorption by surface) increases immediately from zero to a constant level depending on the chosen magnitude of 2 . Moreover, the surface resistance against dislocation absorption increases with increasing 2 . From Fig. 3 d, for smaller 1 and 2 , plastic slips are easier to develop and distribute more homogeneously. The curves predicted by the energetic surface model are bounded by those for microhard and microfree models, and the energetic surface model can readily reduce to the microhard case or the microfree case by increasing the surface parameters 1 and 2 to infinity or decreasing them to zero. From Fig. 3 a, strain hardening after surface yielding is insignificant unless the energetic hardening modulus is impractically large. 2 A similar conclusion was drawn for the energetic surface contribution in the model of Hurtado and Ortiz (2012) . There, an impractically large surface parameter was required in order for the model to capture size-dependent strain hardening of micropillars under compression. The current results imply that the energetic surface effects are not sufficient to account for the size-dependent hardening in small-sized single crystals. For that reason, the role of dissipative surface effects needs to be investigated. Thus, we consider next the case with only dissipative surface effects. To this end, the energetic surface microforces are ignored by setting 1 and 2 to zero. The initial surface slip resistance R s 0 and the dissipative surface hardening modulus H s are varied. In Fig. 4 a-d, the stress-strain curve, the evolution of the total average GND density ρ ge , the evolution of slip rate ˙ γ s at the surface, and the distribution of plastic slip γ for γ ext = 0 . 01 for different values of R s 0 and H s are displayed, respectively. The purely dissipative surface model predicts plastic behavior quite similar to that of the purely energetic one. Particularly, the surface yielding depends on the initial slip resistance R s 0 , while the hardening rate, the accumulation rate of the total GND density after surface yielding, and the slip rate at the surface at steady-state are governed by the dissipative surface hardening modulus H s .
In view of reversibility associated with surface step formation, we explore the influence of surface effects on the plastic behavior during unloading. To this end, the stress-strain curve, the evolution of slip rate at the surface during a loading and unloading cycle predicted by the purely energetic surface model, the purely dissipative surface model and the general surface model with both energetic and dissipative surface effects are displayed in Fig. 5 a and b, respectively, where 1 = 10 N / m , 2 = 100 N / m , R s 0 = 5 N / m and H s = 200 N / m . The loading is applied up to γ ext = 0 . 01 , followed by full unloading. From the figures, the curves predicted by the two models are similar during loading. For unloading, the situation is different. For the purely energetic surface model, during unloading, the energetic surface microforce, which can be regarded as a surface back-stress, acts as the driving force for the recovery of surface steps (inverse to dislocation absorption during loading). Hence, bulk yielding and surface yielding occur at the same time. The hardening rate and the slip rate at the surface immediately after yielding (the second unloading stage) are identical to that of the third loading stage with both bulk and surface hardening. If γ ext decreases to a critical value, the energetic microforces decrease to zero, the surface becomes impenetrable again. After that the hardening rate remains the same as that of the second loading stage with only bulk hardening, and the slip rate at the surface decreases to zero. For the purely dissipative surface model, during unloading, due to the existence of the dissipative surface microforce, the recovery of surface steps does not take place until the driving (d) distribution of plastic slip for γ ext = 0 . 01 . Surface yielding depends on initial surface resistance R s 0 . Surface hardening, the accumulation rate of the total average GND density after surface yielding, and the slip rate at the surface at the steady-state are governed by dissipative surface hardening modulus H s . For smaller R s 0 and H s , plastic slips develop faster and distribute more homogeneously after surface yielding. The results of the dissipative surface are bounded by those of the microhard and microfree models.
force (the net contribution of bulk and surface energetic microforces) is sufficient to overcome the dissipative resistance. Hence, the bulk yields first and, subsequently, the surface yielding occurs only if the driving force overcomes the resistance. The hardening rate and the slip rate at the surface for the second and third unloading stages are identical to those for the corresponding loading stages. For the general surface model with both dissipative and energetic surface effects, the loading and unloading behaviors are similar to those of the purely dissipative surface model. In addition, a pronounced Bauschinger effect is also observed in Fig. 5 a.

Influences of film thickness and orientation of slip systems
To explore the influences of the film thickness h , for different values of h , the stress-strain curve, the evolution of the total average GND density ρ ge , the evolution of slip rate ˙ γ at the surface, and the distribution of plastic slip γ for γ ext = 0 . 01 are plotted in Fig. 6 a-d, where θ 1 = π / 3 and θ 2 = 2 π / 3 , 1 = 5 N / m , 2 = 100 N / m , R s 0 = 5 N / m and H s = 200 N / m . The results displayed in Fig. 6 clearly depict a size-dependent material response. In thinner films, the flow stress and the total average GND density at a specific γ ext are larger, which indicates that "smaller is stronger". As shown in Fig. 6 a and c, in thinner films, surface yielding occurs at a smaller applied shear γ ext , and the slip rate at the surface (measuring the rate of dislocation absorption) at the steady state is larger, which implies that surface yielding occurs more easily at smaller scales. In fact, in the experiments by Gruber et al. (2008) , it is found that for single-crystalline Au thin films with smaller thickness, after bulk yielding, the strain hardening rates decrease suddenly at a specific point which is suspected to be the surface yielding point, while for films with larger thickness, no such phenomenon occurs, which is qualitatively consistent with the predicted results here.   The size effect of plastic slip close to the surface is inverse to that in the middle of the film, see Fig. 6 d. It is attributed to the difference between governing mechanisms for slip development in the middle and near the surface. In the middle, plastic deformation is more difficult to develop in thinner films, resulting in a smaller plastic slip therein. However, near the surface, the slip development is governed by dislocation absorption by surfaces. As in the thinner film, the onset of dislocation absorption occurs more easily, the slip near the surface is larger.
To investigate the influences of the orientation of slip systems, for the case with two symmetric slip systems ( θ 1 = π / 3 , θ 2 = 2 π / 3 ) and the case with two asymmetric slip systems ( θ 1 = 5 π / 18 , θ 2 = 11 π / 18 ), the stress-strain curve, the evolution of the total average GND density ρ ge , the evolution of slip rate ˙ γ at the surface, and the distribution of plastic slip γ for γ ext = 0 . 01 are shown in Fig. 7 a-d. The bulk/surface yielding and the strain hardening rate significantly depend on the orientation of slip systems, see Fig. 7 a and c. For the symmetric case, only slight differences between the two systems exist. In the asymmetric case, the slip systems possess different Schmid factors. Thus, the onset of plastic slip in the bulk and at the surface for the two slip systems occur at different applied shears γ ext , which results in several "yielding points" on the stress-strain curve in Fig. 7 a. As illustrated in Fig. 7 c, the slip rate at the surface of slip system 2 ( θ 2 = 11 π / 18 ) suffers a sudden change after the onset of plastic slip at the surface for slip system 1 ( θ 1 = 5 π / 18 ), which indicates that the interaction between slip systems significantly influences dislocation absorption by surfaces.

Summary
In this contribution, a surface model has been developed within the framework of finite deformation gradient single crystal plasticity theories. Energetic as well as dissipative surface effects are considered via an energetic and a dissipative surface microforce for each slip system. The physically motivated energetic surface microforces result from the change of surface energy due to the formation of surface steps after dislocation absorption by surfaces. Motivated by experimental findings ( Oh et al., 2009;Zhong et al., 2017 ), the dissipative surface microforces are introduced to account for the dissipative nature of dislocation absorption by surfaces. The key feature of the surface model is an evolution equation for the slip rate at the surface. It constitutes the governing equations for dislocation absorption by surfaces.
To realize the numerical implementation of the new model, a dual-mixed FEM is extended by incorporating the new surface model. As an application, the benchmark problem of plane constrained shear of a single crystal thin film is studied. The validation of the FEM implementation is verified by comparing the numerical results with analytical solutions. Energetic and dissipative surface effects dramatically affect the plastic behavior of thin films. The plastic behaviors predicted by the purely energetic surface model and the purely dissipative model are similar during loading but differ considerably during unloading. The energetic surface model predicts insignificant strain hardening after surface yielding unless the surface parameter is impractically large, implying the importance of the dissipative surface effect. The influence of film thickness and the orientation of slip systems on surface yielding and hardening is also found to be significant. Surface yielding is more likely to occur at smaller scales. The predicted film thickness-dependence in terms of the stress-strain response is qualitatively consistent with the experimental results of Gruber et al. (2008) .