Elastic model for stress–tensor-induced martensitic transformation and lattice instability in silicon under large strains

ABSTRACT Recent molecular dynamics simulations for phase transformations (PTs) Si ISi II under multiaxial loading revealed that lattice instability criteria for Si ISi II PTs are linear in the stresses normal to the cubic faces; instability stress depends on and is independent of shears; for some stresses , the instability stresses are the same for direct and reverse PTs; at critical stress , the first-order PTs are substituted with the second-order PT. We demonstrated that these features can be qualitatively reproduced using the simplest model with quadratic elastic energy in terms of Lagrangian strain. IMPACT STATEMENT Various nontrivial features of the martensitic transformation and lattice instability under complex loading in silicon are reproduced using the simplest model with quadratic elastic energy in terms of Lagrangian strain. GRAPHICAL ABSTRACT


Introduction
Displacive stress-induced phase transformations (PTs) under both normal and high pressures are broadly studied from both fundamental and applied points of view. They are important for understanding and improvement of behavior of shape memory alloys and caloric materials [1][2][3][4] and for quantifying the effects of nonhydrostatic stresses on high-pressure PTs [5][6][7][8][9], which can be utilized for the synthesis of new phases. The starting point in studying PTs at low temperature is finding stresses at which the crystal lattice loses its instability. There are several approaches that describe CONTACT Valery I. Levitas vlevitas@iastate.edu Departments of Aerospace Engineering, Mechanical Engineering, and Material Science and Engineering, Iowa State University, Ames, IA 50011, USA; Division of Materials Science and Engineering, Ames Laboratory, Ames, IA 50011, USA lattice instability under stresses. In the most traditional approach [10][11][12][13][14][15] (named below as the zero-moduli condition), lattice instability occurs when some elastic moduli or their combinations reduce to zero. For complex lattices, intracell displacements (shuffles) are included in the instability criteria [13][14][15][16]. A lattice instability criterion with respect to change in the order parameters within a Landau-type theory was developed in [17][18][19] and was coined as the modified transformation work criterion. Also, phonon stability (soft-mode) criteria [13][14][15][16] were applied, which offer additional conditions for short wavelengths.
(a) While the instability stress for Si I→Si II PT is described by both zero-moduli and the modified transformation work criteria, the reverse Si II→Si I PT is described by the modified transformation work criterion only.
(b) Instability criteria for both Si I↔Si II PTs are linear in the true (Cauchy) stresses σ 1 , σ 2 , and σ 3 normal to the cubic faces of Si I; i.e. they represent two planes in the stress space σ i (i = 1,2,3).
(c) The effect of the shear stresses within the considered range ≤ 3 GPa is negligible.
(d) The instability stress σ 3 along the c direction of the tetragonal phase or direction in a cubic lattice, which transforms to the c direction, depends on σ 2 + σ 1 rather than on σ 1 and σ 2 separately.
(e) Instability points for Si I→Si II and Si II →Si I PTs correspond to the maximum |σ max 3 | and minimum |σ min 3 | stresses, respectively, in the stress -Lagrangian strain E 3 curves for different fixed σ 1 and σ 2 .
(f) For some tensile stress σ 1 = σ 2 , the instability line for the direct σ max 3 (σ 1 ) PT and reverse σ min 3 (σ 1 ) PT intersects. For larger tensile stresses, while the straight line σ min 3 (σ 1 ) remains the same, the line σ max 3 (σ 1 ) bends toward and then coincides with σ min 3 (σ 1 ). In this region, the stresses for the instability of the crystal lattice are the same for direct and reverse PTs, which leads to a number of interesting and important phenomena [20].
(g) The magnitude of the strain in the c direction, at which instability starts, |E max 3 |, increases with increasing tensile stresses σ 1 = σ 2 (and reduces with the growing magnitude of the compressive stresses |σ 1 | = |σ 2 |), and jump in strain from the parent to the product phase decreases as well. At some σ 1 = σ 2 the jump becomes zero, which means that the first-order PT is substituted with the second-order PT.
All the above results create the impression of a very sophisticated nonlinear elastic material behavior for large strains at the edge of lattice instability under multiaxial loading. In particular, the coincidence for the Si I→Si II PT of two completely independent criteria, the zeromoduli condition and the modified transformation work criterion, may impose very strict constraints on nonlinear elasticity equations. Even without intending to describe the above features and constraints, the elastic energy for a cubic to tetragonal PT is considered to be a fourth-degree polynomial in terms of the Lagrangian strain tensor E [21].
In this letter, we demonstrate that all the enumerated features of Si I→Si II PT can be qualitatively reproduced using the simplest finite-strain model with quadratic elastic energy in terms of Lagrangian strain E. This results in a linear relationship between the second Piola-Kirchhoff stress T and E. However, under compression, e.g. along axis 3, it produces curves for the first Piola-Kirchhoff (nominal) stress P 3 (force per unit reference area) and the Cauchy (true) stress σ 3 (force per unit actual area) vs E 3 with the maximum (see Figure 1, similar to figures in [19,20]). The maximum of the Cauchy stress corresponds to the lattice instability point and the initiation of the direct PT Si I→Si II [19,20]. The reverse PT is also described in the simplest way, yet still consistent with MD results in [19,20]. Such a surprisingly simple description of multiple complex phenomena in PTs and lattice instabilities reveals that the main reason for all of them is geometric rather than physical nonlinearities.
Below, we designate contractions of tensors A = {A ij } and B = {B ji } over one and two indices as A · B = {A ij B jk } and A : B = A ij B ji , respectively. The superscript T designates transposition and I is the unit tensor. Because tensile (compressive) stresses and strains are positive (negative) and our main interest is in compression, we use the magnitude | · · · | and in many cases change the sign of stresses and strains, similar to [19,20].

