Computing the crystal growth rate by the interface pinning method

An essential parameter for crystal growth is the kinetic coefficient given by the proportionality between super-cooling and average growth velocity. Here we show that this coefficient can be computed in a single equilibrium simulation using the interface pinning method where two-phase configurations are stabilized by adding an spring-like bias field coupling to an order-parameter that discriminates between the two phases. Crystal growth is a Smoluchowski process and the crystal growth rate can therefore be computed from the terminal exponential relaxation of the order parameter. The approach is investigated in detail for the Lennard-Jones model. We find that the kinetic coefficient scales as the inverse square-root of temperature along the high temperature part of the melting line. The practical usability of the method is demonstrated by computing the kinetic coefficient of the elements Na, Mg, Al and Si from first principles. It is briefly discussed how a generalized version of the method is an alternative to forward flux sampling methods for computing rates along trajectories of rare events.


I. INTRODUCTION
Crystal growth is of paramount importance in many branches of condensed matter physics [1,2]. A important parameter in phase field equations [3] describing crystal growth is the kinetic coefficient defined as the proportionality constant between super cooling (or super heating) and the average interface growth (or melting) velocity ẋ 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 super cooling (or super heating): In the spirit of the fluctuation-dissipation theorem [6] we suggests 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 assumes 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 * ulf@urp.dk Sketch of an elongated periodic simulation box containing a two-phase configuration. The projected interface position, ignoring capillary waves, is xs = Ns/2XY ρs.
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 model and combine the method with ab initio DFT to compute growth rates of the elements Na, Mg Al and Si. The paper is finalized with a discussion of the method (Sec. VI). The Appendix gives a detailed analysis of our stochastic model.

