Integrated structural reconstruction of unit structures of the metastable grain boundaries in diamond-structured materials presents � rst-order like phase transition


 One of the important issues of studying grain boundaries (GBs) which has recently attracted increasing interests is to investigate the phase behavior of GBs that one GB with determined disorientation and plane orientation (known as macroscopic parameters) can exist as distinct phases and perform phase transition. While such an issue has been investigated in fcc and bcc metals, GB phases in other elemental materials have not been reported. This work by applying molecular dynamics (MD) simulation explored totally around 7000 meta-stable GB phases of the <110>∑9(22‾1‾) symmetric tilt GB of silicon, germanium and diamond carbon as diamond-structured elemental materials. Meta-stable phases commonly exist in different elements were discovered and some of them were successfully verified to be reasonable by first-principle simulation. The verified meta-stable GBs were subsequently proved to have different capability to transform to the ground-stable GB at elevated temperature under MD simulation and to perform different pre-melting behaviors. We discovered a bi-directional structural reconstruction mechanism of the unit structure belonging to one of the verified meta-stable phases, by which the unit structures can transform to identical unit structures of the ground-stable GB which can present ‘opposite orientation’. Through computing the kinetic barriers by nudged-elastic-band and annealing simulation using MD, the integral behavior of the unit structures’ reconstruction is found to be a first-order like phase transition. Our work extended the research on GB phases from metals to diamond-structured materials and discovered a new GB phase transition mechanism which has not been reported before.


Introduction
Grain boundaries (GBs) play a signi cant role in affecting the properties of many important polycrystalline materials 1 . Recent studies have reported that GBs consisting of elemental materials possess different phases or complexions 2,3 . Evidences for existence of these GB phases not only include indirect evidences from discovered discontinuous variation in properties of polycrystalline materials 4 and successfully explored GB phases of face-centered cubic (fcc) and body-centered cubic (bcc) metals by atomic simulation [5][6][7] . A direct observation of coexistence of binary GB phases of aluminium under atomic-scale microscopy was also achieved recently 8 .
Despite some discoveries in the phase behaviors of elemental fcc and bcc GBs, understanding of GB phase behavior of materials with other lattice structures is still poor. Considering that fcc and bcc metals are mostly of metallic bonds, it is of interest to extend such studies into materials with different nature of chemical bonds, such as those with covalent bonds. One of the signi cant species with this aspect is the diamond-structured materials and their GBs also have notable effects on performances of many functional devices like solar cells and transistors [9][10][11] . Compared with fcc and bcc GBs, few reported studies have applied the term 'GB phases' on diamond-structured GBs though some properties related to phase behaviors have been discussed. Firstly, computational work has suggested that silicon GBs have enormous meta-stable states 12 and their most stable structures can be obtained by reconstruction from those of certain meta-stable GBs 13 . Secondly, molecular dynamics simulations have shown that several different twist GBs of silicon whose structures were optimized at 0K by a molecular dynamic (MD) simulation 14 have different temperature-induced transitions, where some GBs disorder as the temperature raises while some can possess a highly ordered structure even at a high temperature close to pre-melting of the GB 15 . Considering these two aspects, there still exist interesting properties which have not been revealed yet. While previous studies by applying empirical atomic potentials have explored many GB states or phases of silicon at 0K 12,16 , only the reasonability of the most stable GBs were veri ed by rstprinciple computation. As for the explored meta-stable GBs as the by-product, they were often ignored and were thought to be unstable in a real material. However, there has been reported experiments verifying existence of meta-stable GBs. For example, one has observed a meta-stable GB containing a 5fold coordinated atom in a twinning GB of silicon 17 . Also, while the reconstruction from meta-stable GBs to the most stable GBs were predicted by comparing their atomic structures, rare previous studies have deeply investigated this reconstruction process. The integrated behavior of the localized reconstruction of each unit structures belonging to a two-dimensional GB still remains unknown. Therefore, it is of importance to investigate the reasonability of these meta-stable GBs and the transition process from them to the most stable GBs for the diamond-structured GBs.
In this work, we explored some structures at meta-stable states of the GB, which commonly exist in diamond-structured materials including silicon, germanium and carbon, by the molecular dynamics simulation, and veri ed their reasonability by applying density-functional-theory (DFT) simulation. We found that one of the explored meta-stable GB which has not been reported before has an interesting phase transition mechanism by a structural reconstruction of the well-de ned unit structures to become unit structures identical to the well-known GB with the lowest energy. This work discusses the dynamics of GB reconstruction process in the diamond-structured materials and the discovered performances of meta-stable GBs makes new insights in how to understand the importance of the enormous number of meta-stable GBs for diamond-structured materials.

