Interplanar potential for tension-shear coupling at grain boundaries derived from ab initio calculations

Based on ab initio density functional theory (DFT) calculations we derive an analytical expression for the interplanar potential of grain boundaries and single crystals as a function of coupled tensile and shear displacements. This energy function captures even details of the grain boundary behaviour, such as the tension-softening of the shear instability of aluminium grain boundaries, with good accuracy. The good agreement between the analytical model and the DFT calculations is achieved by introducing two new characteristic parameters, namely the position of the generalised unstable stacking fault with respect to the stable stacking fault, and the ratio of stable and unstable generalised stacking fault energies. One of the potentials’ parameters also serves as a criterion to judge if a grain boundary deforms via crack propagation or dislocation nucleation. We suggest this potential function for application in continuum models, where constitutive relationships for grain boundaries need to be derived from a sound physical model.


Introduction
For many modern nanostructured materials the mechanical characteristics of interfaces dominate deformation and failure processes. Hence, mechanical models of interfaces based on the underlying physics are desirable to understand and to quantify the relation between interface characteristic properties and macroscopic material behavior. Establishing such models is challenging, because the deformation mechanisms at interfaces and even scalar properties such as the tensile and shear strength do only depend to a certain extent on the macroscopic degrees of freedom (such as grain boundary misorientation) but predominantly on the details of the atomistic structure [1,2]. These details in turn depend on the nature of the chemical bond between the atoms at the interface [3][4][5]. Nevertheless, a physical quantity which can connect bonding and grain boundary structure on the one hand, and bonding and mechanical properties on the other hand, is the stacking fault energy (SFE) of the preferred dislocation slip system, in fcc metals { }⟨ ⟩ 1 1 1 1 1 2 . This can be made plausible by the fact that grain boundary structures can be modelled as arrays of dislocations [6,7], and that the structure of dislocation cores varies with the stacking fault energy as well. In cases of very low stacking fault energy it is even directly visible from the relaxed structure of the grain boundary. One of the clearest examples is the formation of the 9R structure at the Σ3 (1 1 2)[1 1 0] tilt grain boundary in low-SFE materials, such as Cu and Ag [8,9].
Thus the SFE can provide the starting point for an energy or potential based mechanical model of grain boundary sliding. Together with our previous insight that the energy as a function of tensile displacement at grain boundaries of different structure [10] and chemical composition [11] can be scaled in the spirit of the Rose-Ferrante-Smith universal binding energy relationship (UBER) [12] this enables us to set up a coupled interplanar potential and to derive constitutive relations for combined tensile and shear loads at grain boundaries. Such potential based models, like the Peierls concept for dislocations [13], have a long tradition in materials science and are very successful despite or maybe because of their apparent simplicity.
The Griffith' model of fracture relies on the surface energy γ s as an important physical parameter to predict cleavage or crack blunting as well as to predict the fracture stress in brittle materials [14]. A dislocation emitted from a crack tip leads to blunting of the crack tip and a reduction in the crack-driving force by elastic shielding, which increases the ductility of a material [15]. Rice proposed an analysis of dislocation emission from a crack tip [14] based on the Peierls concept. The characteristic parameter in this analysis is the unstable stacking fault energy γ us , which is the energy barrier describing the sliding of blocks in the direction of the Burgers vector → b, where → b corresponds to a full dislocation or a Shockley partial. For mode II loading, the slip plane coincides with the crack plane, and from Rice's solution for mode II loading the critical stress intensity factor for the emission of a dislocation can be derived as a function of shear modulus and Poisson's ratio [14], again two well defined single crystal parameters. However, Rice's original work neglected the tensionshear coupling effect, which can alter the interplanar potential, and the formation of a ledge which can blunt the crack tip [16][17][18].
For the tension-shear coupling, Sun et al [19] first introduced a convenient analytical representation based on Rice's model. Their atomistic simulations of slip in single crystal solids showed that the normal stress distribution across a slip plane will reduce the critical loading required for dislocation emission from a crack tip [19], which is called the tension-softening effect. Lazar et al [20] modified their coupled potential to describe the nonsinusoidal energy profiles for the { }⟨ ⟩ 111 1 1 0 slip system in NiAl, and also reproduced the tension-shear coupling effect very well. However, the energy profiles for slip at grain boundaries are much more complex than those for single crystals and are not captured yet by this model.
Several ab initio studies also focused on the influence of superimposed normal stress on the theoretical shear strength in single crystals (bcc, fcc metals and covalent crystals) [21][22][23]. The results show that the theoretical shear strength follows a linear relationship with the normal stress in metals. The theoretical shear strength increases with increasing compressive normal stress and decreases with increasing tensile normal stress. However, in covalent crystals the relations are different in different crystal structures because of the directional characteristics of covalent bonding. Notwithstanding the importance of interfaces in polycrystals, little attention was paid to a coupling of various stress tensor components at grain boundaries so far. In our study we investigate the coupling of shear and tensile loads on aluminium single crystal and grain boundaries to fill this gap.
In our previous work [10,24] we tested the universal behavior of tensile and shear processes separately at grain boundaries. In case of tensile loading [10], we found that all the different investigated geometries exhibit energy-displacement curves that can be brought into coincidence with the UBER, as mentioned above. For the shear case [24], there is no simple relationship to describe the theoretical shear strength and sliding displacement, but for periodic grain boundary structures, the energy-displacement curves of a shear processes (the generalised stacking fault energy or γ-surface) can be expressed as an expansion in sinusoidal functions [25]. The objective of the paper at hand is to focus on the tensionshear coupling modes at grain boundaries and to propose an extended potential-based, threedimensional constitutive relationship for the interface plane that can be used to simulate deformation and fracture at grain boundaries as well as the prediction of dislocation nucleation at the interface.

