Migration mechanisms of a faceted grain boundary

We report molecular dynamics simulations and their analysis for a mixed tilt and twist grain boundary vicinal to the 7 symmetric tilt boundary of the type {1 2 3} in aluminum. When minimized in energy at 0 K, a grain boundary of this type exhibits nanofacets that contain kinks. We observe that at higher temperatures of migration simulations, given extended annealing times, it is energetically favorable for these nanofacets to coalesce into a large terrace-facet structure. Therefore, we initiate the simulations from such a structure and study as a function of applied driving force and temperature how the boundary migrates. We find the migration of a faceted boundary can be described in terms of the flow of steps. The migration is dominated at lower driving force by the collective motion of the steps incorporated in the facet, and at higher driving forces by the step detachment from the terrace-facet junction and propagation of steps across the terraces. The velocity of steps on terraces is faster than their velocity when incorporated in the facet, and very much faster than the velocity of the facet profile itself, which is almost stationary. A simple kinetic Monte Carlo model matches the broad kinematic features revealed by the molecular dynamics. Since the mechanisms seem likely to be very general on kinked grain-boundary planes, the step-flow description is a promising approach to more quantitative modeling of general grain boundaries.

The thermodynamic driving force for facet formation and coarsening is reduction in the total GB free energy. It has been argued that GB surface stress acts as an opposing force to facet coarsening [11], thereby stabilizing facets with finite lengths. Hamilton et al. [7] challenged this idea by showing that in metals such as aluminum and for 3 {112} type facets, these interface stresses are very small and do not inhibit facet coarsening. Recent molecular dynamics (MD) studies [9] confirmed his theory and extended the investigation to GBs of various orientations. The authors of Ref. [9] concluded that facet coarsening is thermodynamically favorable for many GB systems, although it can be kinetically sluggish for some boundaries.
In our previous work, we applied MD to investigate the atomic-scale processes of migration in GBs with general plane orientations close to that of a symmetric tilt boundary [12], so-called vicinal GBs. Prior to migration simulations, * r.hadian@mpie.de we ran constant temperature MD simulations for up to 2 ns, to anneal the boundaries, which is a conventional practice. In these calculations, the boundaries all relaxed within one periodic length into nanofaceted structures on the order of 5-6 nm. However, our further analysis showed that over longer equilibration times, the nanofacets tended to coarsen and grow in size and formed a large two-facet structure with the maximum facet size allowed within one periodic length, while the plane inclination was constrained by the periodic boundary conditions to remain unchanged. We can anticipate that such a drastic transition in the equilibrated structure of the GB will also change its response to the presence of any driving force for GB migration, both in terms of its structure and its migration rate.
In the present paper, we study this effect for a particular GB with constant temperature MD simulations and find that, indeed, after a transient period, a new steady-state dynamic structure emerges in response to a constant driving force. This dynamic structure depends on the driving force and the temperature. This dependence is interpreted using geometric, continuum, and stochastic models.

