Computing the crystal growth rate by the interface pinning method

Crystal growth is of paramount importance in many branches of condensed matter physics.^1,2 An important parameter in phase field equations³ describing crystal growth is the kinetic coefficient defined as the proportionality constant between supercooling (or superheating) and the average interface growth (or melting) velocity $〈 x ̇ s 〉$⁠.^4,5 We define the kinetic coefficient M using the difference in chemical potential between the two phases μ_sl = μ_s − μ_l as a measure of supercooling (or superheating),

In the spirit of the fluctuation-dissipation theorem,⁶ we suggest to learn about the interface dynamics by investigating spontaneous fluctuations when a bias potential is added to the Hamiltonian. Specifically, we determine M indirectly in a computation where trajectories are pinned to two-phase configurations (Fig. 1) by adding a harmonic energy coupling to an order-parameter of crystallinity.^7,8 In effect, non-equilibrium crystal growth is converted into a well-defined equilibrium problem. We devise a stochastic model of fluctuations that assume Smoluchowski dynamics for crystal growth and show that the growth rate can be inferred from the terminal exponential relaxation of the order-parameter. The chemical potential difference between the solid and the liquid μ_sl is known from the average force exerted by the bias field on the system, and thus the proportionality constant M can be computed in a single equilibrium simulation.

The paper is organized as follows: Sec. II is a brief introduction to the interface pinning method. In Sec. III, we motivate and give the solution to a stochastic model of thermal fluctuations with and without a pinning potential. In Secs. IV and V, we compute crystal growth rates for the Lennard-Jones (LJ) model and combine the method with ab initio density functional theory (DFT) to compute growth rates of the elements Na and Si. The paper is finalized with a discussion of the method (Sec. VI). The Appendix gives a detailed analysis of our stochastic model.

In the following, we give a brief review of the recently proposed interface pinning method^7,8 for studying the solid-liquid phase transition. Consider a system of N particles in an elongated periodic simulation cell as sketched in Fig. 1. For the configuration R = (r₁,  r₂, …,  r_N), let the system have the Hamiltonian $H 0 ( R , R ̇ )$⁠. Hydrostatic pressure is ensured by allowing the box length in the z-direction to change and constant temperature is ensured by connecting the momenta of the particles to a heat bath. We refer to this as the Np_zT-ensemble. When the initial R is a two-phase configuration, the thermodynamically stable phase will grow on the expense of the other phase. In an interface pinning computation, this is avoided by adding an auxiliary spring-like energy term with a spring constant κ and an anchor point $Q ̄$ to the Hamiltonian so that R is biased towards two-phase configurations,

Here, Q(R) is a measure of crystallinity of the system. In practice, we use long-range order as measured by the magnitude of the collective density field |ρ(k)| (Ref. 8) (unless otherwise stated)

where k is a wavevector that corresponds to a Bragg peak. The averaged force α the system exerts on Q is proportional to the chemical potential difference between the two phases

Here Q_sl = Q_s − Q_l and Q_s and Q_l are the mean values of the order parameter when the system is entirely crystalline or liquid, respectively. When the system is in equilibrium with respect to $H$⁠, the relative position of the interface stops moving up to thermal fluctuations and the force α is balanced by the average force $κ ( 〈 Q 〉 − Q ̄ )$ of the spring-like bias-field. Thus, μ_sl can be computed from the average spring force as follows:

Since, we have used a harmonic pinning potential, the distribution of the order parameter will be Gaussian with variance k_BT/κ. The time evolution of the order parameter Q(t) depends on the trajectory R(t) that itself is determined by $H 0$⁠. The main idea of the method suggested in this paper is that if we can understand the time fluctuations of Q(t) when the interface pinning field is applied, we can deduce how the system evolves in the absence of the interface pinning potential. In Sec. III, we devise a stochastic model that assumes that crystal growth is an over-damped Smoluchowski process. From this model, we deduce the value of the kinetic growth coefficient M.

