Deformation mode analysis by eigenvectors of atomic elastic stiffness in static uniaxial tension of various fcc, bcc, and hcp metals

In order to clarify the physical meaning of the eigenvector of the atomic elastic stiffness matrix, B aij = Δ σ ai / Δ ε j , static calculations of uniaxial tension are performed on various fcc, bcc, and hcp metals with four different embedded atom method (EAM) potentials. Many fcc metals show instability for the constant volume mode, or the eigenvector of ( Δ ε xx , Δ ε yy , Δ ε zz ) = ( ± 1, ∓ 1, 0), under the [001] tension. Bcc also loses resistance against other constant volume mode, ( Δ ε xx , Δ ε yy , Δ ε zz ) = ( ± 1, ± 1, ∓ 2), in the [001] tension. Hcp shows shear modes Δ γ yz and Δ γ zx under the [0001] tension, which correspond to atom migration by dislocation on the slip plane. Similar shear modes appear in the [111] tension of fcc and [110] tension of bcc. Hcp also changes the mode to constant volume and shear in the [1010] tension, which imply the deformation in the pyramidal and prismatic planes.


I. INTRODUCTION
We have proposed a deformation mode analysis based on a 6 × 6 matrix of atomic elastic stiffness (AES), B a ij = Δσ a i /Δεj, for each atom to discuss local deformation such as crack propagation and dislocation glide. [1][2][3][4] Here, Δσ a i is the change in the atomic stress, Δεj is that in the strain, and the subscripts i, j are the Voigt notation 5 of i, j = 1 ∼ 6 = xx, yy, zz, yz, zx, xy. We can construct the eigenequation of B a ij Δεj = η a Δεi for each atom and obtain six eigenvalues η a(1) , . . . , η a (6) . This process is equivalent to the diagonalization of the matrix B a ij , and the six eigenvalues are the diagonal elements aligned in ascending order. If atoms have negative eigenvalues, they can be regarded as "unstable" from the relation of Δσ a i = B a ij Δεj = η a Δεi. The corresponding eigenvector {Δεi} = {Δε 1 , . . . , Δε 6 } T = {Δεxx, . . . , Δγxy} T for η a(k) < 0 atoms is also expected to reflect the deformation path. Here, superscript T means transpose to denote that the column vector is shown in the row vector form. In fact, the crack opening mode and the glide direction of the dislocation core were visualized by the principal axes of the strain tensor [Δεij] (where i, j are the Cartesian indices of x, y, and z) of which components are those of the eigenvector {Δεi} of the η a (1) < 0 atoms. [2][3][4] Nishimura et al. 6,7 also evaluated B a ij in the buckling simulation of a multi-wall carbon nanotube (MWCNT) and reported that the global buckling of MWCNT is triggered by the emergence of λ a (2) < 0 atoms (λ a (2) corresponds to η a (2) ), after all the atoms become λ a (1) < 0 in MWCNT. 6,7 They also revealed that the corresponding eigenvectors are {0, 0, 0, 0, γzr, 0} T for λ a (1) < 0 (in-plane shear in each CNT surface) and {0, 0, 0, 0, 0, γ rθ } T for λ a (2) < 0 (out-of-plane shear) in the (r, θ, z) cylindrical coordinate. Thus, the eigenvector definitely suggests deformation mode; however, the diagonalization is performed in the six-dimensional space of xx, yy, zz, yz, zx, xy so that the eigenvalue may have dependency on the coordinate, in contrast to the coordinate independence of the principal strain/stress of the tensor.
In the present study, the characteristics of the eigenvalue and eigenvector of B a ij are discussed in static uniaxial tension for two different crystal orientation, in order to show the fundamental aspects of these quantities. Eight fcc, four bcc, and four hcp elements were calculated with embedded atom method (EAM) potential by Zhou et al. 8,9 and compared with different EAM potential for Ni and Al by Voter et al., 10,11 Fe by Finnis-Sinclair, 12 and Mg by Liu et al. 13

