Collective Atomic Displacements during Complex Phase Boundary Migration in Solid-Solid Phase Transformations

The A15 to bcc phase transition is simulated at the atomic scale based on an interatomic potential for molybdenum. The migration of the phase boundary proceeds via long-range collective displacements of entire groups of atoms across the interface. To capture the kinetics of these complex atomic rearrangements over extended time scales we use the adaptive kinetic Monte Carlo approach. An effective barrier of 0.5 eV is determined for the formation of each new bcc layer. This barrier is not associated with any particular atomistic process that governs the dynamics of the phase boundary migration. Instead, the effective layer transformation barrier represents a collective property of the complex potential energy surface.

Many properties of bulk materials are determined by internal interfaces.As only small amounts of interface active elements are required to modify the stability and mobility of interfaces, interface properties are a focus in alloy design.In Ni-based superalloys [1], for example, refractory elements such as Re, Mo, or W are added to suppress creep.These elements can induce the formation of topologically close-packed (TCP) phases [2] that are not coherent with the cubic and single-crystalline superalloy material.The TCP phases are detrimental to the mechanical properties of the alloys and thus their formation needs to be avoided or retarded.In addition to a detailed knowledge of the structure and stability of the interfaces between TCP phases and the host material it is key to obtain insight into the kinetics of the migration of the phase boundaries.
Simulating the kinetics of complex phase boundaries in solid-solid phase transformations up to experimental time scales remains one of the great challenges in the atomistic modeling of materials.Molecular dynamics (MD) simulations are limited to time scales that are orders of magnitude shorter than experimental ones.To observe phase boundary migration on such short time scales unrealistically high driving forces are required that can alter the underlying atomistic mechanisms.Another challenge is presented by the collective atomic displacements at the interface that are too complex to use a lattice based kinetic Monte Carlo [3,4] approach to follow the dynamics over extended time scales.To capture the kinetics of complex phase boundaries we have to go beyond standard atomistic simulation techniques.One possibility is to use accelerated MD [5] if a suitable bias potential can be defined (hyperdynamics), very large computational resources are available (parallel replica dynamics), or the important reaction rates can be assumed to follow an Arrhenius form to a high temperature where they can be observed directly with MD (temperature accelerated dynamics).
In this study we present the application of adaptive kinetic Monte Carlo (AKMC) [6] to interface migration between different phases.The AKMC approach allows for simulations of atomic systems over long time scales by focusing on the dynamics of rare events.The method remains flexible by identifying the atomistic processes on the fly during the simulation.This is crucial as the phase boundary migration involves long-range collective motions of atoms across the interface.While the AKMC approach is also used in conjunction with ab initio methods, for the simulation of interface migration these methods are still computationally too expensive and we therefore use an interatomic potential for the evaluation of energies and forces.Along the trajectory the system evolves through a complex energy landscape.The analysis of the topology of the energy landscape indicates the fundamental mechanisms that determine the kinetics of the moving phase boundary.
The system that we have studied is the phase boundary between the body-centered cubic (bcc) and A15 phase in molybdenum.The A15 phase is one of the TCP phases and has been observed as a metastable phase in elemental tungsten [7][8][9] as well as in transition metal alloys [10][11][12][13], in particular also in the context of superconductivity [14,15].The unary bcc=A15 interface provides a suitable initial model system for a complex phase boundary.The A15 phase has a cubic unit cell with eight atoms and two inequivalent lattice sites that are 12-and 14-fold coordinated, respectively.Energies and forces are calculated using an embedded atom method (EAM) potential for molybdenum [16].To check the validity of the EAM potential we performed density-functional theory calculations for the basic bulk properties of the bcc and A15 phase.The EAM results are in good agreement with the density-functional theory data (details are given in the Supplemental Material [17]).The supercell setup used in the AKMC simulations is shown in Fig. 1.For both phases the [001] direction is perpendicular to the interface and periodic boundary conditions are applied in all three dimensions.One interface is fixed during the simulation while the other can move freely.Our choice of cell, which is optimized for the bcc lattice constant, has a compression of the A15 phase by 7% resulting in a driving force for the transformation from the A15 to the bcc phase of 0.36 eV=atom.Despite the large energy gain we assume that the system fully thermalizes after each layer transition.The initial interface was optimized by removal of high energy atoms and subsequent relaxation of the atomic positions and the cell dimension in the z direction.
The AKMC simulations were performed using the EON code [26].In each kinetic Monte Carlo step the potential energy surface (PES) around the current minimum is explored to find all important escape processes and corresponding barriers.The barriers are used to determine the rate for each process within harmonic transition state theory.The corresponding prefactor can be calculated from the vibrational frequencies at the minimum and the saddle point.In the AKMC simulations we used a fixed prefactor of ν 0 ¼ 5 × 10 12 s −1 for all processes.We verified that the results are not affected by this approximation by comparing with simulations where the prefactor was calculated explicitly for each process.The identification of saddle points and respective atomistic processes on the highdimensional PES of this complex system is challenging and requires a robust approach that can find important saddle points with a high confidence.Here, we use hightemperature MD trajectories together with nudged elastic band [27,28] calculations (details are given in the Supplemental Material [17]).This allows us to assign a confidence to the completeness of our rate catalog at each step in the simulation [29].
In Fig. 2 the change in potential energy is plotted as a function of simulation time for a typical AKMC trajectory at T ¼ 300 K.As the simulation progresses the bcc phase grows and the energy drops.The steps in the energy profile result from layer-by-layer growth, where the energy decreases by 0.36 eV=atom for each bcc layer formed (13 eV for our 36 atom=layer simulation).The time scale of the simulation reaches up to microseconds for the entire A15 → bcc transformation to complete.Overall we have performed 20 simulations for each temperature with per layer transformation times ranging between 7 ps and 22 ns at T ¼ 600 K and 0.1 ns and 0.1 ms at T ¼ 300 K. Especially at low temperatures the time scales of the layer transition far exceed the range of classical MD.In between each layer transformation there is an energetic plateau consisting of a large number of states with similar energy connected by relatively small barriers.These states are similar in structure with an almost constant number of atoms in the bcc and A15 phase.They differ, however, in the detailed atomic structure of the interface.During the simulation the interface does not remain sharp; rather, a disordered layer is formed between the two phases.The bcc phase grows out of the disordered interface layer, which appears to be energetically more favorable than maintaining a sharp interface between the two phases.We observed the formation of the disordered interface in all our simulations, in the AKMC simulations as well as in the hightemperature MD simulations.This could be due to the relatively large mismatch and corresponding driving force and the disordered interface structure might be less pronounced or even absent for other cell geometries.Partially, the disordered interface may also be an artifact our simulations due to the short-ranged EAM potential.Changes in the atomic positions in this disordered interface layer do not have a significant influence on the total energy of the system, but they nevertheless constitute new system states and show how the system must explore the energy landscape to find a mechanism for the growth of the bcc phase.These states can be combined into so-called superbasins where transitions between superbasin states are fast.Computationally, these fast transitions significantly slow down the simulation.We therefore apply the Monte Carlo with absorbing Markov chains approach [30,31] to analytically determine an exact escape time from each superbasin.The superbasin states are essential for the overall dynamics and the underlying atomistic mechanisms of the transformation are nontrivial.
To analyze the AKMC trajectories in more detail we visualize the complex energy landscape that is explored during the transformation in a disconnectivity graph [32].An example for a trajectory at T ¼ 300 K is shown in Fig. 3.The disconnectivity graph in Fig. 3 does not represent the entire PES, but only a region that is visited over the course of a single layer transition.In the disconnectivity graph the minimum belonging to a particular basin on the PES can be identified together with the corresponding transition state energy at which the minimum can interconvert.It displays the topology of the energy landscape where the superbasins are clearly visible as groups of states.The processes connecting the different minima are not simple movements of a few atoms but complex rearrangements of entire groups of atoms at the interface, cf.Fig. 3.For each bcc layer to form the system needs to climb out of a superbasin to continue further down in the disconnectivity graph.
All layer transitions were not as simple as shown in Fig. 3.An example of a more complex transition corresponds to the energy plateau containing point b in Fig. 2.This transition is significantly slower than the other transitions for this AKMC run.Also notable is the larger potential energy gain at the end of the plateau of 26 eV compared to the 13 eV potential energy gain for the formation of a single bcc layer.The disconnectivity graph of this slow transition, in Fig. 4, illustrates that this transition is significantly more complex than the typical single-layer transition, in Fig. 3. Here, multiple barriers around ∼0.5 eV are crossed; these barriers are associated with forming bcc layers in a larger disordered interface.Once the system is able to escape the large disordered region, two layer transitions occur in rapid succession.We further quantify the phase transformation by the average time needed to form a new layer of the bcc phase.Using the common neighbor analysis [33] we locally define for each atom whether it is in a bcc or A15 environment.With this we can monitor the position of the interface.The bcc layer closest to the interface is identified and the time τ layer is recorded until the next full bcc layer has formed.The mean first-passage time τlayer is calculated as the average over 90 values and the rate for layer transformation is k layer ¼ 1=τ layer .The layer transformation rate is proportional to the interface velocity v, which is related to the mobility M and the driving force P via ¼ MP.At T ¼ 1000 K the velocity is 8.6 m=s, a typical value in MD simulations.At T ¼ 600 K the velocity is 0.14 m=s, which marks the lower limit of velocities that can be observed within MD simulations.At T ¼ 300 K the interface velocity is only 11 μm=s showing that with our AKMC simulations we can cover interface velocities that are in the range of experimentally observed ones.The interface mobility often shows an Arrhenius relationship and as shown in Fig. 5 we find such a behavior for the layer transformation rate The fit to the AKMC data (black dots in Fig. 5) results in an effective barrier for the layer transformation of ΔE layer ¼ 0.47 AE 0.07 eV (with a coefficient of determination R 2 ¼ 0.98).If we compare the effective barrier to the energy barriers of the individual atomistic processes, we do not find any particular type of process that is associated with a barrier of ΔE ≈ 0.5 eV.The transformation time is not dominated by any rate determining step, but the effective energy barrier extracted from the Arrhenius relationship is a characteristic of the energy landscape.It is correlated with the energetic depth of the superbasins, given by the energy difference between the lowest minimum and the transition state that must be overcome to leave the superbasin, cf.Fig. 3.At high temperatures of T ¼ 600-1000 K, the transformation times are short enough to compare with MD simulations [34] (shown in Fig. 5, red and blue dots).The fit to the MD data yields an effective barrier of ΔE layer ¼ 0.52 AE 0.03 eV (R 2 ¼ 0.97), which is in good agreement with the AKMC data.In Fig. 5 we also show the fit (solid line) to all data points (AKMC and MD), which results in a value of ΔE layer ¼ 0.50 AE 0.02 eV (R 2 ¼ 0.99).The good agreement between the low temperature AKMC and the high temperature MD activation barriers suggests that effectively the same parts of the energy landscape are explored so that the overall mechanism of the transformation is independent of temperature.For simulation cells with larger interface regions the effective barrier might change as we would expect a step nucleation and growth mechanism of bcc layers that is only beginning to show in the smaller cells used in this study.
In this Letter we have shown that we can determine the mechanism of complex phase boundary migration up to experimental time scales using the AKMC approach.The evolution of the phase boundary involves long-range  2. The mechanism has a larger disordered region than the characteristic mechanism in Fig. 3 and a higher barrier is overcome for the phase transition to progress.collective motions of groups of atoms and proceeds via a disordered interface layer.The effective activation barrier for the layer transformation and correspondingly for the mobility of the phase boundary is a characteristic of the complex energy landscape explored during the transformation.The progress of the phase boundary is therefore determined by collective features of the potential energy surface instead of a few simple atomistic mechanisms.