This section suggests a simple model for the dynamics of Q(t) that enables us to determine the value of M from an interface pinning simulation. Assuming Newtonian dynamics of R, the Q(t)-trajectory is deterministic and uniquely defined from the initial values of R and $R ̇$⁠. We are, however, only interested in the statistical behavior of Q(t), e.g., the autocorrelation function 〈ΔQ(0)ΔQ(t)〉. To this aim, we devise a stochastic model with an effective Hamiltonian that is a function of coarse-grained coordinates. The equations of motion of these coordinates include Langevin noise forces representing degrees of freedom that have been neglected.⁹ The essential model parameter for determining M is “the friction constant of crystal growth” γ. We determine the value of this parameter by fitting the model to simulation results.

As mentioned earlier, the interface pinning method requires the definition of an order parameter Q(R) that quantifies how much of the system is crystalline at the configuration R, which gives the positions of all atoms. Ideally, the order parameter Q should be proportional to the number N_c of crystalline particles such that any increase in Q corresponds to a growth of the crystalline region. In practice, however, the crystallinity of the system is measured using the Fourier transform of the density field as order parameter, Eq. (3). In this case, Q changes not only due a decrease or increase of the number of crystalline particles but also due to other spatial fluctuations such as lattice vibrations. Thus, we write the order parameter Q as a sum of two contributions,

where q is assumed to be proportional to the number of crystalline atoms and f takes into account fluctuations that do not change the number of crystalline atoms. The argument t emphasizes that all three variables Q, q, and f evolve in time along the trajectory R(t) of the system.

In an interface pinning simulation, the order parameter Q is computed explicitly, while the variable q, which contains the relevant information about the growth of the crystalline region and the motion of the interface, is not directly accessible. To extract information about the interface dynamics from the simulation, we, therefore, construct a simplified two-dimensional model that separates the effect of q and f on the dynamics of Q. In this model, q and f are the sole dynamical variables. We postulate that their dynamics is governed by the effective Hamiltonian

The first term on the right hand side, αq, arises from the imbalance in chemical potential between solid and liquid, which drives the growth (or decrease) of the crystalline region. This term drives the interface motion in the absence of the pinning potential (κ = 0). The effective force α is proportional to the difference in chemical potential, μ_sl. The proportionality constant depends on the choice of crystallinity order parameter, Eq. (4). The second term, $( κ / 2 ) ( f + q − Q ̄ ) 2$⁠, couples to Q = q + f and holds the interface fixed near a given position. The parameters κ and $Q ̄$ appearing in the pinning potential denote the force constant and the position of the bias, respectively. The third term, (κ_f/2) f², reflects the harmonic potential energy experienced by phonon fluctuations. Finally, the last term is the kinetic energy associated to the time derivative $f ̇$ of the phonon degree of freedom f, which carries the effective mass m_f. A kinetic energy term for the coordinate q is not necessary, because the dynamics of q is assumed to be overdamped.

Based on the effective Hamiltonian of Eq. (7), we next devise stochastic equations of motion to describe the coupled time evolution of q and f. Since, the growth of the interface separating the liquid from the crystalline phase is slow on the molecular time scale, inertial effects in the motion of the variable q are negligible, suggesting to describe its dynamics with an overdamped equation. In contrast, lattice vibrations are fast, which requires an underdamped description for the dynamics of f. Thus, on a coarse-grained level, we expect the dynamics of q and f to be governed by the following pair of coupled Langevin equations,

Here, γ and γ_f are friction constants associated with q and f, respectively, and η_q(t) and η_f(t) are δ-correlated Gaussian random forces. Both the friction and random forces reflect the degrees of freedom neglected in the simplified model.¹⁰ The variance of the random forces and the friction constants are related by the fluctuation-dissipation theorem,⁶ 〈η_q(0) η_q(t)〉 = 2γk_BTδ(t) and 〈η_f(0) η_f(t)〉 = 2γ_fk_BTδ(t), where k_B is the Boltzmann constant, T is the temperature, and the angular brackets 〈⋯〉 indicate a thermal average. Solving the above Langevin equations yields q(t) and f(t) from which the time evolution of Q(t) = q(t) + f(t) can be computed.

To apply the model to simulation data, we first consider the autocorrelation function 〈ΔQ(0)ΔQ(t)〉, where ΔQ(t) = Q(t) − 〈Q〉. As shown in details in the Appendix, this autocorrelation function is a sum of three complex exponentials with rather complicated arguments. To obtain a practical expression, we consider the long-time limit of the model and utilize the “separation of time scales”-assumption γ(1/κ + 1/κ_f) ≫ γ_f/κ_f,