Computational details
To investigate the grain boundary energy under tensile and shear strain we carried out ab initio total energies and electronic structure calculations using the Vienna ab initio simulation package (VASP) [26,27]. The exchange-correlation effects were treated in the local density approximation (LDA) by using the exchange-correlation of Ceperley and Alder [28] as parametrized by Perdew and Zunger [29]. Electron-ion interactions were treated within the ultrasoft pseudopotential (USPP) method [30,31], with the valence configuration (3s2 3p1) for Al. The plane-wave cutoff was chosen to be 450 eV, which corresponds to the value we obtained from our previous convergence tests [10,24]. In that work we also obtained the LDA lattice constant of Al as 3.977 Å, and describe the generation and optimisation of the atomistic models of various grain boundary structures used in the study at hand.
Starting from the fully optimised grain boundary structures, we applied normal and tangential rigid displacements, ∆ n and ∆ t respectively, between the grain boundary and the next crystallographic plane to obtain the energy as a function of a three-dimensional displacement. In a second step, for fixed lateral displacements all atomic positions were relaxed perpendicular to the grain boundary plane.
We use the term cleavage energy G c to denote the work of brittle fracture, which is needed to separate the fully relaxed grain boundary structure into two free surfaces in a rigid manner, i.e. without subsequent relaxation of the atomic positions. The energy difference between the grain boundary and two fully relaxed free surfaces is the work of separation W sep . It was calculated for all grain boundaries presented here in [10,24] and is reprinted for convenience in table 1, but it is not used for the construction of the coupled potential. For simplicity we use the terms unstable and stable stacking fault energy, γ us and γ sf , not only for the single crystal, but also when discussing the energy landscape for grain boundary sliding. In the latter case, the reference energy is the energy of the grain boundary in its equilibrium state, which would have to be added to obtain the total defect energy. However, to investigate the resistance with respect to a shear displacement, only the energy difference is of interest.