II. INTERFACE PINNING
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 1 , r 2 , . . . , r N ) let the system have the Hamiltonian H 0 (R). Hydrostatic pressure is ensured by allowing the box length in the zdirection to change and constant temperature is ensured by connecting the momenta of the particles to a heat bath. We refer to this as the N p z T -ensemble. When the initial R is a two-phase configuration, the thermodynamically stable phase will grow on the expense of the arXiv:1406.2659v2 [cond-mat.stat-mech] 11 Jul 2014 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 pointQ 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)| [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 B T /κ. 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 the following section 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 .

III. STOCHASTIC MODEL
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Ṙ. 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 [9]. 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.

A. Effective Hamiltonian
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 κ andQ appearing in the pinning potential denote the force constant and the position of the bias, respectively. The third term, (κ f /2)f 2 , reflects the harmonic potential energy experienced by phonon fluctuations. Finally, the last term is the kinetic energy associated to the time derivativeḟ 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.

B. Langevin equations of motion
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 reflects the degrees of freedom neglected in the simplified model [10]. The variance of the random forces and the friction constants are related by the fluctuationdissipation theorem [6], η q (0)η q (t) = 2γk B T δ(t) and 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.

C. Solution of the stochastic model
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 γ( where A = κ f /(κ + κ f ) is the strength and τ = γ(κ −1 + κ −1 f ) 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 depictured in Fig. 1, using that Ṅ s = q ∂ N s /∂ q = q N/Q sl and γ q = −α = −µ sl N/Q sl : Since µ sl is known from Eq. (5), the averaged crystal growth rate ẋ s can also be deduced. It is easier to write up the complete solution 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 (see the Appendix) In the following subsection we derive this using an analogy to an electric circuit that is mathematically identical to the stochastic model.

D. Electric network representation
A system of coupled linear equations of motion can be represented as an idealized electrical circuit. The circuit shown 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 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 loop, 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 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 that are summed. With respect to the currentQ 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 we arrive at Eq. (14).

IV. THE LENNARD-JONES SYSTEM
To validate the use of the interface pinning method to compute crystal growth rates, we investigate a Lennard-Jones system [11] 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 for r < 2.5σ and zero otherwise. Trajectories are generated with the LAMMPS software package [12] using the Verlet integrator with a time step of 0.004σ m/ε where m is the particle mass. The N p z T -ensemble is realized using the Nose-Hoover thermostat [13,14] with a coupling time of 0.4σ m/ε and the Parrinello-Rahman barostat [15] 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 rateQ i is computed from cen- 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. 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ε/σ 3 [8]. The solid lines shows µ(ω)/ω 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.8σ √ mε. 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.6 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 ẋ s = Z s /2t s where t s is the average time to crystallize one box length Z s = N/ρ s XY (nega-FIG. 6. Time evolution of Q along 50 statistically independent LJ melting runs (gray) at T = 0.8ε/kB and p = 1.8ε/kB where the initial configurations are taken from equilibrated interface pinning simulations (κ = 4ε;Q = 26). At this state-point the liquid is the thermodynamically stable phase (µ sl = 0.043ε) and the crystal melts. From the average rate of (red dashed) we compute the average growth velocity to ẋs = −0.0716 ε/m corresponding to M = 1.6 1/εm (using that Zs = 32.1σ). The inset shows ẋs along the tive 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 sl . The 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ε/σ 3 (obtained from 50 statistically independent runs). From this result we find that ẋ s = −0.0716 ε/m. At this state point µ sl = 0.0431ε and, thus, we get a kinetic coefficient of M = 1.6 1/εm. This is the same as the value determined by the interface pinning method.

V. FIRST PRINCIPLE COMPUTATIONS
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) [16] and conducted density functional theory (DFT) calculations of the melting temperature of the period three elements sodium (Na), magnesium (Mg), aluminum (Al) and silicon (Si). Without any additional computations, data obtained from these simulations allow us to compute the kinetic coefficients of these elements in the z-direction of the chosen crystal  Table I and Ref. [7] for details). A fit to the terminal exponential relaxation (solid) yields a kinetic coefficient of M = 120Å/(ps eV).

directions.
For computing the trajectories we employed a Verlet integrator with a timestep of 4 and 3 fs for Na and the other elements, respectively. The N p z T -ensemble was realized using a Langevin thermostat with a coupling time of 1 ps and a Parinello-Rahman barostat [15] with a coupling time of 0.33 ps. PBE and LDA functionals [17] were used for the DFT calculations.
We compute the kinetic coefficients from the terminal relaxation times of the autocorrelation function ∆Q(0)∆Q(t) of the order parameter used for the interface pinning calculation. Fig. 7 shows the autocorrelation function obtained in a first principles calculation of Al. From a fit to the terminal relaxation we compute M = 120Å/(ps eV) for the kinetic coefficient in the (100) crystal direction. Table I lists the computed kinetic coefficients for the remaining elements along with the direction of crystal growth and the settings used for the interface pinning calculation.

A. Dependence on stochastic model
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 longtime 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 ] −1 .

B. Other techniques
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 three classes [4]: (i) free solidification simulations (see Fig. 6) [19][20][21], (ii) fluctuations analysis [22][23][24] and (iii) forced velocity simulations [25,26]. 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 twophase configurations are stabilized by simulating the constant N V T ensemble. This can be viewed as a special case of the interface pinning method, where the orderparameter 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 twophase 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 are 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 super cooled or super heated 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.

