Atomistic Simulation on the Twin Boundary Migration in Mg under Shear Deformation

In this paper, the {101¯2} twinning and detwinning was studied by molecular dynamics simulation under different shear directions and strain rates. The results showed that the twin was thickened under [1¯011] shear direction and shrunken with shearing in the opposite direction. The critical resolved shear stress of {101¯2} twin boundary migration increased with the increase of the strain rate. By analyzing the atom’s displacement, it was concluded that the {101¯2} twin migration was achieved by both the shear and the atomic shuffling. Every atom would be affected by the shear, and different shear directions would cause opposite move directions, which led to twinning or detwinning. The atom shuffling was only used for adjusting the glide twin boundary and mirror-symmetric twin boundary structure evolution.


Introduction
Twin boundary (TB) is a special kind of two-dimensional planar defect in crystalline materials. It separates two crystalline regions structurally as the mirror images of each other. The highly symmetrical discontinuity in structure can be produced by deformation or annealing. Inevitably, the structure and the properties of TBs deeply affect the mechanical properties of materials. Therefore, the thermodynamics and kinetics of twin boundary are central to understanding the properties and evolution of many materials systems.
Magnesium (Mg) and its alloys, due to their low density and high specific strength, have potential applications in automotive and aerospace industries [1]. However, low strength and ductility induced by their hexagonal closed-pack (HCP) structure are the major obstacles preventing Mg alloys being widely used. As is well known, the basal slip in Mg is the only preferable slip system to activate at room temperature, which is insufficient to accommodate general homogeneous plastic deformation according to the Von Mises criterion. Twinning is the other important deformation mechanism for Mg at room temperature [2][3][4][5]. In addition, since it can adapt to the deformation along the c-axis of HCP crystalline, the introduction of twinning deformation can significantly enhance the strength and ductility of Mg [6][7][8][9][10][11]. Among all the possible twinning systems, the relatively lower shear critical resolved shear stress (CRSS) of the {1012} extension twinning mode and the high mobility of the extension twin boundary (ETB) make it the most frequently observed during various deformation of Mg alloys [2,3].
Considerable attention has been focused on the improvement of mechanical properties via extension twinning induced by predeformation [6,7,9]. He et al. [6] introduced {1012} tension twin in

Materials and Methods
In the investigation presented here, simulations were performed using molecular dynamics techniques. The simulation represents a small area of extension twin migrating through a deformed matrix under the external shear deformation. Compared to the extension twin microstructures observed in some experiments, the simulation cells are somewhat simplified, while still capturing the essential features of growth/shrinkage of the extension twin and allowing for a thorough analysis of the migration process of TB.
The potential used here is an embedded atom method developed by Sun [38], which has been successfully used to investigate many fundamental issues in Mg, such as the crack growth, transformation of twin boundary, etc. The initial setup consisted of three parts "patched" together in the simulation cell. Periodic boundary conditions were applied along the x-and z-directions; while the free boundary condition was applied along the y-direction. Figure 1 shows the simulation cell set up.
Materials 2019, 12, x FOR PEER REVIEW 3 of 11 transformation of twin boundary, etc. The initial setup consisted of three parts "patched" together in the simulation cell. Periodic boundary conditions were applied along the x-and z-directions; while the free boundary condition was applied along the y-direction. Figure 1 shows the simulation cell set up.
(a) (b) The simulation procedure included three stages: Setup of atoms, relaxation, and deformation, which was the same as other research work [31][32][33]37] on the boundary migration. During setup, atoms were arranged to produce crystallographic microstructure with two coherent ETBs existed between the matrix and twin blocks in cell. Initially, blocks of atoms in matrix arranged in perfect HCP crystal structures were with [1 21 0] direction along the x-direction and [1 011] direction along the z-directions, respectively. The atoms of twin in the middle block were arranged in extension twin orientation with the matrix. The normal direction of the twin boundary between the matrix and twin blocks was along the y-direction. The simulation volume was chosen so that the crystal lattice blocks matched up perfectly at the periodic boundaries along the x-and z-directions. The size of the simulation cell was about 6.4 × 32 × 20 nm 3 . The number of atom layers in x-, y-, and z-direction were 40, 169, and 108, respectively; the total number of atoms was about 180,000. The twin block was in the middle of the simulation cell with about 11.7 nm thick in the y-direction.
The relaxation and deformation procedure were performed using a constant number of atoms, volume, and temperature (NVT) velocity Verlet MD algorithm by the LAMMPS [39] package. A Nose-Hoover thermostat algorithm was used for temperature control. The temperatures in the relaxation and deformation procedure were controlled at 0.1 K. The time step used in the calculation was 1 fs. The simulation cell was firstly relaxed for about 20,000 time steps without any applied stress and strain. Then, as shown in Figure 1b, the atoms with about 1 nm thick layer at the left end were applied displacement along the z-direction; while the atoms with 1 nm thick layer at the right end were fixed during the whole simulation process. The simple shear strain was γyz = d/Ly, here, d and Ly were the loading-layer displacement and the length along y-axis, respectively. Three different strain rates were applied on the simulation cell: γ = 10 7 s −1 , γ = 5 × 10 7 s −1 and γ = 10 8 s −1 . Meanwhile, another simulation along the opposite direction was conducted on the same cell to investigate the influence of loading direction on the twin boundary migration. The microstructural characteristic and its evolution rule were visualized by the OVITO 3.0.0 software [40]. Atoms were colored by their orientation calculated by the polyhedral template matching method [41].