where A = κ_f/(κ + κ_f) is the strength and $τ=γ ( κ − 1 + κ f − 1 )$ is the characteristic time of the terminal relaxation. Thus, the friction constant of crystal growth can be determined from a fit to the terminal relaxation time as

The kinetic coefficient of crystal growth M in Eq. (1) is inversely proportional to the friction coefficient γ. The proportionality constant depends on the specifics of the simulations cell used. It is found by first taking the average time-derivative of the definition of x_s schematically depicted in Fig. 1, using that $〈 N ̇ s 〉=〈 q ̇ 〉∂〈 N s 〉/∂〈q〉=〈 q ̇ 〉N/ Q s l$ and $γ〈 q ̇ 〉=−α=− μ s l N/ Q s l$⁠,

Since μ_sl is known from Eq. (5), the averaged crystal growth rate $〈 x ̇ s 〉$ can also be deduced.

It is easier to write up the solution of Eqs. (8) and (9) in the frequency domain. Let us consider the complex admittance, i.e., the frequency dependent mobility μ(ω) of Q. From the fluctuations of Q(t), it is possible to compute μ(ω) using the fluctuation-dissipation theorem^6,10 (see the Appendix)

The complex admittance of the model is given by

A derivation of this relation is provided in the Appendix. As shown in Subsection III D, this result can also be obtained using an analogy to an electric circuit that is mathematically equivalent to the stochastic model.

A system of coupled linear equations of motion can be represented as an idealized electrical circuit. The circuit shown in Fig. 2 is equivalent to our stochastic model and we can use the rules of electrical networks to derive the solution. Order parameter values of Q, $q ′ =q− Q ̄$⁠, and f are represented by electric charges while voltage drops correspond to forces acting on these order parameters. The electric circuit consists of two capacitors with the admittances iω/κ and iω/κ_f, two resistors with the admittances 1/γ and 1/γ_f, one inductor with the admittance −i/ωm_f, and one battery with a voltage α. Additionally, the resistors include thermal Johnson-Nyquist noise η_q(t) and η_f(t). The Langevin equations of motion, Eqs. (8) and (9), are retrieved by applying Kirchhoff’s voltage law to the left and the right loops, respectively. The effective Hamiltonian Eq. (7) corresponds to the electric energy of the circuit.

The frequency-dependent admittance (Eq. (14)) over the central capacitor can be obtained as follows: The admittances of the unconnected elements from left to right are μ_q(ω) = 1/γ, μ_Q(ω) = iω/κ, and μ_f(ω) = 1/(iωm_f + γ_f − iκ_f/ω). For the latter, we used Ohm’s law: impedances (i.e., inverse admittances) connected in series are summed. When currents are running in parallel, however, the admittances are summed. With respect to the current $Q ̇$ though the central capacitor, the loops on the left and on the right are connected in parallel with each other and in series with the central capacitor. The frequency dependent admittance, we are interested in is, thus,

By inserting the admittances μ_q(ω), μ_f(ω), and μ_Q(ω), we arrive at Eq. (14).

To validate the use of the interface pinning method to compute crystal growth rates, we investigate a Lennard-Jones system¹¹ of N = 5120 particles (8 × 8 × 20 face centered cubic unit cells; solid-liquid interface in the 100 plane; see Ref. 8 for details). The potential part of the Hamiltonian is $∑ i > j N u ( | r i − r j | )$ with u(r) = 4ε((σ/r)¹² − (σ/r)⁶) − 4ε((1/2.5)¹² − (1/2.5)⁶) for r < 2.5σ and zero, otherwise. Trajectories are generated with the LAMMPS software package¹² using the Verlet integrator with a time step of $0.004σ m / ε$ where m is the particle mass. The Np_zT-ensemble is realized using the Nose-Hoover thermostat^13,14 with a coupling time of $0.4σ m / ε$ and the Parrinello-Rahman barostat¹⁵ with a coupling time of $0.8σ m / ε$⁠.