Results
Explored diamond-structured GB phases at 0K GB states of silicon, germanium and carbon were explored by relaxing a series of non-identical initial con gurations obtained by tailoring a bicrystal model of the GB in molecular dynamics (MD) simulation at 0 K. This S9 GB was selected because it not only is a popular GB in polycrystallines 18 of silicon but also has its most stable structures well studied by simulation 13,16 . The tailoring operations include applying rigid body translation (RBT) on one grain shifting in the GB plane and deleting some atoms near the GB within a cutoff distance. All the relaxed GBs by this method are periodic and have the same periodicity as the two dimensional CSL of the GB. The atom deletion operations always eliminate complete atomic planes so the atomic density corresponding to a whole GB plane as discussed in previous studies on fcc and bcc GB phases 19 is constant. Although this is likely to exclude some GB states, we will show that the explored states here have generated important new insights in the phase behaviors of the diamond-structured GBs.
Previous studies on fcc and bcc GB phases assigned the explored GB states into several GB phases by a GB-property-density clustering method upon a series of GB properties [20][21][22] . Here we simply visualized the GB properties of the explored GBs including the GB energy (GBE), excess volume, GB stress and excess order parameters. The excess order parameters were normalized by the lattice parameter of silicon, germanium and diamond carbon to become dimensionless so that they are close for similar GBs of different elements. The de nition of these properties can be seen in the reference 23 . For each element, we screened out non-identical relaxed states and each of them has difference in at least one excess order of Q3, Q6 or Q7 larger than 0.02. The data of the properties and structures of all the non-identical GBs are provided in S1. Figure 1 illustrates the GB energy distribution of the relaxed non-identical GB states, which re ected an apparent diversity of GB meta-stable structure with intermediate energy. Figure 2 shows the results of all the computed GB properties. As can be seen in Figure 2, the excess order parameter distribution of silicon, germanium and carbon almost coincide, with that of germanium expands slightly out of the silicon and carbon regions. This means the meta-stable GBs of the three elements have similar structure though the germanium ones have slightly larger variance. On the other hand, despite the similar GB structures of different elements, the difference in properties of carbon GBs is much more distinctive than that of silicon and germanium GBs, which might derive from the difference in their atomic bond properties. Besides obtaining the statistical insights of the structure and property of meta-stable GBs, the physical signi cance of this result is that considering the complex ensembles of the meta-stable GBs, if these explored metastates are proven to be reasonable to exist, it predicts a complex performance of this GB under elevated temperatures.
Selected Meta-stable GBs with reasonability veri ed by density functional theory (DFT) simulation DFT simulation was conducted to verify the reasonability of some of the explored GB states. We veri ed the GB states of Si shown in Figure 1 (a) including the one with the lowest GBE (denominated as Ground GB) and three meta-stable states: two of them (denominated as Meta-1 and Meta-2 GBs) have all the atoms with 4-fold coordination and the other (denominated as Meta-3 GB) has a 6-fold atom per unit structure. Figure 3 shows the structures of the four selected GBs. While the Ground GB and Meta-1 GB has been reported in previous studies 13,16,24 , the Meta-2 and Meta-3 GBs have not been reported or discussed before. All the four GB states were also found in the explored germanium GBs by the MD simulation. As for the carbon GBs, the Ground, Meta-1, and Meta-3 GBs were explored by the MD simulation while the Meta-2 GB was not and thus not stable under the MD simulation. To verify the reasonability of all the four GB states of these three elements, all these four relaxed GBs by the MD of silicon, germanium and carbon (the Meta-2 GB of carbon was made manually by normalizing the size of Meta-2 silicon GB according to the lattice parameters) were input in the DFT simulation with the xed simulation cell to determine whether they are reliable structures. The comparison of the MD and DFT simulations of these four GBs for silicon, germanium and carbon is shown in Table 1. The DFT results veri ed the reasonability of most the selected meta-stable GBs including the Meta-3 GB of silicon and germanium with coordinate defects, though the Meta-3 GB is not reasonable for carbon which suggests that carbon is less likely to allow coordinate defects in the GB. The structures of all the GBs stabilized in the DFT simulation are provided in S2. Also, different from the MD results, the Meta-2 GB of Carbon were stabilized in the DFT simulation. Therefore, all the three GB states without coordinate defects were veri ed to be reasonable for Si, Ge and C. The relative phase stability of the four GBs were re ected by the GB free energy which is computed by quasi-harmonic approximation 25 as shown in Figure 4. The slopes of the curves present the difference in GB entropy S(GB state) where S(Ground) < S(meta2) < S(meta3) S(Meta-1). According to this computation, the stability of the Ground GB maintains to be the most stable phase up to a temperature of 1900K. Therefore, as what has been taken for granted in previous studies, the meta-stable GBs are non-stable and tend to transform to the Ground GB if enough kinetic energy is provided. All the GBs with computed GBE have identical structures by MD and DFT simulation; None means that the GB phase was not explored by the MD or not stable in the DFT simulation.
Structure variation of meta-stable GBs at elevated temperature (This title also a bit confusing) In order to understand the kinetic effects on the meta-stable to ground-stable transition, we investigated the structural variation of the meta-stable GBs at elevated temperature by heating the four selected GBs at different temperatures in MD simulation holding by 10 ps and computed their averaged excess order parameters of Q3, Q6 and Q7 within the subsequent 10 ps to quantify their structure similarity. As shown in Figure 5, by increasing the holding temperature, the similarity of the annealed structures originating from Meta-2 and Ground GBs is promoted, which was re ected by the converging curves of their excess order parameters which merges at about 1800K. The separate of their curves after 1800K suggests that their structures are not completely identical and ends up with different structural variations before premelting. On the other hand, Meta-1 and Meta-3 GB do not present such increase of similarity with the Ground GB by promoting temperature within the simulated time range. While for the Meta-3 GB with different atomic density, absent of transition to the Ground GB can attributed to the forbidden atomic diffusion due to the periodic boundary condition applied parallel to the GB plane, the fact that the Meta-1 with the same density behaves differently with Meta-2 under elevated temperature can be attributed to their difference in the free energy where Meta-2 has larger thermodynamic driving force to transform to the Ground GB; or to their difference kinetic barriers where, as will be seen latter, the unit structure of Meta-2 has a straightforward transition mechanism to become identical to the unit structure of the Ground GB which is not discovered in Meta-1. Furthermore, the Meta-1 and Meta-3 GBs present earlier decrease of excess order parameters which indicates an earlier pre-melting process than the other two GBs. This difference between the performances of meta-stable GBs at elevated temperature suggests the incompatible preferences of different meta-stable GBs to transform to the Ground GB and to pre-melt.
Structure reconstruction of the Meta-2 GB To understand the mechanism of transition from Meta-2 to Ground GB, we quenched the structure of Meta-2 GB which had been held at 1500K for 20 ps and relaxed it by MD. Figure 6 (a) shows the relaxed structure, where the two-dimensional GB structure comprises unit structures with three different states including one as the original unit structure of the Meta-2 GB (Fig. 6 (b)), and two which has been reconstructed to be identical to the unit structure of Ground GB where the two reconstructed states have certain difference in 'orientation'. The existence of differently oriented unit structures might be the reason that the annealed Meta-2 GB behaves differently with Ground GB before pre-melting although their unit structures are nearly identical after the complete transition of Meta-2 GB. According to Figure 6 (e), the reconstruction of each unit is re ected by rearrangement of bonds between two ' ag' atoms (blue and red balls) and their neighboring atoms. Coupling the two ag atoms as a 'dumbbell', the orientation of each transformed unit is determined by the direction of the rotation of this dumbbell during the reconstruction. Here, we symbolize each unit structure by a gray block plus a green block. In the gray block, the transition happens by rearranging the bonds of ag atoms of the dumbbell. A dumbbell lying horizontally in the z-direction gives a non-transformed Meta-2 GB unit which is represented by a 'minus' symbol in the gray block; and a dumbbell with ag atoms having different heights in the z-direction represents a unit transformed to the Ground GB which is symbolized by an 'up' or 'down' arrows in the gray block depending on whether the blue atom is higher or lower than the red atom. As the bonds in the green blocks are simply distorted during the transition, the green blocks can be thought as nearly non-changed structures. In this way, the two-dimensional GB can be thought as an 'tessellation' of chains of gray blocks with reconstruction separated by chains of green blocks without reconstruction.
Kinetic barrier of the structural reconstruction and the integrated behavior of unit structures To understand the integrated transition behavior of the whole two-dimensional GB by this reconstruction of each unit structures, it is of interest to investigate how the transition of a unit structure is affected by the states of its neighboring unit structures. Here we de ned some different transition routines and computed their barriers shown in Figure 7 by Nudged Elastic Band (NEB) implemented in LAMMPS. Process (1) can be thought as a kind of 'nucleation' process with one Ground unit appearing in a pure Meta-2 environment; process (2) & (5) have 'growth'-like characters generating one Ground unit adjacent to an existing Ground unit; process (3) and (4) can be regarded as certain 'mutation' which ends up with a Ground unit with opposite orientation to a nearby existing Ground unit. The transition barriers of process (1), (2) & (3) are almost identical, which indicates a weak interaction of neighbors along the z-direction due to the separation effects by the nearly non-changed green blocks. On the other hand, the transition barriers of process (1), (4) & (5) follows that (5) > (1) > (4), which illustrates two aspects of the transition behavior for units in a same gray chain. Firstly, once a nucleus is formed, the neighboring structure tends to transform to the same orientation as this nucleus. Secondly, the fact that the barrier of nucleation is larger than the barrier of growth predicts a typical rst-order like transition behavior. Note that the 'nucleation' and 'growth' de ned here of course have different aspects with those de ned in a typical phase transition process, where the nucleation barrier mainly comes from the interface energy and have thermodynamic barrier in the form of free energy. In this case, the neighboring structures are closely bonded with each other without generating an extra 'interface' defect between them and this nucleation & growth de nition relies on the fact that the reconstruction near an already transformed unit structure is easier to happen than to reconstruct a unit from a pure meta-stable environment.
To verify the prediction of rst-order like transition, a large Meta-2 GB of silicon containing 504 unit structures were heated at temperatures from 300K to 1600K (at each temperature holding for 2 ns) and the states of all the unit structures were monitored. Figure 8  the annealing simulation at low (< 700) and high ( 700K) temperatures. At T < 700K, the nucleation rate is promoted by increasing temperature providing higher energy to overcome the nucleation barrier. It also illustrates a simultaneous nucleation and growth process, which results from the fact that the barrier of growth is only slightly lower than that of nucleation. The steeper growth line than the nucleation line means higher growth rate than nucleation rate which is consistent with our kinetic barrier computation that growth has lower barrier than nucleation. At T 900K, the nucleation rate shortly reaches a maximum and then gradually decreases to an equilibrium value. This decrease derives from a merging process where a Meta-2 region separating two Ground regions with identical orientations is transformed as shown in Figure 8 (c). This merging process also limits the maximum nucleation rate to be less than 1/2 of the total number of unit structures which cannot be further promoted by increasing temperature. These results illustrate a rst-order like transition behavior happening in each gray block chain especially at high temperature where the nucleation process which nishes shortly is followed by a notable growth process to form grains with two different orientations.

