Phonon structure of titanium under shear deformation along $\{10\bar{1}2\}$ twinning mode

We investigated phonon behavior of hexagonal close packed titanium under homogeneous shear deformation corresponding to the $\{10\bar{1}2\}$ twinning mode using first-principles calculation and phonon calculation. By this deformation, we found that a phonon mode located at a point on Brillouin zone boundary is drastically soften increasing the shear and finally it triggers a spontaneous structural transition by breaking the crystal symmetry toward twin from parent.


I. INTRODUCTION
The crystallography of deformation twinning has been studied by Bilby et al. 1 in 1965 and Bevis et al. 2 in 1968. Croker and Bevis also reported that for titanium specifically in 1970. 3 These studies were well summarized in the review paper by Christian and Mahajan in 1995. 4 These studies are based on the concept that the deformation twinning is characterized by collective displacement of atoms. 5 Therefore, phonon should play a central role in the deformation twinning. Recently, increase of computational power has enabled to perform first principles phonon calculation routinely. We applied this phonon calculation to the deformation twinning to discuss the structural transition from parent to twin in a quantitative manner as described in the following sections. In this report, we take the {1012} twinning mode of the hexagonal close-packed (HCP) titanium as a simple example to apply our general approach using the phonon calculation.
There are recent reports by using first principles calculations on microscopic mechanism of deformation twinning involving twin boundary. 6,7 However we take the traditional crystallographic approach to construct the calculation model with homogeneous shear without including the twin boundary. The first-principles phonon calculation is applied to this sheared-parent structure model to seek a characteristic phonon mode that exhibits structural instability under a certain shear. We show that the spontaneous structural transition induced by this phonon instability triggers transformation from the sheared-parent to twin. The final collective displacement of atoms predicted by this calculation excellently agrees with those reported in Refs. 1-4 as the shuffling mechanism.
Throughout this report, we employ symmetry constraints to distinguish the parent, sheared-parent, shuffling, sheared-twin, and twin structures. Here the shuffling structure is defined as the instantaneous structure that appears during rearrangement of atomic positions by the transformation from the sheared-parent to shearedtwin structure. More details are provided in the following sections. The aim of use of crystal symmetry is to analyze the mechanism of the deformation twinning in a comprehensible way by limiting the number of degrees of atomistic freedom. Thus an important contribution from this report may be detailed symmetry information provided along the transition pathway from the parent to twin structure.
This report is organized as follows. In Sec. II, a calculation model of the crystal structure is defined using a transformation matrix to represent the homogeneous shear corresponding to the {1012} twinning mode. In Sec. III, computational details of the first principles calculation and phonon calculation are given along with the convergence test against k-point density and smearing width of electronic structure calculation. In Sec. IV, the transformation from the parent to twin structure is discussed from a viewpoint of collective displacement of atoms. Finally, shuffling mechanism of the deformation twinning is associated with symmetry breaking of the sheared-parent structure induced by the unstable phonon mode.