First, we compute μ(ω) using Eq. (13). The integral is evaluated numerically from a discrete time series Q_i representing Q(t). The rate $Q ̇ i$ is computed from central difference of Q_i: $Q ̇ i = ( Q i + 1 − Q i − 1 ) / ( t i + 1 − t i − 1 )$⁠. Since the autocorrelation function $〈ΔQ ( 0 ) Δ Q ̇ ( t ) 〉$ is zero at t = 0 and t → ∞, we can replace the Fourier-Laplace transform with a regular Fourier transform, and use the efficient Fast Fourier-Transform (FFT) algorithm to evaluate the integral numerically (we note that the FFT algorithm assumes a periodic dataset and we would get an erroneous result if the integrand did not have the property of vanishing values in the limits). For the analysis, we ensure that the discrete Q_i trajectory has a high sampling frequency so that aliasing is avoided (⁠$t i + 1 − t i =0.04σ m / ε$⁠), and are sufficiently long so that slow interface fluctuations are represented (⁠$t total =1 0 6 σ m / ε$⁠). The dots in Fig. 3 show the real and imaginary part of the computed μ(ω) at T = 0.8ε/k_B and the coexistence pressure p_m = 2.185ε/σ³.⁸ The solid lines show μ(ω)/ω of our stochastic model with the four parameters γ, γ_f, κ_f, and m_f determined by a least square fit to the real part. The agreement with both the real and imaginary part, which was not used for the fit of the parameters, is excellent. The fit gives $γ=15.6 ( 3 ) σ m ε$⁠. Here, and throughout the paper, numbers in parenthesis indicate statistical uncertainties on the last digit estimated by dividing runs into statistically independent blocks.¹⁶ Fig. 4 validates that the terminal relaxation time of the autocorrelation function 〈ΔQ(0)ΔQ(t)〉 agrees with this γ, see Eq. (10). Using Eq. (12) with N = 5120, X = Y = 12.82σ, $ρ s =0.973 σ 1 3$⁠, and Q_sl = 56.16,⁸ we arrive at a kinetic coefficient of $M=1.67 ( 3 ) 1 / ε m$⁠. We note that the model does not give a perfect fit in the high-frequency part. This is, however, not important for the estimate of γ, as discussed in Sec. VI. Fig. 5 shows that changing κ over three decades yields consistent results with the model parameters determined from the data shown in Fig. 3.

For comparison, we compute M in a direct simulation of crystal melting (κ = 0). The average interface velocity can be computed as $〈 x ̇ s 〉= Z s /2 t s$⁠, where t_s is the average time to crystallize one box length Z_s = N/ρ_sXY (negative values indicate melting). In a simulation, this time can be computed from the average rate of change of the order-parameter $1/ t s =〈 Q ̇ 〉/ Q s l$⁠. The solid (black) line in Fig. 6 is the average trajectory when melting a crystal at T = 0.8ε/k_B and p = 1.8ε/σ³ (obtained from 50 statistically independent runs). From this result, we find that $〈 x ̇ s 〉=−0.0716 ( 4 ) ε / m$⁠. At this state point μ_sl = 0.0431ε and, thus, we get a kinetic coefficient of $M=1.65 ( 9 ) 1 / ε m$⁠. This is within the statistical uncertainty of the value determined by the interface pinning method.

We used the interface pinning method in combination with ab initio simulations to compute the crystal growth rate for real materials. In Ref. 7, we implemented the interface pinning method in the Vienna Ab initio Simulation Package (VASP)¹⁷ and conducted DFT calculations of the melting temperature of the period three elements sodium (Na) and silicon (Si). These simulations allow us to compute the kinetic coefficients of these elements in the z-direction of the chosen crystal directions.

For computing the trajectories, we employed a Verlet integrator with a timestep of 4 and 3 fs for Na and Si, respectively. The Np_zT-ensemble was realized using a Langevin thermostat with a coupling time of 1 ps and a Parrinello-Rahman barostat¹⁵ with a coupling time of 0.33 ps. The Perdew, Burke, Enzerhof (PBE) exchange correlation functional was used for Na, and the local density approximation (LDA) was used for Si.¹⁸ For the analysis, we use a 120 ps trajectory for Na and 53 ps for Si. In Fig. 7, we compute the kinetic coefficients by fitting our model to the autocorrelation function 〈ΔQ(0)ΔQ(t)〉. Table I lists the computed kinetic coefficients.