Discussion
In this work, a series of meta-stable GBs of silicon, germanium and diamond-carbon were explored by molecular dynamics (MD) simulation, from which we selected the ground-state GB (Ground GB) and three meta-state GBs (Meta-1, Meta-2, and Meta-3 GBs) of silicon, germanium and diamond-carbon. The selected GBs were also veri ed to be stable in the density-functional-theory (DFT) simulation except one of the diamond-carbon with coordinate defects, which shows that some meta-stable GBs which were ignored or even thought to be non-reasonable with coordinate defects before are possible to exist. This work promoted the reliability of the assumption on the complexity for diamond-structured GBs at elevated temperature.
By the MD simulation at the elevated temperatures, we found that the unit structures of the Meta-2 GB are capable to reconstruct to become identical to the unit structures of Ground GB under the simulated temperature by MD in an acceptable simulating time range. On the other hand, such reconstruction is not discovered for Meta-1 and Meta-3 GBs and they are more likely to pre-melt at high temperature. This difference challenges the well-accepted convention to ignore the meta-stable GBs based on the assumption that they can always transform to the GB with lowest energy. By visualizing the unit structure of the annealed Meta-2 GB, the reconstruction of each unit structure was found to happen by the bond rearrangement of two ' ag atoms' in each unit structure and can result in two opposite orientations. By combining the results of kinetic barriers of reconstruction and MD simulation at elevated temperature, we presented a method to investigate the integrated behavior of a bunch of reconstructing unit structures, and obtained a rst-order like phase transition as an integrated behavior of localized reconstructions happening in separated chains of unit structures of the Meta-2 GB.
What is also important of our results is that we promote a new insight on how to understand the premelting of the GBs. In the previous study on the metastability of GBs, the meta-stable GBs explored at 0K were thought to be accessible at elevated temperatures. If this is true, we should expect that at certain high temperature with enough energy overcoming the kinetic barriers the meta-stable states with comparable energy but different structures can become simultaneously accessible so that the whole GB will become disordered. This disorder of structure might be the reason for the Meta-1 and Meta-3 GBs to pre-melt. In other words, these two GBs tend to pre-melt possibly because that their capability to transform to other multiple meta-stable GBs with higher energy is much greater than that to transform to the unique Ground GB. On the other hand, the Meta-2 GB can easily to transform to the Ground GB by a straightforward reconstruction mechanism at a relatively lower temperature before obtaining enough energy to jump the barriers to reach other meta-stable states. Once the transition is nished, it becomes close to the most stable Ground GB and no longer capable to transform to meta-stable GBs due to the large free-energy and kinetic barriers and this might be the reason that the Ground GB and Meta-2 GB has much less tendency to pre-melt.
The complex behaviors of the discussed meta-stable GBs here illustrate the importance to consider the kinetic barriers of transitions between different GB states and the entropy effects on their free energy variation at elevated temperature when discussing the ensemble properties of GB structures in real materials because we have shown that starting from different states of this GB annealing can generate quite different structure variation routines. Therefore, instead of taking it for granted that the meta-stable GBs can always transform to the ground-state GB, a more mindful consideration should be taken for subsequent studies on the diamond-structured GBs as well as for GBs of other materials which also possess large number of metastates. Also, this work presented a new possibility for the GB engineering once a meta-stable GB or a transformed GB obtained by annealing a special meta-stable GB with desired properties is discovered in the future.