Results and Discussion
The stress-strain curves are plotted in Figure 2. Stress tensors of the system were derived from atom stress, which was calculated by the compute stress/atom command [42] on LAMMPS. As shown in Figure 2, it showed similar stress-strain behavior for both loading conditions. The stress initially The simulation procedure included three stages: Setup of atoms, relaxation, and deformation, which was the same as other research work [31][32][33]37] on the boundary migration. During setup, atoms were arranged to produce crystallographic microstructure with two coherent ETBs existed between the matrix and twin blocks in cell. Initially, blocks of atoms in matrix arranged in perfect HCP crystal structures were with [1210] direction along the x-direction and [1011] direction along the z-directions, respectively. The atoms of twin in the middle block were arranged in extension twin orientation with the matrix. The normal direction of the twin boundary between the matrix and twin blocks was along the y-direction. The simulation volume was chosen so that the crystal lattice blocks matched up perfectly at the periodic boundaries along the x-and z-directions. The size of the simulation cell was about 6.4 × 32 × 20 nm 3 . The number of atom layers in x-, y-, and z-direction were 40, 169, and 108, respectively; the total number of atoms was about 180,000. The twin block was in the middle of the simulation cell with about 11.7 nm thick in the y-direction.
The relaxation and deformation procedure were performed using a constant number of atoms, volume, and temperature (NVT) velocity Verlet MD algorithm by the LAMMPS [39] package. A Nose-Hoover thermostat algorithm was used for temperature control. The temperatures in the relaxation and deformation procedure were controlled at 0.1 K. The time step used in the calculation was 1 fs. The simulation cell was firstly relaxed for about 20,000 time steps without any applied stress and strain. Then, as shown in Figure 1b, the atoms with about 1 nm thick layer at the left end were applied displacement along the z-direction; while the atoms with 1 nm thick layer at the right end were fixed during the whole simulation process. The simple shear strain was γ yz = d/L y , here, d and L y were the loading-layer displacement and the length along y-axis, respectively. Three different strain rates were applied on the simulation cell: . γ = 10 7 s −1 , . γ = 5 × 10 7 s −1 and . γ = 10 8 s −1 . Meanwhile, another simulation along the opposite direction was conducted on the same cell to investigate the influence of loading direction on the twin boundary migration. The microstructural characteristic and its evolution rule were visualized by the OVITO 3.0.0 software [40]. Atoms were colored by their orientation calculated by the polyhedral template matching method [41].