Model
Let the deformation gradient F represent the linear mapping of the undeformed (stress-free) cubic crystal lattice in the deformed state, in particular, into a deformed tetragonal lattice, and the Lagrangian strain is E = 0.5(F T · F − I). The Jacobian determinant describes the ratio of volumes V in the deformed and the undeformed configurations, J = dV/dV 0 = det F. We will accept the simplest finite-strain model with quadratic elastic energy per unit undeformed volume where C is the fourth-rank elasticity tensor for the cubic lattice. Relationships between the second Piola-Kirchhoff stress T, the first Piola-Kirchhoff (nominal) stress P and the Cauchy (true) stress σ and E and F in nonlinear elasticity are defined as [22] It is clear that the relationship between T and E is linear, i.e. this is indeed the simplest model in nonlinear elasticity. However, expressions for P and σ are nonlinear and for large strains, they will reproduce all of the above features related to lattice instability and PTs. The reverse PT in [19,20] is not related to the elastic instability, but rather is described within a Landau-type theory by the modified transformation work criterion. In order to avoid a complex formulation based on an order parameter [19,20] which is not necessary for this paper, we will utilize some results from [19,20] to develop the simplest description. According to Figure 3 in [20], the elastic behavior of Si II is linear with very large elastic moduli, and stress-strain curves for Si II for several σ 1 = σ 2 represent the single (almost) vertical line. Thus, elastic deformations of Si II are negligible in comparison with deformations of Si I, and we will neglect them. Then strains in Si II are constant and equal to the strains of the stress-free Si II. According to Levitas et al. [19,20], deformation gradient F r for the stress-free Si II has components F r 1 = F r 2 = 1.175 and F r 3 = 0.553, which correspond to the Lagrangian strain E r with components E r 1 = E r 2 = 0.19 and E r 3 = −0.35. Thus, we will assume that PT Si I→Si II completes and that the reverse PT Si II→Si I starts at fixed strain E r , and that the product phase has infinite elastic moduli ( Figure 2). Despite the simplicity, this corresponds well with the MD results in [19,20] and is fully sufficient for our objectives here.

General equations
We will focus on the loading of a crystal of cubic symmetry by three stresses normal to the cubic faces. In this case, tensors F, E, T, F, and σ are coaxial and diagonal and Equation (2) for the Cauchy stress simplifies to or in component form for σ and P: Explicit expressions for the cubic symmetry are where C 11 and C 12 are elastic constants for the cubic crystal. All stresses and C 12 can be normalized by C 11 . To avoid additional notations, we just set C 11 = 1 and designate C = C 12 /C 11 . In addition, we express all F i = √ 1 + 2E i , and obtain from Equation (5): We will also need