A. Grain boundary construction
Most of the GBs we studied in Ref. [12] can be described as kinked boundaries. They had nanofaceted structures, which can be described as a system of steps, comprising straight segments joined at kinks, which can migrate over flat areas of the GB described as terraces, as seen in Fig. 1(a). These structures were shown to migrate by atomic shuffling at the kink sites, which, at the smallest scale, leads to the propagation of the kinks and showing a tendency to form a two-facet structure. (c) As the kinetics is too slow to form large facets, the low-energy faceted structure was manually set up.
hence to the propagation of the steps or nanofacets. We have chosen such a kinked GB as an initial structure for our further study reported here. We use an orthorhombic supercell, which defines Cartesian (X,Y,Z) axes, with the Z axis normal to the GB. We refer to the GB plane as the mean boundary plane (MBP). This plane is constrained by the periodic boundary conditions and does not deviate from its original orientation, whatever facets and steps it may display. The Miller indices of the MBP with reference to the lower grain are as follows: boundary normal lies along [− 7 16 29] in the Z direction and [− 22 23 −18] and [5 4 −1] form the in-plane X-Y directions.
(Capital X,Y,Z is used throughout the paper for the supercell coordinate system.) The GB plane is a mixed tilt-twist, one deviated by a small angle, θ = 5.87 • , from the symmetric tilt boundary plane, (−1 2 3).
The molecular statics and dynamics simulations with this initial structure were done using the Zope-Mishin interatomic potential for Al [13] in the parallel MD code LAMMPS [14]. Details of the GB creation and energy minimization steps are explained in Ref. [12]. In our previous work, we performed simulated annealing for up to 2 ns in an isothermal-isobaric ensemble (NPT ) prior to running migration simulations. Further investigation showed how an increase in the annealing time up to 10 ns results in two fully formed single terraces as shown in Fig. 1(b). Annealing simulations were performed at 400 K. In this figure and all other upcoming snapshots of GBs, bulk atoms are removed and only GB atoms are depicted, characterized by an order parameter ξ ord i that varies from 0 in the perfect crystal environment of the lower grain to 1 for atoms in the perfect crystal environment of the upper grain, where i labels an atom. The facet-coarsening process is associated with a small reduction in GB energy from 495 to 490 mJ/m 2 at T = 0 K. Facet coarsening was verified on a number of other kinked GBs vicinal to the same symmetric 7 tilt boundary. Additionally, the same GB structure of Fig. 1(a) was also annealed using a Cu interatomic potential [15], which revealed similar coarsening behavior.
The time for full thermodynamic equilibration to the large two-facet structure, which we refer to henceforth as a terrace-facet structure, can extend to tens or hundreds of nanoseconds. Therefore, we manually set up a terrace-facet structure as seen in Fig. 1(c), such that the terrace corresponds to the terrace {1 2 3} plane of the nanofacets and the facet to the {0 1 4} plane, which is topologically equivalent to a pileup of the steps that migrate across the terrace plane. As these planes are all coincident site lattice (CSL) planes, a triangle relationship of the following form holds: Accordingly, the terrace plus facet GB was constructed with a ratio of 7:2 repeat units on the corresponding planes. We note that an integer multiple of the same relative combination would produce the same MBP, with a corresponding scaling of the area of the supercell on this plane. The total energy of the system was then minimized by a sequence of rigid grain translations along y and the removal of atoms with anomalously short nearest-neighbour distances-in this case 30% of the lattice parameter-followed by static relaxation of the atomic positions. As a final preparation step, the structure, depicted in Fig. 1(c), was annealed for 2 ns at 400 K. This structure was further used as an initial structure for migration simulations.

B. Migration simulations
We ran migration simulations using a local driving force based on an artificial external potential proportional to the orientational order parameter ξ ord i [16], which resembles physical driving forces such as magnetic field anisotropy [17] or elastic anisotropy [18], at a number of temperatures and driving forces. All MD simulations were performed within an NPT ensemble. The Nosé-Hoover thermostat [19,20] and the Parrinello-Rahman barostat [21] with the pressure set to 0 were used for controlling temperature and pressure, respectively. Snapshots of migration were recorded every 5 ps and were relaxed. The GB velocity was obtained from the mean GB position, which was extracted from the recorded snapshots using the relationship where Z i is the normal component of a GB atom, the nonbulk atoms selected via LAMMPS centrosymmetry parameter, and ξ ord i is the aforementioned order parameter, which is 0 for one grain, 1 for the other, and ranges between 0 and 1 for the boundary atoms. The mean velocity v n is then determined by fitting a straight line to Z t versus time. Figure 2 depicts two snapshots of migration at 400 K, driven by a force of 3 meV/atom, at two different simulation times. In this figure, a periodic image of the supercell has been added for clarity. Evidently, new steps have detached from the lower edge of the facet and are migrating along the terrace. The equilibrated terrace-facet structure is seen to be unstable in the presence of an applied force of this magnitude. These steps will go on to cross the terrace and join to the next facet. Depending on the temperature and applied driving force, the structure evolves differently; specifically, the number of detached steps and their spacing varies. Hence it was necessary to develop a step-detection method to describe the evolution of the terrace-facet structure in terms of its propagating steps. We can then see how, under a sufficiently high-driving force, the boundary evolves from the fully faceted state to a quasisteady state, characterised by a train of migrating steps on the terrace.