II. EIGENVALUE AND EIGENVECTOR IN PERFECT LATTICE
As already indicated in our previous papers, 1-4 the mathematical form of the atomic stress, σ a i , and the atomic elastic coefficients, C a ij , are easily formulated in the EAM potential; thus, both are directly calculated only with the atom configuration at each impact. B a ij can be evaluated with C a ij and the stress contribution, e.g., B a 11 = C a 11 +σ a 1 , B a 12 = C a 12 −(σ a 1 +σ a 2 )/2, and B a 44 = C a 44 +(σ a 2 +σ a 3 )/2. In the no-load perfect lattice, B a ij coincides with Cij or the elastic coefficient of the crystal since all the atoms are same. Supposing isotropic elasticity, the Cij matrix is expressed as follows: (1) ±0.58 ±0.58 ±0.58 0 0 0 η (5) ±0.71 ∓ 0.71 0 0 0 0 η (6) ±0.58 ±0.58 ±0.58 0 0 0 ±0.58 ±0.58 ±0.58 0 0 0 ±0.71 ∓0.71 0 0 0 0 η (6) ±0.58 ±0.58 ±0.58 0 0 0 The characteristic equation of the 3 × 3 submatrix of the upper left is (2) Thus, six eigenvalues of Cij become as follows: The corresponding eigenvector for the eigenvalue η = C 11 + 2C 12 is Δεxx = Δεyy = Δεzz. For the double root of η = C 11 − C 12 , the corresponding eigenvector satisfies the relation Δεxx + Δεyy + Δεzz = 0 or the volumetric strain is zero. Presumable solutions are (Δεxx, Δεyy, Δεzz) = (1, −1, 0) and (− . Except Pt and Al, fcc elements show C 11 − C 12 < C 44 so that the 1st and 2nd eigenvalues are C 11 -C 12 . The corresponding eigenvectors are listed in Table II. Here, a small strain of εzz = 0.001 is given to avoid ambiguity in the solution due to the perfect symmetry in this coordinate. That is, the eigenvectors for the double root η = C 11 − C 12 and triple root η = C 44 are not uniquely defined so that math library (LAPACK) returns an equivalent arbitrary solution. The shear components Δγyz and Δγzx still have rotational symmetry against the z-axis so that the value is different with element by element; however, they are equivalent to the solution of (Δγyz, Δγzx) = (1, 0), (0, 1). The direction for normal strain components coincides with the prediction in the isotropic Cij mentioned above. Many fcc elements of C 11 − C 12 < C 44 show lowest resistance against normal strain components, (Δεxx, Δεyy, Δεzz) = (1, −1, 0), (−1, −1, 2). The former corresponds to proportional elongation and shrinkage in the lateral direction, and  Table III. Contrary to the previous fcc case, only Fe has the relation of C 11 − C 12 < C 44 so that the other elements have the lowest eigenvalue of η (1) = η (2) = η (3) = C 44 in the shear without normal strain components. Table IV shows the eigenvector under a small strain of εzz = 0.001. They are almost the same as in Table II, except the order in the shear mode. That is, the shear Δγxy is always located at the top in the three shear modes, while it is at the bottom in Table II. Although the deformation resistance of C 44 is the same for the three modes, the xy shear is definitely preferable in bcc metals, while it is not in fcc ones against strain perturbation in the z-direction.
The Cij of hcp metals are listed in Table V. Hcp metals have anisotropy about C 33 , C 13 , and C 66 so that the characteristic equation The solution is 13 2 , All hcp elements have the C 11 − C 12 > C 44 , C 66 relation. The eigenvalues and eigenvectors of Zhou's Mg are shown in Table VI. Here, no strain perturbation is applied in the calculation of the eigenvector; however, the ambiguity for Δγyz and Δγzx vanishes due to the anisotropy of the hcp structure. The same deformation mode of (Δεxx, Δεyy, Δεzz) = (1, −1, 0), (−1, −1, 2), and (1, 1, 1) as isotropic fcc and bcc is obtained for the eigenvalues of Eq. (5). The eigenvector is the same for all other elements in Zhou's potential. Eigenvectors for η (1) and η (2) are switched in Liu's potential, but it might have no physical significance due to the same value of η (1) = η (2) = C 44 . implemented with 10 × 10 × 10 unit lattices, and the lateral strain εxx and εyy is controlled to reduce the normal stresses to |σxx|, |σyy| < 0.1 [GPa]. η (6) = C 11 + 2C 12 shows an extremely high value against the other eigenvalues so that it is omitted in the graph. The stress shows a rather unrealistic curve due to the energy cutoff under unrealistic large strain. The components of the eigenvector could be changed with the loading; however, the deformation mode is labeled by the initial eigenvectors. The "instability" is caused by the (Δεxx, Δεyy, Δεzz) = (1, −1, 0) mode in all fcc elements of Zhou's potential and Ni of Voter [first crossing point over the η a = 0 line in Figs. 1(a)-1(c)]. This instability triggers anisotropic shrinkage/expansion in the xy-cross section and leads to local deformation or dislocation nucleation in actual molecular dynamics simulation. 14 On the other hand, the (Δεxx, Δεyy, Δεzz) = (−1, −1, 2) mode becomes lowest and almost zero at εzz = 0.072 in the case of Zhou's Al. The stress plateau is also observed around the point. In fact, we performed molecular dynamics simulation of the Ni/Al bilayer with Zhou's potential and observed amorphization in the Al phase at small strain, which is rather different from our previous report 15 by Voter's EAM. In the case of Voter's Al, the Δγyz and Δγzx modes are kept lowest and reach η = 0 first at εzz = 0.136. However, the (Δεxx, Δεyy, Δεzz) = (−1, −1, 2) mode also becomes negative at εzz = 0.137 and turns to be the lowest with a steep drop of the eigenvalue. Zhou's Al has a local minimum for the mode as mentioned; however, the eigenvalue quickly recovers and the change in the other eigenvalues crosses each other after the strain. These complicated mode changes may cause the amorphization in Zhou's Al.

III. UNIAXIAL TENSION IN [001] AXIS
The results for bcc Fe, Mo, and Ta are shown in Fig. 2. W is similar to Mo and omitted to compare Fe of Zhou and FS. Before the loading, the shear modes are lowest in Mo, Ta, and W; however, the instability is caused by the (Δεxx, Δεyy, Δεzz) = (−1, −1, 2) mode in all elements and potential type. Contrary to the previous fcc case, this mode may cause phase transformation since the unstable deformation occurs simultaneously in three normal directions. The unstable point coincides with the stress peak. The corresponding eigenvector at the unstable strain is of (Δεxx, Δεyy Here, the absolute sign and the magnitude of the strain components do not have any physical significance due to the arbitrary sign of the eigenvector normalized as √ Δε 2 xx + ⋯ + Δγ 2 xy = 1. Figure 3 shows the results of hcp Mg. The tendency is thoroughly different with the potential type. The other hcp elements, such as Co, Ti, and Zr of Zhou's potential, are quite similar to those in Fig. 3(a). That is, the lowest shear modes Δγyz and Δγzx reach η < 0 first, and then, the lowest switches to the mode (Δεxx, Δεyy, Δεzz) = (−1, −1, 2). These points are close to the strain of the stress peak. On the other hand, all the eigenvalues show the sudden change just before the stress-strain peak, and three modes of Δγyz, Δγzx, and (Δεxx, Δεyy, Δεzz) = (−1, −1, 2) simultaneously become negative at the stress-strain peak in Liu's Mg.

ARTICLE
scitation.org/journal/adv  to the slip plane. Returning to Fig. 4(d), the constant volume mode of (Δεxx, Δεyy, Δεzz) = (−1, −1, 2) rapidly drops in the Voter's Al and reaches η < 0 before these shear modes on the slip plane. This tendency is the same as the previous [001] tension of Fig. 1(d). Voter's Al might be rather "isotropic" in deformation resistance compared to the other fcc elements. Despite this difference at εzz = 0, the first instability is caused by the Δγzx mode in all elements and potential, as shown in Fig. 6 (W is omitted). This mode corresponds to the atom migration on the slip plane normal to tension, as schematically shown in Fig. 5(b). We also refer to that the stress-strain peak coincides with the point where η (2) becomes negative, under this "static" condition.  Table IX. The eigenvector is the same for all the elements including Liu's Mg. That is, all hcp elements show the lowest mode in Δγxy and Δγyz in this coordinate before tension. However, the instability is caused by the (Δεxx, Δεyy, Δεzz) = (1, 0, −1) mode in Zhou's potential, as shown in Fig. 7(a), while the eigenvector changes to (Δεxx, Δεyy, Δεzz) = (−0.03, 0.34, −0.94) at the instability point. This mode may correspond to the deformation on the pyramidal plane, as schematically illustrated in Fig. 5(c). Zhou's other elements show a similar change as in Fig. 7(a). Liu's Mg shows simultaneous instability of (Δεxx, Δεyy, Δεzz) = (0.34, −0.07, −0.94) and Δγyz shear, as shown in Fig. 7(b). The former also may reflect the pyramidal plane slip. The latter represents the shear parallel to the prismatic plane, which is normal to the tensile direction.

Mg Co Ti Zr
The eigenvalues of Liu's Mg show a discontinuous curve, especially in Fig. 7(b). The cusp of η a(1) definitely corresponds to the change in the number of surrounding atoms within the cutoff radius; that is, the eigenvalue rapidly drops just before the top and bottom neighbor atoms go out of the cutoff sphere by tension. Then, the eigenvalue recovers with fewer local density ρ and corresponding embedding energy F(ρ) given by the spline function. Liu's EAM also represents the 2-body pair potential ϕ(r) and the contribution of atomic densityρ(r) by the spline function so that it is sensitive against the change in the surrounding atoms. Voter's EAM also uses a spline function for F(ρ) so that we can see the same tendency in Fig. 1(d).

V. CONCLUSION
Static calculations on the eigenvalues and eigenvectors of a 6 × 6 matrix of atomic elastic stiffness, B a ij = Δσ a i /Δεj, are implemented for fcc, bcc, and hcp elements with different EAM potentials. First, the eigenvalues and eigenvectors of the elastic coefficients Cij are discussed as the characteristics before tension. The isotropic fcc and bcc have the eigenvalues of double root C 11 − C 12 , triple root C 44 , and C 11 + 2C 12 . The corresponding eigenvectors are constant volume modes, (Δεxx, Δεyy, Δεzz) = (±1, ∓1, 0) and (±1, ±1, ∓2), for the double root C 11 − C 12 , each shear Δγyz, Δγzx, Δγxy for triple root C 44 , and dilatancy mode Δεxx = Δεyy = Δεzz for C 11 + 2C 12 . Anisotropic hcp also shows the same modes. Then, the eigenvalues and eigenvectors of B a ij (=Cij except the σ 3 related term) are discussed under a static uniaxial tension of fcc [001], bcc [001], and hcp [0001]. Fcc elements basically show the instability for the mode of (Δεxx, Δεyy, Δεzz) = (±1, ∓1, 0), which leads to proportional expansion/shrinkage in the lateral direction. Bcc shows the instability for the mode of

ARTICLE
scitation.org/journal/adv (Δεxx, Δεyy, Δεzz) = (±1, ±1, ∓2) that may correspond to phase transformation. Hcp shows the instability of Δγyz and Δγzx shear, which reflects the atom migration on the basal slip plane. Another calculation of the uniaxial tension for fcc [111], bcc [110], and hcp [1010] reveals that the instability changes to shear mode on the slip plane in fcc and bcc. Hcp also changes its mode to constant volume mode and shear; the former may reflect the deformation on the pyramidal plane, and the latter may reflect the deformation on the prismatic plane.
Due to the potential problem and static analysis, the obtained results may not have significance for real material behaviors; however, the present paper could be a milestone for the precise discussion by ab initio density functional theory calculation (revision of our previous report 16 ) or for the basic knowledge against complicated local phenomena in molecular dynamics simulations.