Uniaxial loading of a weakly compressible material
To obtain simple results independent of the material parameters, we assume that J−1 is small and that the volumetric part of geometric nonlinearity can be neglected. Thus, we set J 1, and consider uniaxial tensioncompression along axis 3 at σ 1 = σ 2 = 0. Then it follows from Equations (5) and (7) A plot of all three stresses normalized by B is shown in Figure 1 In [19,20], lattice instability was identified at the maximum of the Cauchy stress, independent of which stress or strain was prescribed. This was proven analytically [18] for the Landau-type description of a PT in terms of the order parameter. If the prescribed Cauchy stress |σ 3 | is greater than |σ max 3 |, then the system cannot be equilibrated and the magnitude of strain will increase toward  |. In terms of the order parameter η, which is 0 in the parent phase and 1 in the product phase [19], the parent phase loses its stability and PT starts for | under prescribed P 3 means that this state is an equilibrium intermediate homogeneous state with 0 < η < 1. However, for heterogeneous perturbations, e.g. in MD simulations, multiple two-phase stationary solutions can be observed for |E Cmax 3 | < |E 3 | < |E Nmax 3 | and prescribed P 3 . Thus, the maximum Cauchy stress |σ max 3 | will be considered as the stress for lattice instability and initiation of direct PT. The same criterion was used in most of previous publications for uniaxial loading, e.g. in [11,12].

3D loading with σ 1 = σ 2
In this case, also E 1 = E 2 and analytical solution can be obtained as well. Equation (6) simplifies to Solving Equation (9) for E 1 , and substituting Equation (12) in Equations (10) and (11), we obtain Plots of the stresses σ 3 and P 3 vs. strain E 3 for stress σ 1 = 0 and C = 0.3855 are shown in Figure 1(b). Although compression magnitude of both stresses possess maxima, tensile curves are monotonous (similar to MD simulations in [19,20]). The maximum of the true stress, The effect of different fixed stresses σ 1 on the stress σ 3 -strain E 3 curves is shown in Figure 2. Qualitatively, it is similar to what is found in MD simulations (see item (g) in the Introduction). Thus, tensile stresses σ 1 increase the instability strain |E Cmax 3 | and decrease the instability stress |σ max 3 |, lowering stress-strain curves. Compressive σ 1 acts in the opposite way. For the tensile E 3 , stress-strain curves are monotonous for all σ 1 . One of the interesting features of Equation (13) is that for E 3 = C, the stress σ 3 = E 3 √ 1 + 2E 3 is independent of σ 1 and all curves for different σ 1 pass through this point (Figure 2(a)).
The instability stress |σ max 3 | for direct PT vs. σ 1 is shown in Figure 3. (ii) The curve |σ max 3 |(−σ 1 ) bends toward the intersection point of the straight lines for the direct and reverse PTs, and then there is a significant region σ 1 ⊂ [−0.2, 0.02], where both curves practically coincide (item (f)). In this region, the stresses for the instability of the crystal lattice are the same for direct and reverse PTs, i.e. the energy barrier between phases disappears and the Gibbs energy possesses a flat portion. This leads to unique homogeneous, hysteresis-free and dissipation-free first-order PTs, for which each intermediate lattice along the transformation path is in indifferent thermodynamic equilibrium and can be arrested and studied at fixed strain [20]. The above properties are of great interest from the fundamental and applied points of view, especially for various PT-related applications [1][2][3][4].  Figure 3 in [20] and item (g) above) to zero at σ 1 = 0.088. At this point, the first-order PT is substituted with the secondorder PT, as reported in [20] and item (g).