C. Step detection
To identify the step structure, it is convenient to rotate the Cartesian coordinate frame about the Y axis, taking the {X,Y } plane from the MBP to the plane of the terraces, as depicted in Fig. 3(a) by the red Cartesian axes. We refer to this choice of axes as the terrace frame. As the steps all have the same height h, they can be located and categorized in the terrace frame. The line traces are determined from binning the atomic positions in each step category [ Fig. 3(b)] in the y direction and recording the position of atoms with the lowest x as seen in Fig. 3

(c).
By recording individual step traces in the {x,y} plane as a function of time, we have all the data we need to monitor and reconstruct the entire migration history of the GB in three dimensions. It is worth noting that, in addition to their height, these steps are characterized by a small Burgers vector of the displacement shift lattice that was identified using a defect analysis method [22] to be of the type: 1/10 1 1 1 . Thus the step regarded as a line defect is more precisely a disconnection [23,24] of which we only consider here the part played by the step component. The role, if any, of the small Burgers vector in migration of the steps, is a subject of ongoing investigation.
For each step line trace, the average along the y direction is recorded at regular time intervals, and this average is what we refer to as the position of the step. Figures 4(a), 4(b) and 4(c) show the position of an individual step, marked by a thick line, as the simulation advances. The position of this individual step as a function of time is plotted in Fig. 4(d). At the end of a run when the position of all steps are plotted, a characteristic step/time graph is produced as shown in Fig. 4(e). This final graph provides us with an overview of the evolution of the step distribution projected onto the x axis during the course of a simulation.

A. General observations
The simulations all started from the periodically repeated system of a single terrace plane (−1 2 3). This terrace is terminated at each end by a facet of orientation (0 1 4) (2 2 3), which is formed of eight closely packed steps on the terrace plane. The total number of steps in the supercell is fixed at eight by the orientation and supercell size. At any given time, a vertical line drawn on Fig. 4(e) crosses all eight steps. Figure 4(d) identifies different stages of the motion of each step, in which the step begins to detach from the bottom of the facet with a frequency f d , migrate across the terrace at a velocity we shall refer to as v 1 , and as it arrives at the top of a facet, its velocity decelerates to v 2 , and at this slower velocity it moves down through the facet.
The history of all the steps is similar and can be followed in the step/time graph of Fig. 4(e). The average values of all step properties can be obtained from this graph. f d is determined from the average number of detached steps per unit time and the two characteristic velocities of v 1 and v 2 are loosely defined by the two sets of slopes of the corresponding lines on the step/time graph at stages (2) and (3) of the step motion, respectively. There are, of course, local fluctuations in these velocities.  (2) corresponds to its fast progress on the terrace plane, and stage (3) corresponds to its motion when it has arrived at the facet in front. The graph in (e), which depicts the average position vs time for all of the steps during the course of the simulation, is referred to in the text as a characteristic step/time graph. This particular simulation was done at 400 K, using a driving force of 3 meV/atom. At any given temperature and driving force, after a short transient period the system reaches an approximate steady state, in which the number of steps migrating across the terrace and the number remaining in the facet is, on average, unchanged. One extreme in behavior was observed at the lowest