Ab initio Results
One slip system in the single crystal and several in special symmetrical tilt grain boundaries (STGBs) as well as one twist boundary were selected: the [1 1  We now present some ab initio results on the energy surface ( ) Ψ ∆ ∆ , t n for sliding and opening.  (1 1 1) planes. Due to the rotational symmetry of the γ-surface for these planes, no additional vector is needed to cover sliding in these slip planes [24]. Figure 1(a) shows the energy ( ) Ψ ∆ ∆ , t n as a function of the slip displacement ∆ t at various opening displacements ∆ n for the (1/6)[1 1 2 ](1 1 1) slip in fcc Al. Symbols represent data obtained from ab initio calculations. This slip process exhibits a nonsinusoidal energy profile; the maximum energy γ us occurs at , and the local minimum at t p corresponds to the stacking fault energy. In figure 1(a) we see that with finite and positive opening displacements the energy at zero slip displacement increases, but the energy barrier and the energy difference between unstable stacking fault energy and stable stacking fault energy decreases, which is the tension-softening effect. Similarly, the offset at zero slip displacement increases with a negative displacement, however, in this case the energy barrier increases. It can be noticed that the position of the unstable stacking fault energy slightly shifts towards larger slip displacements as the opening displacement increases. The ratio of the positions of the unstable and stable stacking faults as a function of opening displacement is shown in figure 2. Figure 1(b) shows the energy ( ) Ψ ∆ ∆ , t n as a function of separation ∆ n with the slip displacement as a parameter. It is clearly seen that the minimum energy point shifts to larger normal separation with increasing slip displacements, a process that lowers the cleavage energy. All curves coincide at larger normal separation at the value of two times the surface energy: 2γ s .
In the Σ3 STGB, the direction perpendicular to the easy sliding direction [1 1 0] in the grain boundary plane is the [1 ¯1 1] direction. Starting from the absolute minimum of the γ surface, sliding along this direction initially would lead through an elongated valley of the γ-surface but then encounter the maximum peak of this energy landscape [24]. So we assume that the system would rather deviate to a parallel path along [1 ¯1 1] and study the energy barrier of this parallel path to check how the tension affects shear properties. The corresponding stacking fault is located at [¯¯]a 0.35 1 1 1 0 , thus we define the partial displacement as The ab initio calculation results for the energy ( ) Ψ ∆ ∆ , t n as a function of the slip displacement ∆ t at various opening displacements ∆ n are displayed as symbols in figure 3. In this case, the position of the energy barrier is not located at the middle point of the partial Burgers vector, however it remains fixed, i.e. it is independent of the normal displacement. Compared to the previous case, the energies at zero slip displacement increase only slightly with normal displacement, but the unstable stacking fault energies and stable stacking fault energies, taken w.r.t. to the zero-point energy, decrease significantly.

Development of tension-shear coupling model
From the ab initio results the demands on the required potential become clear. They are visualised in a set of schematic diagrams of energy versus displacement (opening and sliding) in figure 5. Energy versus displacement graphs for tensile displacements are placed on the left side and for shear displacements are placed on the right side. Figure 5(b) illustrates that for a given tensile displacement ∆ * n , the energy landscape for shear is expected to become more flat (the tension-softening effect). For a given shear displacement ∆ * t , the minimum of the tensile E-∆ n curve is expected to move to higher energies and larger displacements, as shown schematically in figure 5(c).
To derive a potential which can describe these coupling effects, we start from the uncoupled stresses that occur in a single crystal upon either shear or normal displacement. The energydisplacement curve for a purely tensile loading is assumed to be given by a Rose-type UBER, where l is the interface characteristic length corresponding to the position of the inflection point in the tensile-displacement curve. The initial assumption of the Peierls model for the response to a shear displacement of two half-crystals is that it can be described by a simple sinusoidal function,  To take this into account we use a combination of quadratic sine functions for our potential (see below) and introduce two new parameters: the slip ratio of the unstable stacking fault to the Burgers vector, s, and the ratio of the stacking fault energy γ sf to the unstable stacking fault energy γ us , K: x us is the slip displacement of the (generalised) unstable stacking fault, and | | → b is the length of the full/partial Burgers vector.
To describe the influence of normal (tangential) displacements on the shear (tensile) stress we introduce coupling terms in the expressions for the pure tensile or shear stresses which where x 0 is the opening displacement corresponding to zero tensile stress for the shear displacement x us , i.e. it is the position of the minimum of the UBER-like tensile curve for ∆ = x t u s . With these and our additional parameters, plus the demand that the mixed partial derivatives be equal, we construct the common potential:     The empirical term '1 − s' in equation (5) is used to enhance the tension-softening effects. Equation (6) takes care of the abovementioned fact that the energy profile for sliding is not always simple sinusoidal, by using a combination of sinusoidal functions. Furthermore this expression enables us to introduce a phase shift under normal loading. If the initial position of the unstable stacking fault is not located at | | → b /2, this position is likely to shift either to larger or to smaller displacements if a normal loading is applied. For instance in single crystalline Al the ratio s in figure 2 increases with increasing normal displacement. This means that the energy barrier for sliding along the [1 1 2 ] direction is shifting to a larger tangential displacement for increasing normal displacements. The sign of the shift can also be negative, depending on the ratio: K/s. To reproduce this effect, we introduce an exponential term to modify the argument of the sinusoidal function in the correct direction, depending on whether K/s < or K/s > 1. For K = 0 and s = 1/2 expression (6) is reduced to the original formulation of Sun et al.
The analytical potential (5) has two important benefits: firstly, it can describe energy profiles for sliding beyond simple sinusoidal profiles for ≠ s 1/2. Secondly, athough the analytical expression might look cumbersome, it only requires a few input parameters which have physical meanings and can be calculated using ab initio methods: the generalised unstable stacking fault energy γ us and the stable stacking fault energy γ sf of the grain boundary for no opening displacement, ∆ = 0 n . The scaled length l can be obtained from the tensile energydisplacement curve for ∆ = 0 t , x 0 for ∆ = ∆ t u s , the position of the energy barrier γ us . The characteristic parameters of the tension-shear coupling obtained from our ab initio results are summarized in table 2. Sun et al already described the single crystal slip system with their original analytical formulation [19]. Their stacking fault energy calculated with an embedded atom method (EAM) potential, 14.1 mJ m −2 , is an order of magnitude lower than the experimental values from Hirth and Lothe [32], 110-280 mJ m −2 , and DFT results [24,33], 128 mJ m −2 and 124 mJ m −2 . Nevertheless, their q-value of 0.085 is close to our ab initio result of 0.101. However, the potential of Sun et al imposes that γ us occurs at ∆ = | | → b 0.5 t p , which influences the tension-softening behaviour and can't describe the sliding curves accurately (see dashed line in figures 1(a)). To verify the practicality of our extended potential, the modeling curves obtained from the potential are plotted to compare with the ab initio calculation results.

Case I: stacking fault energy γ ≠ 0 sf
The analytic expressions for our potential are shown as solid curves in figure 1, together with the ab initio data. The coupled potential describes the calculated results accurately, including the position and shift of the energy barriers, governed by the ratio s, the barrier heights, and the shift of the minimum energy points under normal loading. With increasing opening displacement, the heights of energy barriers are lowered. This tension softening effect is related to the parameters s and q. Note that also the compressive regime For the case of the [1 ¯1 1] direction in the Σ3 STGB, the energy barrier is not located at the midpoint of the partial dislocation, but at sf us approximately equals s, so the exponential term in equation (6) is close to 1: . Thus equation (6) was simplified as: ( ) displacement ∆ t at various opening displacements ∆ n for slip along the direction of [1 ¯1 1] of the Σ3 STGB is shown in figure 3. Again, the solid curves are obtained from our potential and reproduce the ab initio data very well.  (6) is p . For these cases the potential is simplified to Sun's potential apart from the new factor '1 − s' which enhances the tension-softening effect. As visible in figure 4 this effect is reproduced extremely well by our potential.

Predicted constitutive behaviour
Using this potential, we can further analyse the mechanical properties of grain boundaries. The tensile stress is the first derivative of the potential with respect to normal displacement: , , n t n t n n (7) and the maximum tensile strength, or fracture strength, is obtained as the slope in the inflection point: , .
Similarly, the theoretical (critical) shear strength with fixed opening displacement is obtained from the first derivative of the energy-shear displacement curves in the inflection point: Note: The interface characteristic length l, the opening displacement x 0 which is corresponding to zero tensile stress for the shear displacement at us γ , the relative position of us γ , s, and the materials parameters K, p and q.
t c t n t n t , 0 t n (9) Combining equation (1), (7) and (8), an expression for the normalised tensile stress can be derived: Furthermore, from equation (5) and (9), we obtain the normalised relation between theoretical shear strength and the opening displacement: With equation (10) and (11), we can plot the theoretical shear strength as a function of normal stress (see section 5), to probe the postulated linear relationship.

τ σ / behavior
Several ab initio studies of the relationships between the tensile and shear strength in single crystals can be found in the literature [21][22][23]. For metals, a linear relation is widely assumed. From our new physical potential equation (5), a relation between the theoretical shear strength and the opening displacement can be derived as equation (11), and it shows an exponential decay of shear strength with increasing normal displacement. The relation between the scaled theoretical shear strength at each fixed opening displacement and normal separation are shown in figure 6(a). As one can see, the derived formula fits the ab initio reasonably very well. Next, we combine equations (5) and (7)

Crack propagation versus shear instability
The competition between crack propagation and grain boundary dislocation nucleation leading to grain boundary sliding can be evaluated following the approach of Rice [14]. His original criterion to judge whether a cleavage crack in a single crystal is blunted or not relies on the ratio of the surface energy γ S that would be needed to propagate the crack, and the unstable stacking fault energy in the slip plane, which would have to be overcome to emit a dislocation. Taking into account the geometric relation between crack propagation plane and slip plane the criterion becomes  (12) x and z denote fractional shear loadings, = x K K / II I , and = z K K / III I . K I , K II and K III are the stress intensity factors for different crack modes. θ represents the angle between slip plane and crack propagation plane, and φ is the angle between the direction of the Burgers vector and the crack propagation direction. To generalise the criterion to grain boundaries we replace γ s by G /2 c , so the term on the left hand side is transformed into 1/2q. Restricting ourselves to the pure K II mode (the grain boundary plane coincides with the crack plane, and the slip direction with the crack propagation direction), we set θ = 0 , φ = 0, x = 1, and z = 0. Then the criterion is simplified to 1/2 q > 2 which is equal to for dislocation nucleation. This means that if the q value is smaller than 0.25, a shear instability (which can be identified with dislocation emission) should occur at the crack tip, rather than crack propagation. From table 2, we can see that all easy sliding directions satisfy the criterion ⩽ q 1/4, so we would expect a shear instability to occur instead of crack propagation. This confirms our previous categorisation [24] of easy-sliding and non-easy sliding directions in the grain boundary planes, which was merely an intuitive choice based on the topology of the generalised stacking fault energy surface. Especially for the tilt grain boundaries this demonstrates again the strong anisotropy of the mechanical response with respect to the shearing direction. In all cases a shear instability is predicted parallel to the direction of the tilt axis, while fracture should occur when shearing the interface perpendicular to the tilt axis.
The constitutive behaviour which is derived in this section, is agreement with the literature, and the predicted deformation comes up to our expectations based on our previous work. However, a verification of the simple criterion for q to distinguish between crack propagation and dislocation nucleation, e.g. via molecular dynamics simulations, would of course be desirable. Such simulations would also serve to clarify the nature of the shear instability, which can be dislocation emission from the interface, rigid grain sliding, shear via grain boundary migration, or a combination of several of these processes. However, in general it is reasonable to assume that the lower the q value, the more likely is a shear instability.

Discussion
Our new physical model, equation (5), gives a sufficiently general description of the tensionshear coupling for combined sliding and opening displacements in the [1 1 2 ](1 1 1) slip system in single crystalline Al and in the shear directions in different grain boundaries. The closed form expression describes the tension-shear coupling very well, and has the new feature that the maximum of the energy barrier can be freely located. In the following we give a perspective on a unified description for all slip systems, and we discuss the transferability of our model to other metals.
Up to now the important parameters, the stable and unstable stacking fault energies and their corresponding displacements, as well as the parameters of the UBER, have to be determined individually for each slip system. It would of course be desirable to have just one common potential for each material. Therefore one would have to relate the energy barriers along different directions in the grain boundary planes to the intrinsic stacking fault energies of the metal, and the critical displacement and the cleavage energy of individual grain boundary planes to that of the single crystal (i.e. for zero misorientation) in the spirit of structure-energy models of grain boundaries. Progress has been made recently in finding a closed analytic expression for the energy of a grain boundary based on its geometric degrees of freedom, by fitting a Read-Shockley-type expression to atomistic data for hundreds of grain boundaries [2,7]. A similar approach could be tried for the characteristic quantitites in our potential, but is left for future work.
The mentioned structure-energy relationship seems to work well for different metals. Similarly, the transferrability of our potential is only limited by the limitation of the underlying approaches, the UBER and the expansion of the energy versus shear displacement in terms of trigonometric functions. The UBER has been shown to work well for interplanar potentials of rather different metals [10,12], including bcc metals with segregated impurities [11,34]. Only in the compressive region (which is described well for Al by our potential) the original analytic expression used by Rose [12] shows slight discrepancies [35]. Since our main interest is in crack propagation, i.e. tensile and shear opening displacements, we condone this shortcoming for the benefit of the simplicity of our expression.
The energy as a function of translation parallel to the grain boundary or cleavage plane can always be expressed in terms of the Fourier transform of the interatomic potential, because the two half crystals that meet at the boundary are periodic. A practical expression can be obtained if only pairwise interactions are assumed [25]. A comparable analytic function, which is at the same time an expansion of the approach of Sun et al, is our expression in terms of sine functions. The rather strong approximation of only pairwise interactions is partially compensated by the fact that we introduced the parameters K and s, which act as shape parameters to allow deviations from the ideal trigonometric funtions, and which implicitly include many-body effects if they are obtained by means of ab initio calculations.

Summary
We have established an analytical potential that describes the elastic response of grain boundaries to coupled tensile and shear mechanical loads, including even the so-called tensionsoftening effect. This interplanar potential can be used to derive constitutive relationships for continuum models of deformation and fracture at grain boundaries, as well as to quantify the interdependence of normal and shear strength. Furthermore, the deformation mechanisms (shear instability versus crack propagation) can be anticipated based on one of its parameters. Doing so we predict that for the tilt grain boundaries that have been investigated, a shear instability would occur along the previously identified easy-sliding directions parallel to the tilt axes, rather than crack propagation.
All parameters of our potential are physical quantities that can be determined directly on the basis of ab initio calculations and that enable an application to single crystals as well. In order to verify the practicality of the potential, we performed ab initio calculations of the single crystal (1 1 1) plane and four different grain boundaries in Al. The tension-shear coupling results obtained from ab initio calculations are in good agreement with the predictions of the analytical potential. The theoretical shear strength values decay exponentially with opening displacement for all systems, and more or less linearly with increasing tensile stresses.