Effect of shear stress
In addition to the normal stresses, shear stress τ = σ 13 can be applied. To exclude rigid-body rotation, we accept that the components F 31 = 0. Then the deformation gradient and Lagrangian strain have the following form: Such an expression for F (even if we additionally include F 12 and F 23 ) does not change J = det F. In the MD simulations [19,20], the effect of the shear stress σ 13 = 3 GPa on the lattice instability did not exceed ±2. Since for Si I shear modulus C 44 = 80 GPa [23], an estimate for the shear strain is E 13 = σ 13 /(2C 44 ) = 0.01875 1. Then, Equation (2) for the Cauchy stress has the following form: We kept the term with F 13 in the expression for σ 13 because T 3 T 13 . To estimate E 12 for the geometrically nonlinear formulation, we express T 13 from Equation (16), i.e. from σ 13 = J −1 F 3 (F 1 T 13 + F 13 T 3 ), while allowing for T 3 = σ 3 J/F 2 3 and F 13 = 2E 13 /F 1 = T 13 /(C 44 F 1 ) (see Equation (15)): For uniaxial compression (σ 1 = 0) at instability point, 58 GPa, and, according to Equation (17), T 13 = 5.92 GPa, which is significantly larger than σ 13 = 3 GPa used in the linearized theory. Still E 13 = T 13 /(2C 44 ) = 0.037 1 and F 13 = 2E 13 /F 1 = 0.068 1, and we will see that their effect on the critical stress |σ max 3 | can be neglected. Indeed, an exact expression for σ 3 in Equation (16) is independent of F 13 . An exact equation for σ 1 in Equation (16) differs from the approximate by three terms equal to 0.004 (normalized by C 11 ). According to the expression |σ max 3 | = 0.128 − 0.382σ 1 , this changes |σ max 3 | by 0.0016, i.e. by 1.2%, which is negligible. Also, in expression for E 3 in Equation (15), the term due to shear F 2 13 = 0.0046 is also just 1.5% of the value 0.31 without shear.
Thus, the effect of shear stresses of 3 GPa on the instability stress for uniaxial compression is negligible in the current model, similar to the MD simulation results in [19] (item (c)).

Concluding remarks
In the letter, we demonstrated that all numerous nontrivial features found in MD simulations [19,20] for cubic-tetragonal PTs Si I↔Si II under general multiaxial loading can be qualitatively reproduced using the simplest model with quadratic elastic energy in terms of a Lagrangian strain. These phenomena include: (a) lattice instability criteria for both Si I↔Si II PTs are linear in the Cauchy stresses σ i (i = 1,2,3); (b) the instability stress σ 3 depends on σ 2 + σ 1 rather than on σ 1 and σ 2 separately; (c) the negligible effect of the shear stresses; (d) coincidence of the instability stresses σ 3 for direct and reverse PTs for some range of tensile stresses σ 2 = σ 1 , which results in a variety of interesting phenomena [20]; (e) a jump in strain E 3 from the parent to the product phase reduces with growing tensile stresses σ 2 = σ 1 until it approaches zero, which substitutes the first-order PT with the second-order PT.
Such a success in a description means that all of the enumerated phenomena result from geometric nonlinearity rather than any physical mechanisms or effects. This model can be included in the phase field simulation of stress-induced PTs and microstructure evolution, as well as analyses of the effect of nonhydrostatic stresses on high-pressure PTs.
The results obtained allow us to answer the question raised in [19,20]: why two different lattice instability criteria, one based on the elastic instability and another one that is derived within Landau theory, coincide for direct Si I to Si II PT and both are linear in terms of the components of the Cauchy stress? MD simulations in [19,20] show that instability for direct Si I to Si II PT occurs at zero elastic modulus and the corresponding criterion is linear in the components of the Cauchy stress. When applying the Landau approach, we eliminated terms that give nonlinear contribution; two material parameters in the remaining linear criterion were calibrated utilizing MD results. Here, we obtained the linear in stress instability criterion based on elastic instability and the simplest finite-strain elastic model.
Introducing the third degree contribution to the energy in terms of E may be sufficient for the quantitative description of the elastic behavior observed in MD simulations; this will result in quantitative coincidence with MD results on instability stresses and, consequently, with results based on Landau theory. Similar approaches can be developed for other materials that experience large strains during stress-induced transformations.

Disclosure statement
No potential conflict of interest was reported by the authors.