A. Unit cell representations
The HCP crystal structure is shown in Fig. 1 (a). The set of the basis vectors is given as (a, b, c). To represent the {1012} twinning mode characterized by η 1 and η 2 directions, K 1 and K 2 planes, and the plane of shear P, it is convenient to retake the HCP unit cell by its extended unit cell as depicted in Fig. 1 (b). The set of the basis vectors of the extended unit cell, (a ′ , b ′ , c ′ ), is chosen to satisfy the conditions η 1 a ′ , η 2 c ′ , and b ′ ⊥ P. The change-of-basis is easily represented using an integer matrix Q: where the bars on the matrix elements denote the negative numbers. The conditions of b ′ ⊥ a ′ and b ′ ⊥ c ′ are easily confirmed finding (2a

B. Lattice under homogeneous shear
Following the usual crystallographic definition of deformation twinning, 4 we impose homogeneous shear corresponding to the {1012} twinning mode to the extended unit cell. In Fig. 2, the extended unit cell is illustrated being projected on the plane of shear P. The basis vectors a ′ and c ′ are parallel to the P plane, and the basis vector b ′ is perpendicular to the P plane.
The homogeneous shear is written by an affine transformation: where S is the matrix representation of the shear, and u and v are the vectors. By this homogeneous shear S, the basis vectors (a ′ , b ′ , c ′ ) are transformed to (a ′ s , b ′ s , c ′ s ) as follows: Following the convention by Christian and Mahajan in Ref. 4, to represent the {1012} twinning mode, we take a Cartesian axis l the unit vector parallel to η 1 and m the unit normal to the K 1 plane as shown in Fig. 2. Choosing the orientation by l = (1, 0, 0) T and m = (0, 1, 0) T , the basis vectors (a ′ , b ′ , c ′ ) are given in the Cartesian coordinates as The basis vectors a ′ and b ′ are unchanged applying this S, whereas c ′ is transformed to c ′ s = (c ′ x + sc ′ y , c ′ y , 0) T . It is more convenient for us to represent the homogeneous shear as a linear combination of the basis vectors since the choice of the orientation in the Cartesian coordinates becomes arbitrary. This is written by a transformation matrix T: To have c ′ s = c ′ + ta ′ as shown in Fig. 2, we take In the Cartesian coordinates (5), the relationship between s and t is given as This is found by combining Eqs. (4) and (7) as SL = LT, where L is the 3 × 3 matrix whose columns are a ′ , b ′ , and c ′ in the Cartesian coordinates (5). From Eqs. (1) and (7), we obtain the basis vectors under the homogeneous shear, (a s , b s , c s ), by Here Eq. (10) is independent on the orientation of Cartesian coordinates. The combined transformation matrix QTQ −1 is used to construct the sheared HCP unit cells for the first-principles calculation and phonon calculation.
C. Sheared-parent structure and shuffling At a certain twinning shear s = s t in the matrix (6), the sheared lattice gets to have the same lattice as the parent lattice with a different orientation, and thus the same point group symmetry as the parent lattice. The twinning shear of the {1012} twinning mode is known Sicne the ratio s/s t = t/t t is a good measure of the homogeneous shear in this study, we employ s/s t as the representative parameter in the following part of this report.
Even at s/s t = 1, the space group symmetry of the sheared-parent structure is unnecessarily the same as that of the parent structure. For the sheared-parent structure to have the same space group symmetry as the parent, i.e., to form the twin, usually a certain rearrangement of the atomic positions is required. This rearrangement may be called "shuffling", and for the use of this terminology clearly in our calculations, we define it in a little more detail as follows.
The space group type of the HCP crystal structure is P 6 3 /mmc (No. 194) and its site symmetry type (Wyckoff position) of all atoms is 6m2 (c). Obviously, the homogeneous shear for the {1012} twinning mode breaks these symmetries. Keeping crystallographic coordinates of atoms unchanged under the lattice shear, these symmetries are reduced to C2/m (No. 12) and m (i), respectively, except at s/s t = 1, Cmcm (No. 63) and m2m (c), respectively.
For simplicity, we assume that sheared structure is either sheared-parent structure or sheared-twin structure in the interval of {0 ≤ s/s t ≤ 1}, and the sheared-parent and sheared-twin structures may transform each other by rearranging their atomic positions. In this rearrangement, we expect that the atoms are displaced collectively and their displacement distances are much shorter than the distance between the nearest neighboring atoms.
We fix the crystal symmetry of the sheared-parent structure as C2/m. Unless breaking this symmetry, we allow to relax the positions of the atoms, by which the atoms can be displaced within their planes parallel to the plane of shear P. The atomic displacement u(s/s t ) by this relaxation is written as where r(s/s t ) is the position of the atom after the relaxation at s/s t and x(s/s t ) is the crystallographic coordinates corresponding to r(s/s t ). The second line of Eq. (11) is the convenient expression for us to construct the crystal structure models. This relaxation is not considered as a part of the shuffling in this study. We require that the shuffling arises in association with temporal symmetry breaking by the rearrangement of the atomic positions.

A. Computational details
The phonon calculations were performed by the finite displacement method with the supercell approach, where the atoms were displaced. The supercells were created by 4 × 4 × 3 multiplication of the sheared HCP unit cells obtained by QTQ −1 in Eq. (10). For each phonon calculation, atomic displacements of 0.02Å were introduced to the perfect supercells along a s and c s in plus and minus directions individually. For the phonon calculation, we employed the phonopy code. 8 For the first-principles calculations, we employed the plane-wave basis projector augmented wave method 9 within the framework of density functional theory (DFT) as implemented in the VASP code. [10][11][12] The generalized gradient approximation of Perdew, Burke, and Ernzerhof revised for solids 13 was used as the exchange correlation potential. A plane-wave energy cutoff of 300 eV was employed. The radial cutoff of the PAW dataset of Ti was 1.323Å. The 3p electrons for Ti were treated as valence and the remaining electrons were kept frozen. The smearing method was used with the k-point sampling on a uniform mesh for the Brillouin zone integration of the electronic structure. The 8 × 8 × 6 and 2 × 2 × 2 k-point sampling meshes with half grid shifts along c * directions were used for the sheared HCP unit cells and their 4×4×3 supercells, respectively, in conjunction with the smearing width σ = 0.4 eV in the Methfessel-Paxton scheme. 14 These parameters were chosen after convergence check of the lattice parameter and phonon band structure as presented in the next section.
To perform systematic calculations presented below, we employed the AiiDA environment 15 with the AiiDA-VASP 16 and AiiDA-phonopy 17 plugins.

B. Choices of calculation parameters
To choose the k-point sampling mesh density and smearing width σ, we examined convergence of the lattice parameter of the HCP conventional unit cell and the phonon band structure with respect to these values. In general, use of denser k-point sampling mesh provides better calculation accuracy at a constant smearing width σ. However, to save the computational demand, we expect to employ a sparser k-point sampling mesh. For this purpose, we may choose a larger σ value although sacrificing fine detail of electronic structure. We have to find a good compromise between the k-point sampling mesh and smearing width σ against the required accuracy of our computational research.
The calculated lattice parameters a and c with respect to the k-point sampling mesh of 8 × 8 × 6, 12 × 12 × 9, 16 × 16 × 12, and 20 × 20 × 15 for the HCP unit cell and σ = 0.1, 0.2, 0.3, and 0.4 eV are shown in Fig. 3. Clearly, the lattice parameter c is more sensitive to both of the k-point sampling mesh density and σ than a. At σ = 0.2, 0.3, and 0.4 eV, the lattice parameters are getting to converge increasing the k-point sampling mesh density. Use of σ = 0.1 eV requires an even denser k-point sampling mesh, for which the phonon calculation with the supercell approach can be computationally too demanding for us. Therefore, we decided not to consider the use of σ = 0.  The phonon band structure calculated with the 8×8×6 k-point sampling mesh (2×2×2 for supercell) and σ = 0.4 eV is compared with that with the 20 × 20 × 15 k-point sampling mesh (5 × 5 × 5 for supercell) and σ = 0.2 eV as shown in Fig. 4. These phonon band structures reasonably agree with the experiment in Ref. 19 and calculation in Ref. 20. Although we can see noticeable difference in phonon frequency between them, we consider this difference is small enough to discuss characteristic behavior of the phonon mode that we are interested in. Therefore we chose the k-point sampling mesh of 8 × 8 × 6 along with the smearing width of σ = 0.4 eV.

A. Optimization of sheared-parent structures
The crystal structure models of the parent, shearedparent, and twin structures are shown in Fig. 5. The parent and twin structures have the space group type of P 6 3 /mmc whereas the twin is oriented in a different direction from the parent after the shuffling. To create the sheared-parent structure models, we sampled 21 homogeneous shears in the interval of {0 ≤ s/s t ≤ 1} uniformly, and the first-principles calculations were performed to optimize their internal atomic positions under the symmetry constraint. The displacement distances |u(s/s t )| (see Eq. (11)) obtained as the results of the structure optimizations are shown in Fig. 6. Since |u(s/s t )| of the all atoms in the unit cell at each s/s t are the same due to the crystal symmetry, only one value at each s/s t is shown. The approximate directions of the displacements are shown in the inset. The atoms having the same displacement directions in this figure are symmetrically equivalent by the lattice translation, i.e., only two directions can exist and they are directed in the opposite directions as proposed by the crystal symmetry. The displacement distances |u(s/s t )| increase as increasing s/s t . Even at s/s t = 1, the displacement distance is only about 4% of the nearest neighbor distance. Therefore, for the schematic analysis, performing this optimization is considered unimportant, although it is necessary for accurate phonon calculation.  Fig. 5 (b). See the caption of Fig. 2 for the different symbols of atoms.

B. Energy increase by homogeneous shear
In Fig. 7, we show electronic total energy increase of the sheared-parent structure with respect to the homogeneous shear s/s t under the symmetry constraint. Increasing the shear, the energy increases harmonically. From the definition of the deformation twinning, we know that the sheared-twin structures have to give the same energy curve with respect to 1 − s/s t . This is drawn by the dotted curve as the mirror image of that of the parent. The crossing of these energy curves at s/s t = 0.5 implies that their energy surfaces are disjointed in the atomic configuration space under the symmetry constraints of C2/m.

C. Spontaneous symmetry breaking due to shear
At s/s t > 0.5, breaking the symmetry constraint may allow the sheared-parent structure to transform to the sheared-twin structure spontaneously. This is guessed due to the fact that the {1012} twinning mode exists. Therefore, we investigated this systematically as follows.
We applied a series of the phonon calculations to the sheared-parent structures at the 21 homogeneous shears that are the same shears as those chosen in Secs. IV A and IV B. The phonon band structures calculated at the selected shears are shown in Fig. 8. One phonon mode at the wave vector q = (1/2, 0, 0) exhibits imaginary frequency at the larger shears. This indicates that spontaneous structural transformation should occur by breaking symmetry. Note that the coordinates of q are represented in the basis vectors given by Eq. (10).
This behavior is similar to the structural phase transition induced by changing an intensive parameter such as pressure or temperature. 22 Assuming the second-orderlike structural phase transition, the squared frequency of the characteristic phonon mode is considered to be related to the order parameter. Hence we plotted it as a function of the homogeneous shear as shown in Fig. 9, and the roughly linear trend was found as expected. From this plot, the critical shear was estimated as s/s t ∼ 0.68.  Fig. 5 (b).

D. Shuffling
The polarization vector (eigenvector) of the characteristic phonon mode at q = (1/2, 0, 0) contains necessary information to properly break the crystal symmetry by introducing collective displacement of atoms that is illustrated in the inset of Fig. 7. In this figure, the arrows directed from the atoms in the sheared-parent structure depict directions of the atomic displacements with the same relative amplitude. This symmetry breaking doubles the primitive cell and the space group type is reduced to P 2 1 /c (No. 15).
Small finite displacements along the polarization vector were introduced to the sheared-parent structure models so that the first-principles calculation code, VASP, can detect the broken symmetry correctly. Then, we performed the first-principles calculation to optimize these structure models at the homogeneous shears with their basis vectors fixed. After the structure optimizations, we obtained the sheared-twin structures at s/s t = 0.7, . . . , 0.95 and the twin structure at s/s t = 1, i.e., the structures did not fall into local minima. This structural transformation is the shuffling that we consider.
The atomic displacements by this transformation at s/s t = 1 is shown by small (black) arrows in Fig. 10. The displacement distance is roughly constant by ∼0.5 A (∼17% of the interatomic distance) at all s/s t larger than the critical shear. The minimum structural unit of the collective displacement is composed of four atoms.
We can see that two types of the parallelogram units are arranged to fill the crystal structure by alternately changing their rotation directions depicted by the (red) circular arrows. The spaces surrounded by these rotation units behave like breathing. In this manner, it is considered that the internal structural distortion required by the shuffling is minimized. This illustration is almost the same as that obtained from the crystallographic discussion of the shuffling for the {1012} twinning mode by Crocker and Bevis in Ref. 3. What was found in this study is that this shuffling necessarily occurs at the shear larger than the critical shear.
Energy difference between the sheared-parent and sheared-twin structures is released by the shuffling as depicted by the arrows with the vertical dashed lines in Fig. 7. In between s/s t = 0.5 and the critical shear s/s t ∼ 0.68, the potential energy barrier should prevent the shuffling from being initiated. In real materials, imperfection of crystal may lower the potential energy barrier locally, and the shuffling initiated from this point will propagate over the whole crystal body at the velocity of sound by the help of the released energy to result in the macroscopic transformation. It would be difficult to observe how the twin forms from the parent microscopically, since the sheared parent relaxes toward the twin instantaneously. What we can measure by experiments as an evidence of this shuffling mechanism like a structural phase transition is the systematic frequency change of the characteristic phonon mode at s/s t < 0.5 as the precursor effect, for which inelastic neutron or Xray scattering techniques can be employed. The small (black) arrows depict atomic displacements induced by the structural optimization from the sheared-parent structure (Fig. 5 (b)) to the twin structure (Fig. 5 (c)). The circular (red) arrows show how the displacements are arranged as structural units.

V. SUMMARY
We investigated the microscopic mechanism of the {1012} twinning mode of the HCP titanium using the first principles calculation and phonon calculation. As the calculation model, the sheared-parent structures were defined by change of basis and transformation matrices at the specific homogeneous shears corresponding to the {1012} twinning mode. At the homogeneous shears of s/s t > 0.68, the phonon mode instability under the symmetry constraint was exhibited at a Brillouin zone boundary of q = (0.5, 0, 0). Breaking the crystal symmetry following the polarization vector (eigenvector) of this characteristic phonon mode, the sheared-twin structure was obtained after optimizing internal positions of atoms. This indicates the mechanism of the formation of the deformation twinning of this system. These calculations provide the structural transition pathway from the parent to twin, which may be widely accepted as the "shuffling". The collective displacement of atoms obtained from the phonon calculation agrees with the shuffling mechanism reported by Crocker and Bevis 3 from the crystallographic argument. We applied the same approach to the other twinning modes of HCP metals and found that the shuffling of the {1012} twinning mode is the simplest case. For the others, it may be hard for the simple crystallographic approach to find their reasonable collective displacements of atoms since those atoms are rearranged in more complicated ways three dimensionally than that of the {1012} twinning mode. These results of the other twinning modes will be published elsewhere. In this report, we provided a systematic calculational approach to explore microscopic mechanism of deformation twinning, and the technical details were written in this report. The computational resource and tools necessary to perform this calculation are usual. Therefore any researcher who is interested in the deformation twinning can achieve the similar calculations immediately and find specific details of their twinning modes. This is probably the largest contribution of this report to the field of deformation twinning.

ACKNOWLEDGMENTS
This work was supported by MEXT Japan through ESISM (Elements Strategy Initiative for Structural Materials) of Kyoto University.