Results and Discussion
The stress-strain curves are plotted in Figure 2. Stress tensors of the system were derived from atom stress, which was calculated by the compute stress/atom command [42] on LAMMPS. As shown in Figure 2, it showed similar stress-strain behavior for both loading conditions. The stress initially increased linearly during the elastic deformation for each curve, and then the stress fluctuated around an average value at all the strain rates. The atoms' configurations with the different shear directions under the shear strain 0.015 and shear rate 10 7 s −1 are plotted in Figure 3. It can be obviously seen that the twin grew up with shearing along the [1011] direction and contracted with shearing along the opposite direction. With the twin boundary migrating, the elastic energy was dissipated, leading to the stress fluctuating around a constant value. Like the dislocation motion, TB migration was also an energy release process, corresponding to the process jumping from a Peierls valley to the next with stress drops [43]. Referring to the computing method of CRSS in dislocations, the CRSS was calculated by averaging the stress peaks during the fluctuation process [43][44][45]. When it was shearing along the [1011] direction, the CRSS was about 34 MPa, 37 MPa, and 42 MPa, respectively, at strain rates of 10 7 , 5 × 10 7 , and 10 8 s −1 . For the opposite shear direction, the CRSS was 34 MPa, 37 MPa, and 41 MPa, respectively. The CRSS increased with the strain rates increasing, irrespective of the loading direction. increased linearly during the elastic deformation for each curve, and then the stress fluctuated around an average value at all the strain rates. The atoms' configurations with the different shear directions under the shear strain 0.015 and shear rate 10 7 s −1 are plotted in Figure 3. It can be obviously seen that the twin grew up with shearing along the [1 011] direction and contracted with shearing along the opposite direction. With the twin boundary migrating, the elastic energy was dissipated, leading to the stress fluctuating around a constant value. Like the dislocation motion, TB migration was also an energy release process, corresponding to the process jumping from a Peierls valley to the next with stress drops [43]. Referring to the computing method of CRSS in dislocations, the CRSS was calculated by averaging the stress peaks during the fluctuation process [43][44][45]    increased linearly during the elastic deformation for each curve, and then the stress fluctuated around an average value at all the strain rates. The atoms' configurations with the different shear directions under the shear strain 0.015 and shear rate 10 7 s −1 are plotted in Figure 3. It can be obviously seen that the twin grew up with shearing along the [1 011] direction and contracted with shearing along the opposite direction. With the twin boundary migrating, the elastic energy was dissipated, leading to the stress fluctuating around a constant value. Like the dislocation motion, TB migration was also an energy release process, corresponding to the process jumping from a Peierls valley to the next with stress drops [43]. Referring to the computing method of CRSS in dislocations, the CRSS was calculated by averaging the stress peaks during the fluctuation process [43][44][45].    By measuring the migration distance, we could get the proportion relationship between the loading displacement and TB migration distance. The ratios are presented in Table 1. They were approximately equal under all shear rates. It could be concluded that the TB migration distance was only related to the loading displacement and had nothing to do with the strain rate. In HCP metals, the theoretical {1012} twinning shear is given by v = c/a, where a and c are the crystal constant of material. As for Mg, it could be that v = 1.628, indicating that the theoretical twinning shear was about 0.124. The ratio in Table 1 was close to the twin shear, which meant that the loading displacement might be applied to accommodate the shearing deformation during the TB migration process. In addition, the atom strain tensor was calculated by the atom strain method [46] in Ovito. The atomic deformation gradient tensor F was calculated for each particle based on the particle displacement vectors, and then atomic Green-Lagrangian strain tensor E was derived for each particle from The E yz distribution curves along the y-axis are plotted in Figure 4, when the amount of γ yz was 0.015. The curves showed that the main strain was concentrated in TB migration region. The E yz was −0.0616 in the main migration region with shear along [1011] and it was 0.0626 with the opposite shear direction. The twin shear was calculated by s = 2E yz , they were about −0.123 and 0.125 in the dominating migrated region. By determining from the deflection of marker lines experimentally, Molodov [28] also got the amount of the twin shear of 0.126 in the twin region. The twin shears in simulation and experiment were both in excellent agreement with the theoretical value, which indicated that the shear was an indispensable part in the twin migration.
The TB migration mechanism under different shear directions were studied. The microstructural evolution of atoms were analyzed under the shear rate . γ = 10 7 s −1 ( Figure 5). The atoms in twin and matrix were marked in red and blue color, respectively. The atoms located on different {1210} planes were represented by the solid circles and solid squares, which were labeled as L 1 and L 2 , respectively. The twinning boundary on different {1210} planes were expressed by the black and red dotted line, respectively. By observing the microstructural characteristics in different layers (Figure 5a), the atoms on the twin boundary of L 1 layer were nearly along a straight line and in a mirror-symmetric twin relationship. Therefore, the twin boundary in this layer was recorded as a mirror-symmetric twin boundary (MSTB). Meanwhile, the atoms on the twinning boundary of L 2 layer were in a zigzag line and in a glide twin relationship [47,48], which was recorded as glide twin boundary (GTB). As shown in Figure 5b,d, the displacement vector arrows of atoms during the twinning boundary migration process were plotted. It could be found that the displacement evolution of the atoms in the two types of TB was different during the migration process.  The TB migration mechanism under different shear directions were studied. The microstructural evolution of atoms were analyzed under the shear rate γ = 10 7 s −1 ( Figure 5). The atoms in twin and matrix were marked in red and blue color, respectively. The atoms located on different {1 21 0} planes were represented by the solid circles and solid squares, which were labeled as L1 and L2, respectively. The twinning boundary on different {1 21 0} planes were expressed by the black and red dotted line, respectively. By observing the microstructural characteristics in different layers (Figure 5a), the atoms on the twin boundary of L1 layer were nearly along a straight line and in a mirror-symmetric twin relationship. Therefore, the twin boundary in this layer was recorded as a mirror-symmetric twin boundary (MSTB). Meanwhile, the atoms on the twinning boundary of L2 layer were in a zigzag line and in a glide twin relationship [47,48], which was recorded as glide twin boundary (GTB). As shown in Figure 5b,d, the displacement vector arrows of atoms during the twinning boundary migration process were plotted. It could be found that the displacement evolution of the atoms in the two types of TB was different during the migration process.   The TB migration mechanism under different shear directions were studied. The microstructural evolution of atoms were analyzed under the shear rate γ = 10 7 s −1 (Figure 5). The atoms in twin and matrix were marked in red and blue color, respectively. The atoms located on different {1 21 0} planes were represented by the solid circles and solid squares, which were labeled as L1 and L2, respectively. The twinning boundary on different {1 21 0} planes were expressed by the black and red dotted line, respectively. By observing the microstructural characteristics in different layers (Figure 5a), the atoms on the twin boundary of L1 layer were nearly along a straight line and in a mirror-symmetric twin relationship. Therefore, the twin boundary in this layer was recorded as a mirror-symmetric twin boundary (MSTB). Meanwhile, the atoms on the twinning boundary of L2 layer were in a zigzag line and in a glide twin relationship [47,48], which was recorded as glide twin boundary (GTB). As shown in Figure 5b,d, the displacement vector arrows of atoms during the twinning boundary migration process were plotted. It could be found that the displacement evolution of the atoms in the two types of TB was different during the migration process.  For a better illustration of the atoms' motion during migration process, the movement of six atoms, labeled in Figure 5a, in the twin boundary of L 1 layer and L 2 layer was tracked. In Figure 6a, it is shown that in the L 1 layer, the movements of the atoms A M , B M , and C M were nearly perpendicular to the normal of {1210} plane, which was transmitted from the matrix state (blue circle) to the TB state (green state) and then the twin state (red state). However, the movements of the atoms A G, B G , and C G in L 2 layer presented a different motion path (Figure 6b) it is shown that in the L1 layer, the movements of the atoms AM, BM, and CM were nearly perpendicular to the normal of {101 2} plane, which was transmitted from the matrix state (blue circle) to the TB state (green state) and then the twin state (red state). However, the movements of the atoms AG, BG, and CG in L2 layer presented a different motion path (Figure 6b). The atoms AG in the matrix state (blue square) firstly moved along the [4 043] T direction to the TB state (green square), and then moved along the [2 023] M into twin state (red square). Furthermore, the relative displacement vectors of atoms in different TB structure during the TB migration process under different shear directions is presented in Figure 7. The solid and hollow symbols represent the atoms' position before and after migration, respectively. Five atoms A, B, C, D, E were labeled to observe the atoms' motion during the TB migration process, in which atoms A, B, E were in L2 layer beside the GTB structure and atoms C, D were in L1 layer beside the MSTB structure. When the external shear strain was along [1 011] direction (Figure 7a Atom B, E moved to their proper position by the shuffling, their displacements were both smaller than bTD. As seen from the inset of Figure 7a, the displacement vectors of atoms C and D were smaller than bTD and could be resolved into two components. The component paralleled to the TB was both along [1 011] direction, whose value was approximate to its theoretical shear value. The component perpendicular to the TB of atoms C and D moved along opposite direction by shuffling movement to accommodate the microstructural transformation of the MSTB during the migration process. In Figure 7b, when shearing along [101 1 ], the movement was totally opposite to that shearing along [1 011], but the move mode of atoms A-E was in the same process. Furthermore, the relative displacement vectors of atoms in different TB structure during the TB migration process under different shear directions is presented in Figure 7. The solid and hollow symbols represent the atoms' position before and after migration, respectively. Five atoms A, B, C, D, E were labeled to observe the atoms' motion during the TB migration process, in which atoms A, B, E were in L 2 layer beside the GTB structure and atoms C, D were in L 1 layer beside the MSTB structure. When the external shear strain was along [1011] direction (Figure 7a), atoms A, B, E moved along [1011] direction, and the displacement value of atom A was approximately equal to the value of twin dislocation According to the moving characteristic on GTB, it could be concluded that Atom B, E moved to their proper position by the shuffling, their displacements were both smaller than b TD . As seen from the inset of Figure 7a, the displacement vectors of atoms C and D were smaller than b TD and could be resolved into two components. The component paralleled to the TB was both along [1011] direction, whose value was approximate to its theoretical shear value. The component perpendicular to the TB of atoms C and D moved along opposite direction by shuffling movement to accommodate the microstructural transformation of the MSTB during the migration process. In Figure 7b, when shearing along [1011], the movement was totally opposite to that shearing along [1011], but the move mode of atoms A-E was in the same process.
In the simulations, the magnitude of the shuffling was all smaller than b TD , which was quite different to the shuffling-dominated mechanism [27] proposed by Li et al. They thought the shuffling was much greater than the twin dislocation and the twin dislocation was too small to impact on TB migration. However, they only considered the direct transformation from the parent to the twin without thinking about the TB's contribution. Actually, the {1012} TB was not a single layer of {1012}, it contained several {1012} layers. As shown in Figure 6, the atoms' movement could not adjust immediately in a single step, the {1012} TB transformed gradually. Thus, the closer the {1012} plane to TB, the more similar the configuration to TB's, the smaller the shuffling magnitude would be. But in the process of TB migration, the shearing magnitude was constant, when TB migrated two {1012} layers, the shear displacement was equal to the b TD . That was to say, the transformation from parent to twin mentioned in shuffling-dominated mechanism was not the result of one slip of the TD, but twice, thrice, or even more times. It meant that the TB migration could not ignore the TD, on the contrary, the effect of shear caused by TD accounted for an important part. From the analysis of Figure 7 and the shear strain of the migration region, it was concluded that the TB migration was relative to both shear and atom shuffling. During TB migration, every atom was affected by the shear and different shear directions caused opposite move directions, leading to twinning or detwinning. The shuffling was only used for adjusting the GTB and MSTB structure evolution. In the simulations, the magnitude of the shuffling was all smaller than bTD, which was quite different to the shuffling-dominated mechanism [27] proposed by Li et al. They thought the shuffling was much greater than the twin dislocation and the twin dislocation was too small to impact on TB migration. However, they only considered the direct transformation from the parent to the twin without thinking about the TB's contribution. Actually, the {101 2} TB was not a single layer of {101 2}, it contained several {101 2} layers. As shown in Figure 6, the atoms' movement could not adjust immediately in a single step, the {101 2} TB transformed gradually. Thus, the closer the {101 2} plane to TB, the more similar the configuration to TB's, the smaller the shuffling magnitude would be. But in the process of TB migration, the shearing magnitude was constant, when TB migrated two {101 2} layers, the shear displacement was equal to the bTD. That was to say, the transformation from parent to twin mentioned in shuffling-dominated mechanism was not the result of one slip of the TD, but twice, thrice, or even more times. It meant that the TB migration could not ignore the TD, on the contrary, the effect of shear caused by TD accounted for an important part. From the analysis of Figure 7 and the shear strain of the migration region, it was concluded that the TB migration was relative to both shear and atom shuffling. During TB migration, every atom was affected by the shear and different shear directions caused opposite move directions, leading to twinning or detwinning. The shuffling was only used for adjusting the GTB and MSTB structure evolution.