B. A geometrical model
In this section, we make the above qualitative description of the processes more quantitative. Most importantly, we relate the step and facet migration rates to the steady-state migration rate of the GB plane, v n , which is the time-averaged normal velocity of the MBP.
As the steps move laterally, the MBP progresses upward. The schematic Fig. 6 shows the relationship between the three stages of the step motion and their corresponding rates extracted from the terrace frame in relation to the average velocity of the MBP as a single step progresses. A full period of motion for each single step of the height h includes detachment with the frequency of f d , followed by a propagation with the velocity v 1 to reach the upper terrace-facet corner, and then a collective progression with a velocity of v 2 once the steps are on the facet until they reach a full period of length: L cos θ , in the terrace frame. The height of the step shown in Fig. 6 has been magnified for clarity.
It is clearly visible at all simulation temperatures and driving forces that v 1 > v 2 . From the fluctuation-dissipation theorem, which has been successfully applied to mean GB profiles [25,26], we expect that in the absence of a driving force the mean-squared displacement of a position variable over some time-in this case, the position of a step-will be proportional to its mobility when driven. A step that is within a facet is restricted in its displacements by the proximity of its close neighbours, so it is to be expected that its velocity in response to a driving force will be less than that of a free step on a terrace. Indeed, we can study a GB consisting only of a facet orientation, in this case (0 1 4) (2 2 3), to calculate v 2 in isolation. We have done such a calculation and our results confirm this notion.
Relative to the velocities v 1 and v 2 , the facet itself always appears to migrate much more slowly. The facet velocity, which we will refer to as v f , indicated by the slope of the dashed lines in Fig. 7, remains almost zero in the terrace frame. This implies that the facet during the course of migration is almost stationary in the MBP frame as well because the difference in the facet velocity in the two frames is proportional to sin (θ ) and in our case is very small.
As it can be seen in Fig. 6, there are two distinct geometrical contributions to the upward movement of the boundary. The terrace is moving upward with every passing step at a velocity of It is worth noting that the propagation velocity of the steps on the terrace v 1 does not play a role in the average upward velocity of the GB. The second contribution comes from the facet that is moving due to the motion of the steps at v 2 within it. Figure 6 should make clear, if v 2 were vanishingly small, the facet would be moving backward along the terrace plane at a velocity f d W , due to the steps detaching from its bottom corner and reattaching at its top corner, where W is the width of one single step in the mean boundary coordinate system: If the facet is, on average, stationary on the MBP, v f ≈ 0, this expression simplifies to v n = f d h/ cos θ.
On the other hand, if it is exactly stationary on the terrace plane, the corresponding result is v n = f d h cos θ , differing from Eq. (7) by only 1.05%, which is within the error of our estimates. Figure 8 verifies the agreement between the MD velocities extracted directly from the simulations using Eq. (2) and the geometric prediction from Eq. (6). At the low driving force of 0.5 meV/atom depicted in Fig. 5(a), we observed a limiting behavior, i.e., no step detachment has occurred in our simulation time of 15 ns, hence f d ≈ 0. Then the collective motion of the steps in the facet with the velocity v 2 entirely determines v n . For higher applied forces, v n can be obtained from the simplified Eq. (7). The physical reason why the facet velocities in these cases are so small in magnitude, i.e., why f d and v 2 are balanced, even at different temperatures and driving forces, is unclear, and will require an explicit dynamic model of the interacting steps in future studies. The step detachment rate f d , which was the ratedetermining factor in the majority of our MD cases, had not been anticipated in our previous work due to the incomplete thermal equilibration. Our observation shows that the overall migration rates of our terrace-facet structure drops by at least a factor of 5 at higher driving forces of 3 and 5 meV and by about 2 orders of magnitude for the smallest driving force of 0.5 meV. Figure 8 implies a deviation from the Arrhenius behavior for the terrace-facet structure as opposed to our previous nanofaceted boundary. We would need a denser temperature/driving-force mesh to characterize the source of this deviation. To run long simulations to ensure that the steady state has been reached for a dense set of temperatures and driving forces would be computationally too costly.
In the following section, we illustrate the above formulas within a scenario in which the three variables are stochastic.