We have devised a stochastic model of Q(t) fluctuations in an interface pinning simulation and used it to compute the crystal growth velocity. The critical reader might rightfully ask the question: How dependent is the method on the validity of the stochastic model? To answer this, we emphasize that it is the long-time ω → 0 limit of the model that is of importance. If we use a more complicated model of the fast f(t) fluctuations, the long-time solution is the same, except that κ_f is replaced by an effective spring constant. To motivate the use of a more complicated model, we note that the high-ω part of the μ(ω) spectrum (inset on Fig. 8) indicates a superposition of two peaks. The reason for this is that fluctuations in the bulk of the liquid and the crystal are different. Thus, a better model of Q(t) fluctuations would be to split the fast contribution into contributions of the two phases, f(t) = f_l(t) + f_s(t), and write the effective equations of motion as

The frequency dependent mobility (admittance) of this model is

The long-time limit of this model, however, is the same as in the simpler model we used above, except that κ_f is replaced by an effective spring constant [1/κ_l+1/κ_s]⁻¹.

How does the interface pinning method compare to other methods for computing M? To answer this we note that suggested techniques can be grouped into classes:⁴ (i) free solidification simulations (see Fig. 6),^19–21 (ii) fluctuations analysis of a crystallinity order-parameter,^22–24 (iii) capillary wave analysis,^4,25–27 and (iv) forced velocity simulations.^28,29 This paper’s method is of class (ii). The interface pinning method can be viewed as a generalization of the approach suggested by Briels and Tepper.^22,23 In the Briels-Tepper method two-phase configurations are stabilized by simulating the constant NVT ensemble. This can be viewed as a special case of the interface pinning method, where the order-parameter is the volume, and the bias potential has an infinitely large spring constant. Fluctuations in the number of crystalline particles can be monitored by pressure fluctuations or an order-parameter of crystallinity (e.g., Eq. (3)). An advantage that the interface pinning method inherits from the Briels-Tepper method is that it involves well-defined equilibrium computations that can be done ad infinitum. The following challenges of the Briels-Tepper method are solved with the interface pinning method: (a) In the Briels-Tepper method, the two-phase configurations are stabilized by keeping the volume constant. In effect, computations can only be done near coexistence. In the interface pinning method, the stabilization of two-phase configurations is done by connecting the system with a harmonic field that couples to an order-parameter Q that distinguishes between the crystal and the liquid, and simulations can be performed far into the supercooled or superheated regimes. Also, the size of the fluctuations in the number of crystalline particles is determined by the compressibility of the solid and liquid. In effect, they will be large for large systems, leading to long correlation times. With the interface pinning method, the size of the N_s fluctuations can be controlled by the value of κ and by choosing a good order-parameter making κ_f large. (b) μ_sl is directly computed with the interface pinning method. This is an essential property that otherwise had to be computed in separate calculations.

The interface pinning method also bears some similarities with the capillary wave analysis method^4,25–27 as both approaches relate the interface mobility to the time correlations of spontaneous interface fluctuations based on the fluctuation-dissipation theorem. More specifically, in the capillary wave analysis method, one tracks the position and shape of the crystal-liquid interface and determines its k-space fluctuation spectrum. By combining the decay time of the correlations with the interface stiffness obtained from static correlations, one can then compute the interface mobility. The interface pinning method and the capillary wave analysis method, however, also differ in important ways. While in the interface pinning method, one considers only fluctuations of the global parameter Q(R) that quantifies the fraction of the system in the crystalline state, the capillary wave analysis method requires an exact localisation of the crystal-liquid interface. To do that one needs to be able to differentiate locally between solid-like and liquid-like particles, a distinction that is not necessary in the interface pinning method. Also, the capillary wave analysis method can be applied only to coexisting phases, thus, requiring knowledge of the melting temperature, while the interface pinning method provides the interface mobility away from coexistence too.

In this paper, we consider the hydrostatic Np_zT-ensemble^7,8 near coexistence. Growth or melting of a solid results in latent heat and latent volume that must be removed to stay in the Np_zT-ensemble. To avoid that the growth rate is trivially determined by the characteristic time of the thermostat and barostat, we investigated crystal growth dynamics with a strong coupling of the thermo- and barostat, i.e., short coupling times. We note that in a weak coupling or experimental situation latent heat and volume disperse away from the interface. In effect, the temperature and pressure at the interface are different from that at the boundaries of the system,^24,29 and in general the temperature and strain tensor fields should be taken explicitly into consideration. Addressing this issue goes beyond the scope of this paper.