Conclusions
In this paper, an Mg crystal with two coherent twin boundaries under different shear directions and strain rates was studied by molecular dynamics simulation. The results showed that the twin structure would be thickened under [1 011] shear direction and shrunk with shearing in the opposite direction. The CRSS of {101 2} TB migration increased with strain rate increasing. The ratios of TB migration distance and loading displacement ranged from 0.123 to 0.125 at all the strain rates, they were in excellent agreement with the twin shear value theory. According to the TB structure and atom displacement, twin dislocation was resolved by bTD → [4 043] M +[4 043] T . TB migration was achieved by both the shear and the atomic shuffling. Every atom was affected by the shear, different

Conclusions
In this paper, an Mg crystal with two coherent twin boundaries under different shear directions and strain rates was studied by molecular dynamics simulation. The results showed that the twin structure would be thickened under [1011] shear direction and shrunk with shearing in the opposite direction. The CRSS of {1012} TB migration increased with strain rate increasing. The ratios of TB migration distance and loading displacement ranged from 0.123 to 0.125 at all the strain rates, they were in excellent agreement with the twin shear value theory. According to the TB structure and atom displacement, twin dislocation was resolved by b TD → [4043] M + [4043] T . TB migration was achieved by both the shear and the atomic shuffling. Every atom was affected by the shear, different shear directions caused opposite atom move directions, leading to twinning or detwinning. The atom shuffling was only used for adjusting the GTB and MSTB structure evolution. Frankly speaking, the twinning boundary migration mechanism was discussed based on the coherent twin boundary. Whether it is suited for the incoherent TB would be a further research target in the subsequent works.