C. A one-dimensional stochastic model
We employ a one-dimensional (1D) kinetic Monte Carlo (kMC) algorithm to simulate the possible steady-state conditions numerically. Figure 9 demonstrates schematically how forward-moving walkers on a 1D periodic lattice can be represented to model the motion of a terrace-facet structure. In this 1D picture, when the site at position x is occupied, denoted by σ x = 1, a step is present, and when the site is empty, σ x = 0, there is no step. Figure 9 shows how this 1D construction maps to the motion and pileup of steps in our terrace-facet structure. Similar discrete lattice models, with explicit rules for movement in one dimension, have been used in the simulations of traffic flow [27].
A discrete step forward is represented by Details of the rules and the kMC algorithm are given in the Appendix. We have used the step/time graphs from our MD results at a driving force of 3 meV/atom to estimate the transition rates r i for f d and v 1 using: which is linear in the applied driving force ε; A i is the prefactor, and E i is an effective energy barrier. This model assumes that the energy barrier and prefactor for the two rates are independent of the driving force. Our fitted expressions, although a rough estimate, suffice to demonstrate the theoretical steady-state conditions of Sec. III B. Alternative models of step flow could be constructed based on interactions between the steps, cellular automata, or hydrodynamics. With this model, we have simulated the steady-state behavior under the described conditions in Sec. III B. The motion of steps was recorded and plotted as step/time graphs as we used previously to represent the MD results. The results are shown in Fig. 10 at 3 meV/atom driving force and T = 400 K. Figures 10(a) and 10(c) represent the limiting conditions in which f d and v 2 are vanishingly small, respectively. The direction of motion of the facet is seen to be as predicted by the previous geometric considerations. Figure 10(b), in which the steps move at a velocity of v 2 = f d W , marked by the white arrow, while the facet remains stationary, marked by dashed lines, is the closest to the observation in our MD results. Note that in the 1D lattice, the v 2 cos (θ ) that was used in the geometric model to refer to the velocity of the collective motion steps is reduced to v 2 .
With the parameters appropriate to a stationary facet, we then ran the kMC model to generate step/time graphs at the low driving-force limit of 0.5 meV at 400 K. It is worth noting that in the kMC model, if we fix the temperature, the only effect of varying the driving force is to scale the absolute rates of the events, not their relative rates, as the rates are all linear in driving force according to Eq. (9). Thus Figs. 5(b) and 10(b) are essentially the same model, apart from the stochastic variations, with different time scaling. Figure 5 shows a comparison between MD results and the kMC prediction. The kMC results in Figs. 10(b) and 5(b) are in broad agreement with the two corresponding MD simulations, but the interesting physics lies in their differences.
The stochastic variation in the detachments seen in the MD is much less than that in the kMC, suggestive of a correlation between the arrival event of a step at the lower corner of the facet with the detachment event there. On the other hand, in the MD, there is variation in the individual velocities of steps within a facet v 2 , neglected in the kMC model, and evidence of variations on a longer timescale than represented in the 043601-6 kMC. We can speculate that these features are related to the complexity of the migration process in the MD. Complexities such as the dynamics of kink propagation along steps and the interactions between steps, which are not captured in the more simple geometric and kMC models-although the mentioned models have captured the essential physics.
Unlike its MD counterpart, the kMC simulation in Fig. 5(b) predicts detachment events for steps, which is an indication that the energy barriers for the detachment process at low driving force are higher than predicted by our assumption of linear extrapolation.