Methods
Exploring GB states by MD GB states of the GB of diamond-structured materials including silicon, germanium and diamond-carbon were explored by MD simulation using LAMMPS 26 where the interaction between atoms were modeled with a modi ed Tersoff potential by Pun and Mishin for silicon 27 , and Tersoff potential for germanium 28 and a long-range bond-order potential for diamond-carbon 29 . The atomic structure of all the GBs were visualized by OVITO 30 . To obtain different initial con gurations to be relaxed by MD for GB states exploration, rigid body translations of the two crystals forming the GB were sampled with con nement in the cell of non-identical displacement (c.n.i.d) which is the minimum cell including all the non-identical translations in the GB plane 1 . The sampling grids are as ne so that each vector of one grid is around 0.03 of the lattice parameter of the diamond conventional lattice for each material (around a grid size of 0.106Å, 0.163Å, and 0.170Å for carbon, germanium and silicon GBs respectively). After each rigid body translation operation, atoms closer than a cutoff distance were deleted from one crystal and this crystal is subsequently shifted normal to the GB plane to avoid broken of interface by this atom deletion. For GBs of different materials, both the size of grids in CNID and cutoff distance to delete atoms are proportional to their lattice parameters to make the sampled initial structures of different materials nearly identical. For previous studies on bcc and fcc GBs [20][21][22] , varying the atomic density respect to a complete GB plane is necessary to explore GB phases while for diamond-structured GBs, many different GB phases can be explored without doing so.