Auxiliary harmonic potentials as the ones employed here are routinely used to compute free energies along reaction coordinates for crystal nucleation or other rare events. Specifically, the umbrella sampling method^30,31 uses a series of auxiliary potentials referred to as umbrellas. By reweighting³² probability distributions generated in the presence of the umbrellas, it is possible to compute free energies along a given reaction coordinate. As we have demonstrated here, one can also extract crystal growth rates without further computations. By extending the bias potentials to the regime where only a small crystallite is present in the system, one may be able to compute crystal nucleation rates from IP-simulations. We leave such studies for future investigations.

Fig. 8 shows the kinetic coefficient M along the coexistence line. In the region T < 1.2, M is nearly constant. At higher temperatures, M decreases indicating that the crystal growth (or melting) rate is lower. The inset shows that M in reduced units, using the thermal energy for the energy scale, is roughly invariant at T > 1.2. In other words, M scales as $1/ T$ along the high temperature (and pressure) part of the melting line. This scale invariance is expected if a mapping to the hard-sphere system is valid²⁶ and is predicted by the “isomorph theory”³³ of simple liquids.³⁴ This theory states that there is a class of systems,^35–37 including the LJ system,³⁵ that have “isomorph”-lines in the dense and/or high temperature part of the phase diagram (i.e., not near the critical point). Along these lines structure, dynamics and some thermodynamic properties are invariant in reduced units. One prediction is that the melting line is an isomorph,³³ and indeed this was shown for the LJ system in Refs. 8 and 38. Thus, we expect M to be invariant in reduced units along the melting line. This explains why M scales as $1/ T$⁠. The scaling law does not apply near the triple point (note that the lowest temperature data-point on Fig. 8 is at negative pressure). Consistent with this, the melting line itself deviates from an isomorph near the triple point.⁸ This deviation is probably due to long-ranged attractive interactions.⁸ We will leave further investigation of isomorph-scale invariance of the crystallization to future studies.

This work was financially supported by the Austrian Science Fund FWF within the SFB ViCoM (F41). U.R.P. was funded by the Villum Foundation’s Grant No. VKR-023455. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). The authors are grateful for valuable comments and suggestions from Georg Kresse, Jeppe C. Dyre, Lorenzo Costigliola, Tage Christensen, and Peter Harrowell.

In this appendix, we provide a detailed analysis of the stochastic model suggested in the main part of the paper.

The static averages and fluctuations of q and f (and hence also of Q) can be computed as averages over the equilibrium distribution

where β = 1/k_BT and the partition function

normalizes the distribution. The moments of this multivariate Gaussian distribution can be determined analytically yielding the averages

and variances

Here, Δq = q − 〈q〉, Δf = f − 〈f〉, and $Δ f ̇ = f ̇ −〈 f ̇ 〉$ are the deviations of q, f, and $f ̇$ from their respective averages. While $f ̇$ is uncorrelated to the other variables, q and f are correlated with covariance

From the above expressions, the average and variance of Q = q + f follow

Thus, the static fluctuations of Q depend only on the temperature and the force constant κ of the pinning potential.

To quantify the average dynamics of the model, we now introduce the time correlation functions

which correlate the state of the system at time 0 to its state at a time t later. It follows from the microscopic reversibility of the equations of motion that ϕ_qf(t) = ϕ_fq(t), such that we do not need to consider ϕ_qf(t) and ϕ_fq(t) separately. Thus, the time autocorrelation functions of the observable order parameter 〈ΔQ(0)ΔQ(t)〉 are

One can show, that the following equations hold for the time correlation function 〈α(0) β(t)〉, where α and β are two arbitrary functions of q and f,

Upon time reversal, the time correlation function transform as

where ε_α and ε_β take values +1 or −1, depending on whether α and β, respectively, are even or odd under reversal of the momenta.

By averaging over the equations of motion, the following differential equations for the time correlations can be derived:

where we have omitted the argument of the time correlation functions for simplicity. Introducing the auxiliary function $ψ ff ( t ) = ϕ ̇ ff ( t )$⁠, we write these differential equations in matrix notation as

In a more compact form, this system of homogeneous linear first order differential equations with constant coefficients is expressed as