IV. CONCLUSIONS
We have investigated how an asymmetric GB of mixed tilt and twist character equilibrates its structure and migrates in response to an external driving force. Molecular static simulations show that the initially kinked GB is relatively high in energy. It could achieve a lower potential energy by forming nanofacets that would coarsen with time. Since this process proceeds very slowly to completion on the timescale of MD, we have manually set up the final structure, comprising a single terrace plus facet within a large periodic supercell with the specified mean tilt-twist inclination. Starting from this structure, we have performed MD simulations over a range of temperatures and with various external driving forces for migration. In the presence of a driving force, the fully faceted structure is observed to move by releasing steps from the terrace-facet junction and/or by the collective motion of steps on its facet. To trace the motion of the steps, we have developed an automatic step-detection method. This method proved useful for describing the migration mechanisms of the GB solely on the length scale of its steps, without having to reproduce the atomistic level of detail. By this means, we showed how the steps are conserved as they move from free motion on the terrace to flowing more slowly as a collective bunch of steps on the facet.
By constructing graphs of step position versus time, we could identify three stages of the GB step flow, namely detachment of the steps, their propagation on the terrace and their collective motion within a facet. At higher driving forces, above 0.5 meV/atom, the rate-controlling mechanism appears to be the step detachment from the facet junction. At lower driving forces, our MD results suggest that the relative barrier to step detachment increases and the boundary moves only by the collective motion of the steps on the facet.
We used a geometrical model to show the relationship between the characteristics of lateral step flow and the forward motion of the MBP. We introduced, in addition, a 1D kMC model, which reproduces the steady-state conditions of motion exhibited by MD. A detailed comparison of these models with the MD simulations reveals that the step-flow parameters, that is, the step-detachment frequency and the velocity of the steps in the facet, are highly correlated. This correlation indicates that the motion of a step is affected by step-step interactions.
On the limited timescale of a MD run, initial GB structures do not fully equilibrate to enable us to simulate real processes at GBs. The mechanisms of migration we have described here could only have been discovered by initially preparing the boundary in a structure close to being equilibrated, in which just two inclinations, namely the terrace and the facet, occupy the supercell. As shown in a previous study [12], when starting from the same orientation but with a structure comprising only well-separated steps, without any bunched together to form a facet, we could identify only the propagation stage of the free steps on the terrace, which occurs via the atomic shuffling around the kink sites with a migration barrier of 0.1 eV. For the fully faceted and energetically more stable initial structure of this work, the steps need to detach from the facet they represent before they can propagate on the terrace. As a consequence, the effective energy barrier to boundary migration increases. Indeed, at the lowest driving forces, where, in the present study, the fully faceted boundary is seen to move only via the collective motion of steps on the facet, the barrier to migration might even be higher. Thus, to explore with MD the migration processes likely to be important in a real GB, having decided upon a mean inclination for the boundary, we had to take care to set up an initial structure close to the equilibrium one before running the long MD simulations.
Faceted GBs are ubiquitous in metals, and the step-flow processes described in this paper are likely to be very general. Our modeling and analysis of the migration of a mixed tilt and twist GB has demonstrated quantitatively the importance of step release from facets and the propagation of steps along terraces. The insights gained from this work suggest that a dynamic model for the mobility of steps can be developed further to include an explicit interaction between them. We hope this will be the basis of a versatile and quantitative tool for simulating the migration of a wide range of general GBs.

ACKNOWLEDGMENTS
R.H. was funded by an Alexander von Humboldt postdoctoral fellowship and would like to thank A. Stukowski for fruitful discussions. M.W.F. is grateful to the Alexander von Humboldt-Stiftung for an award that enabled his extended stay in Düsseldorf and consequent collaboration on this project. This project has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 639211).

APPENDIX
As the first configuration, we choose M steps (occupied sites) on L lattice sites such that the ratio M/L matches that of our MD simulation cell, depicted in Fig. 9(a). We then assume that three types of events can occur independently at each simulation time step. The steps can detach, or propagate on the terrace, or move collectively in the facet. The following rules describe the possible configurations in our 1D cell: (1) If N M/2 + 1, where N is the number of contiguous steps, these N steps are defined as a facet and are moved collectively with a probability corresponding to the velocity v 2 .
(2) The rightmost step of a facet can detach with a probability corresponding to the rate f d .  (3) If N M/2 + 1, these steps are considered as free and they are moved with a probability corresponding to the velocity v 1 .
At each time step, all the possible events and their corresponding transition rates r i are enumerated, and a cumulative rate function R = n i=1 r i is constructed, where n is the current total number of possible events.
Following the standard kMC algorithm, a random number is generated from the uniform distribution, ∈ (0,1), which decides which event i will take place according to After an event is executed, the time step is updated with t 2 = t 1 + t, where t = R −1 ln 1/ , where is another random number from the uniform distribution.