DFT simulation
We conducted computations based on density-functional-theory (DFT) using projector augmented wave method (PAW) implemented in VASP 31,32 with a cutoff energy for the planewave basis set for silicon, germanium and carbon as 390eV, 250eV, 530eV respectively. One repetition of the GB unit structure was applied with K-point grids as one grid in the direction normal to the GB, two in the rotation axis direction and four in the third direction. We applied similar convention to that in the study on the bcc GB phases 22 that the shape and size of the simulation cell is xed and simply relax the position of atoms to make the GBE of different GBs comparable. The convergence criterion was that all the force on all the atoms were less than 0.05 eV/angstrom.

Kinetic barrier of localized reconstruction computed by nudged-elastic-band (NEB) in MD
Kinetic barrier of localized reconstruction of Meta-2 GB unit structures was computed by NEB in MD simulation with functions implemented in LAMMPS 33,34 . The process includes two routines. Firstly, 14 replicas from the initial structure to the nal structure were relaxed to nd the minimum-energy-path of transition by a force spring applied on each replica. Here, a parallel spring of 1 ev/angstrom 2 was applied to compute the parallel nudging force separating the replicas and an extra arti cial spring of 1 ev/angstrom 2 was applied in a direction perpendicular to the tangent direction to prevent the paths from forming acute kinks 35 . Secondly, the replica with highest energy climbs to the maximum point of this path to reach the saddle-point to obtain the transition barrier. The convergence criteria is that the force on each replica is lower than 0.04 eV/angstrom.