where $x= ϕ qq , ϕ qf , ϕ ff , ψ ff$ and A is the constant 4 × 4 matrix on the right hand side of Eq. (A22). The formal solution of this equation is given by

with initial conditions

As can be seen in Eq. (A22), the time evolution equations for ϕ_qf(t), ϕ_ff(t), and ψ_ff(t) are independent of ϕ_qq(t). Hence, one can obtain the time correlations functions ϕ_qf(t), ϕ_ff(t), and ψ_ff(t) by solving

As a consequence, the three time correlation functions ϕ_qf(t), ϕ_ff(t), and ψ_ff(t) can be written as superpositions of three exponentials,

where x_i(t) is component i of the time correlation function vector x(t). (Note that we start our numbering at 0 such that the 0-th component of x(t) is the time correlation function ϕ_qq(t)). The time constants λ_i are the eigenvalues of the 3 × 3 matrix on the right hand side of Eq. (A26) and the specific values of the constants a_ij can be determined by diagonalising this matrix and applying the initial conditions ϕ_qf(0), ϕ_ff(0), and ψ_ff(0). The eigenvalues λ_i, which are also eigenvalues of the matrix A and can be complex for certain parameter combinations, are the roots of the cubic equation

Real roots correspond to exponential decay, while complex roots lead to oscillatory behaviour. In general, the time correlation functions ϕ_qq(t), ϕ_qf(t), ϕ_ff(t), and ψ_ff(t) are superpositions of exponential and oscillatory terms. Whether the roots of this equation are real or complex can be determined by computing the discriminant D of the above polynomial in λ. This polynomial discriminant, as well as the roots of the equation, can be expressed explicitly in terms of the constants κ, κ_f, γ, γ_f, and m_f, but the expressions are omitted here because they are complicated and do not provide much insight. If D > 0, one root is real and two are complex conjugate, if D < 0 all roots are real and different, and if D = 0, all roots are real and at least two of them are equal.

Once the time correlation function ϕ_qf(t) is known, ϕ_qq(t) can be determined by solving the remaining differential equation

Since ϕ_qf(t) is already given, this equation is an inhomogeneous linear differential equation with constant coefficients, which can be solved, for instance, by variation of constants yielding

Here, λ₀ = κ/γ and the coefficient a₀ depends on the initial conditions. Like λ₁, λ₂, and λ₃, also λ₀ is an eigenvalue of the matrix A. Thus, in general, the time correlation function ϕ_qq(t) is a sum of four exponentials, two of which can lead to oscillatory behaviour. However, explicit solutions of the differential equations for specific parameter sets indicate that the coefficient a₀ is close to zero, such that also ϕ_qq(t) is, effectively, a sum of three exponentials.

While it is possible to compute the roots of Eq. (A28) analytically, the resulting expressions are exceedingly complicated. For the root closest to zero, which dominates the behaviour of the correlation functions for long times, a simple but accurate approximation can be derived. Neglecting the cubic and quadratic term, the root $λ ̃$ with the smallest magnitude is given by

As a further simplification we will use a “separation of time-scales” approximation $γ ( κ − 1 + κ f − 1 ) ≫ γ f κ f − 1$ to eliminate the last term. Thus, the terminal exponential relaxation time is

justifying Eq. (10) in the main part of the paper. We note that the term τ_f = γ_f/κ_f that we neglect is the characteristic time of the uncoupled (κ = 0) fast f vibrations. For the typical system studied in this paper, this term is less than 1% of τ.

According to the fluctuation-dissipation theorem,^6,10 the response of a system to a weak perturbation can be related to the fluctuation properties of the system in equilibrium. To apply this theorem to the pinned interface, imagine a system with the Hamiltonian $H ( z )$ perturbed by a time dependent external field K(t) that couples linearly to the variable A(z). Here, z denotes a set of variables describing the state of the system, for instance, the positions and momenta of all particles in our simulation or the variables Q, q, f, and $f ̇$ of our stochastic model. The time-dependent Hamiltonian of the perturbed system is then given by

The reaction of the system to the external perturbation can be monitored through the change in the variable B(z), which, like A, is also a function of z, with respect to its equilibrium average,