FIG. 1 .
FIG.1.Illustration of the phase boundary simulation geometry.(a) A 6 × 6 bcc phase and a 4 × 4 A15 phase were paired in the x-y direction with fixed cell lengths of 18.9 Å in the x and y direction, and 37.8 Å along the z direction.The black atoms represent fixed atoms, the blue atoms are the bcc phase, the red atoms are the A15 phase, and the gray atoms are in the interface.(b) Energy versus lattice constant L x ¼ L y , with fixed L z .The relaxed bcc phase is 0.06 eV=atom more stable than the A15 phase; our chosen supercell, near the equilibrium bcc lattice constant, has a driving force of 0.36 eV=atom favoring the bcc phase.

FIG. 2 .FIG. 3 .
FIG. 2. (a) Potential energy gained during an AKMC simulation at T ¼ 300 K and snapshots indicating the transformation mechanism from the A15 (red spheres) to the bcc phase (blue spheres).Gray spheres indicate atoms at the interface.(b) Enlarged portion of the energy plateau circled in (a) indicating a superbasin width of approximately 0.5 eV.C1 C2

FIG. 4 .
FIG.4.Disconnectivity graph of the slow process shown in Fig.2.The mechanism has a larger disordered region than the characteristic mechanism in Fig.3and a higher barrier is overcome for the phase transition to progress.

FIG. 5 .
FIG.5.Arrhenius plot of the formation of a bcc layer.An activation barrier of ΔE layer ¼ 0.50 AE 0.02 eV was found by fitting the MD (red and blue dots) and AKMC (black dots) data over the entire temperature range.The MD data were obtained with two different thermostats, an Andersen (red dots) and a Langevin (blue dots) thermostat.