MD simulation at elevated temperature
The time step of the MD simulations at elevated temperature was 1fs. The canonical ensemble distribution at different temperatures were sampled by applying a Langevin thermostat. To compute the excess order parameter of the selected four GBs, the GBs were rstly relaxed at each temperature for 10 ps and then maintained for another 10 ps to compute the average excess order parameter 23 and the standard variation. To count the number of transformed unit structures and number of transformed nuclei for the Meta-2 GB during the heating process, the Meta-2 GB were heated at a series of temperatures from 300K to 1600K for 2ns. Each temperature was shortly reached within the initial about 1ps. The average differences of the z-direction coordinates of the two ' ag atoms' in intervals of every 1 ps were computed and the unit structure with absolute value of this average difference larger than 1.6 angstrom will be thought to have transform to the Ground state, and the sign of this average difference is used to judge the orientation of the transformed structure.    The four selected GB states. The color represent the potential energy of each silicon atom.

Figure 4
The GB free energy computed by quasi-harmonic approximation of the selected GB states.

Figure 5
Excess order parameters of the four selected GBs at elevated temperatures. Red frames indicate the temperature where the transition from Meta 2 to Ground GB is completed within the simulated annealing time.

Figure 6
Transition mechanism from the Meta-2 to Ground GB. (a) the GB structure of the Meta-2 GB which was quenched from 1500K to 0K and relaxed by MD; (b) the unit structure which has not transformed yet; (c), (d) the unit structure which transformed to the Ground GB state (left) and its identical unit structure in a Ground GB (right); (e) the side view of the GB core shown in (a), ag atoms represented by large red and blue balls, and the symbolized unit structure according to the relative positions of the ag atoms. xdirection is perpendicular to the GB plane; y-direction is along the rotation axis of this symmetric tilt GB. Except for the ag atoms, the color represents the atomic potential energy.

Figure 7
Kinetic barrier of reconstructions of a unit structure with its neighboring structures at different states.