where the angular brackets denote an equilibrium average and B(t) is a shorthand notation for B[z(t)]. Assuming that the external force K(t) has been acting on the system from the infinite past, i.e., from t = − ∞, the average deviation 〈ΔB(t)〉_ne at time t can be written in its most general form as

where the response function R_BA(t) denotes the response of the system to an impulsive force, i.e., K(t) ∝ δ(t) applied to the system at time t = 0. The angular brackets 〈⋯〉_ne imply an average over many realisations of the process. According to the fluctuation dissipation theorem, the response function is related to the fluctuation properties of the system by

where ΔA(t) = A(t) − 〈A〉 denotes the deviation of A(t) = A[z(t)] from its equilibrium average. Hence, in the linear regime, the response of a system to an external perturbation is related to the correlation of spontaneous fluctuations at different times in the equilibrium system.

The significance of the fluctuation dissipation theorem becomes particularly clear in the frequency domain. In this case, the linear relation between perturbation and response is expressed as

where

and

are the Fourier transforms of the response and the force, respectively. In the frequency representation, the fluctuation-dissipation theorem then links the complex admittance χ_BA(ω) to the Fourier-Laplace transform of the time correlation function $〈Δ A ̇ ( 0 ) ΔB ( t ) 〉$⁠,

In the following, we will mainly use this frequency-dependent formulation of the fluctuation-dissipation theorem.

We next apply the fluctuation-theorem to our stochastic model with the Hamiltonian of Eq. (7) and evolving according to the Langevin equations (8) and (9). For this model, we will compute the frequency dependent complex admittance that describes the reaction of the system to an external perturbation acting on the variable Q = q + f. The response of the system to the perturbation is monitored using $Q ̇$⁠, the velocity associated with the variable Q. So what we would like to determine is the complex admittance μ(ω) that relates the average velocity $Q ̇$ to the strength of the external perturbation

By comparing the computed admittance μ(ω) with results of computer simulations carried out for the atomistic system, we will then determine the values of the model parameters. In particular, we will determine the parameter γ which describes the mobility of the interface driven by the difference in chemical potential between the phases.

For the perturbation QK(t) = (q + f) K(t), the Langevin equations of motion (8) and (9) turn into

Since these equations are linear, we can compute the response of the system analytically for an arbitrarily strong external force K(t). To do that, we carry out a Fourier transformation on the Langevin equations above and average over many realizations of the stochastic process. For ω ≠ 0, we obtain

where $Δ q ̃ ( ω )$ and $Δ f ̃ ( ω )$ are the Fourier transforms of Δq(t) = q(t) − 〈q〉 and Δf(t) = f(t) − 〈f〉, respectively. Note that in the above equations we have omitted the argument ω in the averages to simplify the notation. Since, we are interested in the response of the system in terms of the generalized velocity $Q ̇$⁠, we rewrite these equations for the Fourier transforms of the time derivatives of Δq and Δf,

Here, we have exploited that in frequency space taking a time derivative simply amounts to multiplication with iω, such that $〈Δ q ̇ ̃ 〉 ne =iω〈Δ q ̃ 〉 ne$ and $〈Δ f ̇ ̃ 〉 ne =iω〈Δ f ̃ 〉 ne$⁠.

To simplify the notation in the following, we introduce the complex admittances μ_q(ω) and μ_f(ω) for $Δ q ̇$ and $Δ f ̇$ separately in the absence of the coupling term $−κ ( q + f − Q ̄ ) 2$ in the Hamiltonian of Eq. (7), i.e., without pinning potential. In this case, which corresponds to κ = 0, the equations of motion yield

such that

Using μ_q(ω) and μ_f(ω), Eqs. (A46) and (A47) for the system with bias can be written as

By equating the right hand sides of the above equations, one obtains

implying that

The response of the system in terms of $Q ̇ = q ̇ + f ̇$ is then given by

such that the complex admittance μ(ω) is given by

equivalent with Eq. (14) in the main part of the paper. The fluctuation-dissipation theorem links the complex admittance to the Fourier-Laplace transform of the equilibrium time correlation function of $Δ Q ̇$⁠,

Thus, the complex admittance μ(ω) can be obtained by Fourier-Laplace transformation of the time correlation function $〈Δ Q ̇ ( 0 ) Δ Q ̇ ( t ) 〉$ and vice versa with an inverse Fourier-Laplace transformation.