C. Latent heat and volume
In this paper we consider the hydrostatic N p z Tensemble [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 N p z T -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 an weak coupling or experimental situation latent heat and volume disperse away from the interface. In effect, the temperature and pressure at the interface is different from that at the boundaries of the system [24,26], and in general the temperature and strain tensor fields should be taken explicitly into consideration. Addressing this, goes beyond the scope of this paper.

D. An alternative to forward flux sampling
We have demonstrated that it is possible to compute the rate of a non-equilibrium process from an equilibrium simulation where an auxiliary harmonic potential is added the Hamiltonian. Auxiliary harmonic potentials are routinely used to compute free energies along reaction coordinates of rare events. Specifically, the umbrella sampling methods [27,28] uses a series of auxiliary potentials refereed to as umbrellas. By reweighting [29] probability distributions generated by the umbrellas it is possible to compute free energies along order parameters. As we have demonstrated here, it is possible to extract also the rates without further computations. We investigated the growth velocity of a flat interface, but presumably the same approach can be used for studying other reaction paths. As an example, by setting up a system with a crystallite [30][31][32], it should be possible to compute the growth velocity of a curved solid-liquid interface. We leave such studies for future investigations.
A generalized version of the method may be used for studying dynamics along reaction paths of other rare events, including those of biological systems [33][34][35]. Thus, our method provides an alternative to the forward flux sampling methods [36]. The original version of these methods [37] use a series of interfaces along the path of interest. The rate along a given path is then computed by generating short trajectories between the interfaces. Forward flux sampling is limited to stochastic process; this restriction is not needed with the harmonic pinning approach presented in this paper.
E. Scaling of the kinetic coefficient along the Lennard-Jones melting line 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 show 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 predicted by the "isomorph theory" [38] of simple liquids [39]. This theory states that there is a class of system [40][41][42], including the LJ system [40], 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 [38], and indeed this was shown for the LJ system in Refs. [8,43] the triple point (note that the lowest temperature datapoint on Fig. 8 is at negative pressure). Consistent with this, the melting line itself deviates from an isomorph near the triple point [8]. This deviation is probably due to long-ranged attractive interactions [8]. We will leave further investigation of isomorph-scale invariance of the crystallization to future studies.

ACKNOWLEDGMENTS
This work was financially supported by the Austrian Science Fund FWF within the SFB ViCoM (F41). URP was in part funded by the Villum Foundation's grand VKR-023455. Simulations were carried out on 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.

Appendix A: Analysis of stochastic model
In this Appendix we provide a detailed analysis of the stochastic model suggested in the main part of the paper.

Static averages and variances
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 B T 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 ∆ḟ =ḟ − ḟ are the deviations of q, f , andḟ from their respective averages. Whileḟ 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.

Time correlation functions
To quantify the average dynamics of the model we now introduce the time correlation functions φ qq (t) = ∆q(0)∆q(t) , 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) = φ f q (t), such that we don't need to consider φ qf (t) and φ f q (t) separately. Thus the time auto correlation functions of the observable order parameter ∆Q(0)∆Q(t) is 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 (A22) 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 Equ. (A22). The formal solution of this equation is given by with initial conditions (A25) As can be seen in Equ. (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 Equ. (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 (A28) 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, γφ qq (t) = −κ(φ qq (t) + φ qf (t)). (A29) Since φ qf (t) is already given, this equations is an inhomogeneous linear differential equations with constant coefficients, which can be solved, for instance, by variation of constants yielding Here, λ 0 = κ/γ and the coefficient a 0 depends on the initial conditions. Like λ 1 , λ 2 , and λ 3 , also λ 0 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 0 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 Equ. (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 γ( 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 τ .

The fluctuation-dissipation theorem and the solution in the frequency domain
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 simulation or the variables Q, q, f andḟ 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 the process. According to the fluctuation dissipation theorem, the response function is related to the fluctuation properties of the system by 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 an response is expressed as where 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 ∆Ȧ(0)∆B(t) , In the following, we will mainly use this frequencydependent 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 usingQ, the velocity associated with the variable Q. So what we would like to determine is the complex admittance µ(ω), that relates the average velocityQ to the strength of the external perturbation, ∆Q(ω) ne = µ(ω)K(ω). (A41) 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 γq + α + κ(f + q −Q) = K(t) + η q (t), (A42) 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 (κ + i ωγ) ∆q ne + κ ∆f ne =K(ω), (A44) 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 velocityQ, 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 .
The response of the system in terms ofQ =q +ḟ is then given by such that the complex admittance µ(ω) is given by .
(A57) 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.