text:
# Langevin dynamics of the Coulomb frustrated ferromagnet: a mode-coupling analysis.

## I Introduction

Models with a competition between a short-range ordering interaction and a long-range frustrating interaction have been recently introduced to explain the slowing down of relaxation in supercooled liquids and the resulting glass transitionKKZNT95 . The underlying picture is as follows: any given liquid possesses a locally preferred structure which is different than that of the actual crystalline phase, and this local arrangement of the molecules in the liquid cannot propagate at long distances to tile the whole space and form an “ideal crystal” because of ubiquitous frustration. It has been arguedKKZNT95  that, for weak enough frustration, this phenomenon can be described via effective interactions acting on very different length scales: a short-range term describing the tendency to extend the locally preferred structure and a long-range, Coulomb-like term describing the frustration-induced free-energy cost associated with this spatial extension. Such Coulomb frustrated systems have been shown, both by scaling argumentsKKZNT95 ; VTK00  and by Monte Carlo simulationGrTV01 , to display the generic features observed in fragile glassforming liquids, most notably the super-Arrhenius temperature dependence of the relaxation time and the two-step, non-exponential decay of the correlation function.

These Coulomb frustrated models have also been used in quite different contexts to describe the formation of modulated spatial patterns on mesoscopic length scales, such as lamellar and cubic phases in diblock copolymer meltsOK86 ; L80 ; FH87 , microemulsions in water-oil-surfactant mixturesS83 ; WCS92 , or stripe phases in high temperature superconductorsEK93 . In all these cases, slow relaxation is usually observed, and it has been recently arguedScWo00  that high-temperature superconductors could indeed form a “stripe glass” in which glassiness is self generated, i.e., does not result from the presence of quenched disorder. This latter result has been obtained through an investigation of the properties of the free-energy landscape of the Coulomb frustrated $`\varphi ^4`$ scalar field theory: by using a thermodynamic approach combining the replica method proposed for the study of structural glassesM95 ; MP99  and a particular approximation, the self-consistent screening approximation (SCSA)AJBa74 , for calculating the pair correlation functions, Schmalian and WolynesScWo00  have derived that the free-energy landscape of the Coulomb frustrated model becomes non-trivial below some temperature $`T_A`$ at which an exponentially large number of metastable states appears; the associated configurational entropy decreases with further decrease of the temperature and vanishes at a lower temperature $`T_K`$ScWo00 ; SW201 ; endnote1 .

Motivated by these results giving evidence for fragile glassforming behavior in Coulomb frustrated models, we have studied the Langevin dynamics of the Coulomb frustrated $`\varphi ^4`$ scalar field theory within the mode-coupling and related approximations. Mode-coupling approaches have been widely used to study glassforming liquidsgotze91 ; Cummins99 , and the dynamical ergodicity-breaking singularity predicted to occur in the weakly supercooled liquid regions, albeit “avoided” in real systems, is taken by many as a canonical feature of fragile glassforming systems. It is therefore tempting to investigate whether Coulomb frustrated models also display this feature.

The paper is organized as follows. We first present the model and summarize the equilibrium phase behavior and the results previously obtained by Schmalian and Wolynes. We also introduce the Langevin equation describing the relaxational dynamics of the system. In section III, we derive the evolution equations followed by the equilibrium time-dependent correlation function obtained within two resummation schemes of perturbative expansions: the mode-coupling approximation and the dynamical SCSA. Section IV is devoted to the search for an ergodicity-breaking transition. We find that such a phenomenon is indeed observed with the two approximations considered. We also show that the dynamical singularity predicted by the dynamical SCSA coincides with the temperature $`T_A`$ at which the replica analysis of RefsScWo00 ; SW201  predicts the occurrence of an exponentially large number of metastable states. In section V, we present the full numerical solution of the mode-coupling equations, thereby obtaining the time-evolution of the equilibrium correlation function; this latter is similar to that obtained in the idealized mode-coupling theory of supercooled liquidsgotze91  and in mean-field generalized spin-glass modelsKT87 ; BCKM96 ; CHS93 . In the last two section, we address the question of sensitivity of the results to the level and the details of the approximation scheme and we give some concluding remarks.

## II Model

We consider the field-theoretical version of the $`3`$-dimensional Coulomb frustrated Ising ferromagnet defined by the Hamiltonian

$`H[\varphi ]=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3𝐱\left\{(\varphi (𝐱))^2+r_0\varphi ^2(𝐱)+\frac{u}{2}\varphi ^4(𝐱)\right\}}+{\displaystyle \frac{Q}{8\pi }}{\displaystyle d^3𝐱d^3𝐱^{}\frac{\varphi (𝐱)\varphi (𝐱^{})}{|𝐱𝐱^{}|}}`$ (1)
$`=`$ $`{\displaystyle \frac{V}{2}}{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\left(r_0+k^2+\frac{Q}{k^2}\right)\varphi _𝐤\varphi _𝐤}+{\displaystyle \frac{uV}{4}}{\displaystyle \frac{d^3𝐤_\mathrm{𝟏}}{(2\pi )^3}\frac{d^3𝐤_\mathrm{𝟐}}{(2\pi )^3}\frac{d^3𝐤_\mathrm{𝟑}}{(2\pi )^3}\varphi _{𝐤_\mathrm{𝟏}}\varphi _{𝐤_\mathrm{𝟐}}\varphi _{𝐤_\mathrm{𝟑}}\varphi _{𝐤_\mathrm{𝟏}𝐤_\mathrm{𝟐}𝐤_\mathrm{𝟑}}},`$ (2)
where $`\varphi (𝐱)`$ is a real scalar field ($`\varphi _𝐤`$, the associated $`𝐤`$Fourier component), $`V`$ is the volume, $`u`$ is a strictly positive coupling constant, $`Q`$ is the frustration parameter, and all momentum integrations are performed up to a cut-off $`\mathrm{\Lambda }`$, i.e., $`|𝐤|\mathrm{\Lambda }`$; $`r_0`$ is a temperature-dependent mass which is proportional to the deviation, $`TT_{c,MF}^0`$, from the mean-field transition temperature of the unfrustrated $`(Q=0)`$ model.

The equilibrium partition function is

$$Z=𝒟\varphi e^{H(\varphi )/T}.$$
(3)
In what follows, we take $`\mathrm{\Lambda }=1`$, and $`T`$ is set equal to $`1`$ in Eq. (3) so that the whole temperature dependence is contained in $`r_0`$. We are interested in the weak-frustration region for which $`Q<<1`$.

In the absence of frustration $`(Q=0)`$, the model defined by Eqs. (1-3) reduces to the usual $`\varphi ^4`$ theory. It undergoes a second-order transition at a finite temperature $`T_c^0`$ to a broken-symmetry phase characterized by a non-zero value of $`<\varphi _{𝐤=\mathrm{𝟎}}>`$, where $`<>`$ denotes an equilibrium average. For $`Q>0`$, an ordered phase with $`<\varphi _{𝐤=\mathrm{𝟎}}>0`$ is forbidden, but the system can still undergo a phase transition at a temperature $`T_{DO}(Q)`$ to a phase with long-range modulated order. This transition has been studied by Monte Carlo simulation for the case of the Coulomb frustrated Ising ferromagnet on a cubic latticeGTV01b  and via the self-consistent Hartree approximation for a Hamiltonian similar to Eq. (1) describing microphase separation in diblock copolymer meltsFH87 : it has then been shown that, whereas the mean-field theory predicts a second-order transition, the fluctuations change the order of the transition and induce a first-order transition. Such a fluctuation-induced first-order transition was first discussed by Brazovskii for a related modelB75 . Within the self-consistent Hartree approximation, the equilibrium (connected) correlation function

$$VC(𝐤)=<\varphi _𝐤\varphi _𝐤><\varphi _𝐤><\varphi _𝐤>$$
(4)
is obtained via a self-consistent equation; for instance, in the paramagnetic phase where $`<\varphi _𝐤>=0`$, this equation is

$`C^1(𝐤)=r_0+k^2+{\displaystyle \frac{Q}{k^2}}+3u{\displaystyle \frac{d^3𝐪}{(2\pi )^3}C(𝐪)}.`$ (5)
The renormalized mass, $`r=r_0+3u\frac{d^3𝐪}{(2\pi )^3}C(𝐪)`$, is then given by

$$r=r_0+3u\frac{d^3𝐪}{(2\pi )^3}\frac{1}{r+q^2+\frac{Q}{q^2}}.$$
(6)
Since $`k^2+\frac{Q}{k^2}`$ is minimum for non-zero wave vectors with modulus $`k_m=Q^{1/4}`$, a value characterizing the incipient modulated order, one easily checks that $`r`$ only goes to zero when $`r_0\mathrm{}`$. This means that the paramagnetic phase is (meta)stable at all finite “temperatures”, its spinodal being depressed to $`r_0\mathrm{}`$. The Hartree approximation allows one to calculate the free energy of the paramagnetic phase and that of the phase with modulated order. One then obtains the temperature $`r_0^{DO}(Q)`$ of the first-order transition at the point at which the two free energies are equal. The details are given in Appendix A and the resulting phase diagram is shown in Fig. 1.

In their recent work, Schmalian and WolynesScWo00  have applied the thermodynamic approach of non-random glassforming systems developed by Mézard and ParisiMP99  to the Coulomb frustrated model. The basic idea, originally motivated by the behavior of a class of mean-field generalized spin glasses such as $`p`$spin and Potts glassesKT87 ; BCKM96 , is that glassiness arises because of the occurrence of an exponentially large number of metastable states. This occurrence, and the associated emergence of a non-zero complexity or configurational entropy, can be more conveniently studied within the replica formalismM95 . Approximations are of course necessary to solve the corresponding many-body problem and obtain all relevant correlation functions. Schmalian and Wolynes have shown that within the SCSAAJBa74 , an approximation that goes beyond the Hartree result in that it explicitly includes more diagrams of the perturbative expansionendnote2 , there is a temperature $`T_A`$ at which an exponentially large number of metastable states emerges, as signaled by a non-zero configurational entropy. The configurational entropy decreases with further decay of the temperature and it vanishes at a temperature $`T_K`$ at which the system undergoes a random first-order transition to an “ideal glass”ScWo00 .

In the present work, we focus on the dynamics of the Coulomb frustrated model defined by Eq. (1). The starting point is the Langevin equation,

$$\frac{\varphi _𝐤(t)}{t}=\frac{\delta (H[\varphi _𝐤(t)]/V)}{\delta \varphi _𝐤(t)}+\eta _𝐤(t),$$
(7)
that describes the purely relaxational dynamics of the system; $`\eta _𝐤(t)`$ is a gaussian thermal noise with $`<\eta _𝐤(t)>=0`$ and $`<\eta _𝐤(t)\eta _𝐤^{}(t^{})>=2T\delta (𝐤+𝐤^{})\delta (tt^{})`$. Eq. (7) can be explicitly written as

$`{\displaystyle \frac{\varphi _𝐤(t)}{t}}`$ $`=(r_0+k^2+{\displaystyle \frac{Q}{k^2}})\varphi _𝐤`$
$``$ $`u{\displaystyle \frac{d^3𝐤_1}{(2\pi )^3}\frac{d^3𝐤_2}{(2\pi )^3}\varphi _{𝐤_1}\varphi _{𝐤_2}\varphi _{𝐤𝐤_1𝐤_2}}+\eta _𝐤(t).`$ (8)
Solving this set of coupled non-linear dynamical equations is a daunting task, and virtually all available approximations amount to performing some self-consistent resummation of perturbative expansions, e.g,. expansions in powers of the coupling constant $`u`$ or of the inverse of the number of components of the field, $`1/n`$, for an $`O(n)`$ model. In what follows, we shall consider two such self-consistent resummation schemes, the mode-coupling approximation and the dynamical SCSAendnote3 .

## III Dynamical self-consistent approximations

To introduce the mode-coupling approximation, we first define the time-dependent correlation function $`C(𝐤,t,t^{})`$ and the associated response function $`G(𝐤,t,t^{})`$:

$`\delta (𝐤+𝐤^{})C(𝐤,t,t^{})`$ $`=`$ $`<\varphi _𝐤(t)\varphi _𝐤^{}(t^{})>,`$ (9)
$`\delta (𝐤+𝐤^{})G(𝐤,t,t^{})`$ $`=`$ $`<{\displaystyle \frac{\varphi _𝐤(t)}{\eta _𝐤^{}(t^{})}}>`$ (10)
$`=`$ $`{\displaystyle \frac{1}{2T}}<\varphi _𝐤(t)\eta _𝐤^{}(t^{})>.`$
As in the preceding section, we set $`T=1`$ in the following.

The perturbative expansion of $`C(𝐤,t,t^{})`$ and $`G(𝐤,t,t^{})`$ in powers of $`u`$ is more conveniently expressed by introducing the zeroth-order correlation and response functions,

$`C_0(𝐤,t,t^{})=`$ $`{\displaystyle \frac{1}{\mu (𝐤)}}\mathrm{exp}\left[\mu (𝐤)(tt^{})\right],`$ (11)
$`G_0(𝐤,t,t^{})=`$ $`\mathrm{exp}\left[\mu (𝐤)(tt^{})\right],`$ (12)
where

$$\mu (𝐤)=r_0+k^2+\frac{Q}{k^2},$$
(13)
and the two kernels $`\mathrm{\Sigma }(𝐤,t,t^{})`$ and $`D(𝐤,t,t^{})`$ defined though the standard Dyson equations,

$`C(𝐤,t,t^{})`$ $`=`$ $`{\displaystyle _0^t}𝑑t_1{\displaystyle _0^t^{}}𝑑t_2G(𝐤,t,t_1)\left[2\delta (t_1t_2)+D(𝐤,t_1,t_2)\right]G(𝐤,t^{},t_2),`$ (14)
$`G(𝐤,t,t^{})`$ $`=`$ $`G_0(𝐤,t,t^{})+{\displaystyle _t^{}^t}𝑑t_1{\displaystyle _t^{}^{t_1}}𝑑t_2G_0(𝐤,t,t_1)\mathrm{\Sigma }(𝐤,t_1,t_2)G(𝐤,t_2,t^{}).`$ (15)
The diagrammatic representation of the perturbative expansion and a detailed derivation of the mode-coupling approximation can be found in Ref.BCKM96 ; the only difference with the cases considered in BCKM96  is the presence of the frustration term $`Q/k^2`$ in the expression of $`\mu (k)`$, and we merely sketch here the main steps of the derivation.

The mode-coupling approximation amounts to expanding the kernels $`D(𝐤,t,t^{})`$ and $`\mathrm{\Sigma }(𝐤,t,t^{})`$ to second order in $`u`$ and replacing the bare (zeroth-order) functions $`C_0(𝐤,t,t^{})`$ and $`G_0(𝐤,t,t^{})`$ that appear in the resulting expressions by their renormalized counterparts $`C(𝐤,t,t^{})`$ and $`G(𝐤,t,t^{})`$. This leads toBCKM96 :

$`D(𝐤,t,t^{})6u^2{\displaystyle \frac{d^3𝐤_1}{(2\pi )^3}\frac{d^3𝐤_2}{(2\pi )^3}C(𝐤_1,t,t^{})C(𝐤_2,t,t^{})C(𝐤𝐤_1𝐤_2,t,t^{})},`$ (16)
$`\mathrm{\Sigma }(𝐤,t,t^{})18u^2{\displaystyle \frac{d^3𝐤_1}{(2\pi )^3}\frac{d^3𝐤_2}{(2\pi )^3}C(𝐤_1,t,t^{})C(𝐤_2,t,t^{})G(𝐤𝐤_1𝐤_2,t,t^{})}.`$ (17)
At the same time $`\mu (𝐤)`$ is renormalized to include the so-called tadpole diagrams, which replaces Eq. (13) by an expression similar to that obtained within the static Hartree approximation.

In this work, we are interested by the dynamical properties of the system at equilibrium: therefore, the fluctuation-dissipation theorem and the time-translation invariance apply, which reduces the dependence upon the two times $`t`$ and $`t^{}`$ to the mere dependence on the difference $`tt^{}`$ and gives

$$G(𝐤,t)=\mathrm{\Theta }(t)\frac{C(𝐤,t)}{t}$$
(18)
and

$$\mathrm{\Sigma }(𝐤,t)=\mathrm{\Theta }(t)\frac{D(𝐤,t)}{t},$$
(19)
where $`\mathrm{\Theta }(t)`$ is the Heaviside step function ($`\mathrm{\Theta }(t)`$ is equal to $`0`$ for $`t<0`$ and to $`1`$ pour $`t>0`$). Applying the operator $`G_0^1(𝐤)=\left[\mu (𝐤)+\frac{}{t}\right]`$ to both sides of Eq. (15) yields

$$\frac{G(𝐤,t)}{t}=\delta (t)\mu (𝐤)G(𝐤,t)+_0^t𝑑t^{}\mathrm{\Sigma }(𝐤,tt^{})G(𝐤,t^{}),$$
(20)
which, when combined with the time derivative of Eq. (18),

$$\frac{G(𝐤,t)}{t}=\delta (t)\left(\frac{C(𝐤,t)}{t}\right)_{t=0}\mathrm{\Theta }(t)\frac{^2C(𝐤,t)}{t^2},$$
(21)
gives

$$\left(\frac{C(𝐤,t)}{t}\right)_{t=0}=1.$$
(22)
For $`t>0`$, the equation for the response function thus reads

$$\frac{G(𝐤,t)}{t}=\mu (𝐤)G(𝐤,t)+_0^t𝑑t^{}\mathrm{\Sigma }(𝐤,tt^{})G(𝐤,t^{})$$
(23)
with the initial condition $`G(𝐤,t=0^+)=1`$. By Laplace transforming Eqs. (18), (19) and (23) and using the initial condition $`\left(\frac{C(𝐤,t)}{t}\right)_{t=0}=G(𝐤,t=0^+)=1`$, one finally obtains

$`1=`$ $`(z\widehat{C}(𝐤,z)+C(𝐤,t=0))\times `$
$`\times `$ $`\left(\mu (𝐤)izz\widehat{D}(𝐤,z)D(𝐤,t=0)\right),`$ (24)
where $`\widehat{C}(𝐤,z)=i_0^{\mathrm{}}𝑑te^{izt}C(𝐤,t)`$, and a similar expression holds for $`\widehat{D}(𝐤,z)`$. Going back to the time dependence leads to

$`{\displaystyle \frac{C(𝐤,t)}{t}}=`$ $`\left(\mu (𝐤)D(𝐤,t=0)\right)C(𝐤,t)`$
$``$ $`{\displaystyle _0^t}𝑑t^{}D(𝐤,tt^{}){\displaystyle \frac{C(𝐤,t^{})}{t^{}}}`$ (25)
with the initial condition $`C(𝐤,t=0)=\left(\mu (𝐤)D(𝐤,t=0)\right)^1`$ that follows from Eqs. (22) and (III).

The mode-coupling approximation finally results in the following self-consistent equation for the time-dependent correlation function at equilibrium:

$`{\displaystyle \frac{C(𝐤,t)}{t}}=`$ $`C(𝐤,t=0)^1C(𝐤,t)`$
$``$ $`{\displaystyle _0^t}𝑑t^{}D(𝐤,tt^{}){\displaystyle \frac{C(𝐤,t^{})}{t^{}}}`$ (26)
with

$`D(𝐤,t)=`$ $`6u^2{\displaystyle }{\displaystyle \frac{d^3𝐤_1}{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d^3𝐤_2}{(2\pi )^3}}C(𝐤_1,t)C(𝐤_2,t)\times `$
$`\times `$ $`C(𝐤𝐤_1𝐤_2,t).`$ (27)
Except for the absence of inertial term, $`\frac{^2C(𝐤,t)}{t^2}`$, in the purely relaxational dynamics associated with the Langevin equation and the cubic dependence of the memory kernel $`D(𝐤,t)`$ on the correlation function $`C(𝐤,t)`$, the above equations are similar to the mode-coupling equations used to describe the time-dependent density fluctuations in supercooled liquidsgotze91 ; they are also analogous to those derived for the mean-field spin glass with $`4`$-spin interactionsBCKM96 ; CHS93 .

The necessary input for solving the self-consistent Eqs. (III) and (III) is the knowledge of the equilibrium static correlation function $`C(𝐤,t=0)`$. Treating the statics and the dynamics of the system on an equal footing, as for instance done in the above derivation, leads to considering a mode-coupling approximation for the static correlation function, $`C(𝐤,t=0)=\left(\mu (𝐤)D(𝐤,t=0)\right)^1`$. However, we shall rather introduce more flexibility in the mode-coupling scheme (a flexibility that goes with the many ways to implement the self-consistency at the second order of the perturbative expansion) by allowing $`C(𝐤,t=0)`$ to be computed with several approximations, such as the Hartree approximation and the SCSA, that are a priori better behaved than the mode-coupling approximation as far as the static properties are concerned.

A somewhat refined resummation scheme is provided by the dynamical SCSA. (As mentioned in endnote2 , it consists in using an $`n`$-component vector field $`\varphi `$, resumming self-consistently all the diagrams of order $`1/n`$ in the large $`n`$ expansion, and, eventually, for the problem considered here, setting $`n`$ equal to $`1`$.) Details on the derivation of the approximate equation for the time-dependent correlation function $`C(𝐤,t)`$ can be found in Refs.BCKM96 ; MCJPB97 . A convenient way to proceed is to introduce a complex auxiliary field $`\sigma (𝐱)`$ such that the partition function, Eq. (3), can be rewritten $`Z=𝒟\varphi 𝒟\sigma e^{H[\varphi ,\sigma ]}`$ (here $`T=1`$) with

$`H[\varphi ,\sigma ]=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d^3𝐱\{(\varphi (𝐱))^2+r_0\varphi ^2(𝐱)\sigma ^2(𝐱)`$
$`+`$ $`\sqrt{2u}\sigma (𝐱)\varphi ^2(𝐱)\}`$
$`+`$ $`{\displaystyle \frac{Q}{8\pi }}{\displaystyle d^3𝐱d^3𝐱^{}\frac{\varphi (𝐱)\varphi (𝐱^{})}{|𝐱𝐱^{}|}}.`$ (28)
The field $`\sigma (𝐱)`$ is such that $`<\sigma (𝐱)>=\sqrt{u/2}<\varphi ^2(𝐱)>`$ and its connected pair correlation function $`C_\sigma (𝐱,𝐱^{})=<\sigma (𝐱)\sigma (𝐱^{})>_c`$ is equal to

$`C_\sigma (𝐱,𝐱^{})=`$ $`\delta (𝐱𝐱^{})`$
$`+`$ $`u/2(<\varphi ^2(𝐱)\varphi ^2(𝐱^{})><\varphi ^2(𝐱)><\varphi ^2(𝐱^{})>).`$ (29)
One can then apply to the dynamics of the coupled fields $`\sigma (𝐱)`$ and $`\varphi (𝐱)`$ a treatment similar to that sketched above. Defining the equilibrium time-dependent correlation function $`C_\sigma (𝐤,t)`$ via $`\delta (𝐤+𝐤^{})C_\sigma (𝐤,t)=<\sigma _𝐤(t)\sigma _𝐤^{}(0)>_c`$, the associated kernel $`D_\sigma (𝐤,t)`$ obtained through the Dyson equation (see above), and similar functions for the response properties, one can perform a mode-coupling approximation to the coupled dynamical equations for the fields $`\sigma (𝐱)`$ and $`\varphi (𝐱)`$. This leads to Eq. (III) with $`D(𝐤,t)`$ now given byMCJPB97

$$D(𝐤,t)=2u\frac{d^3𝐩}{(2\pi )^3}C_\sigma (𝐩,t)C(𝐤𝐩,t);$$
(30)
the auxiliary-field correlation function $`C_\sigma (𝐤,t)`$ is solution of the equation

$`{\displaystyle \frac{C_\sigma (𝐤,t)}{t}}=`$ $`C_\sigma (𝐤,t=0)^1C_\sigma (𝐤,t)`$
$``$ $`{\displaystyle _0^t}𝑑t^{}D_\sigma (𝐤,tt^{}){\displaystyle \frac{C_\sigma (𝐤,t^{})}{t^{}}}`$ (31)
with

$$D_\sigma (𝐤,t)=u\frac{d^3𝐪}{(2\pi )^3}C(𝐪,t)C(𝐤𝐪,t).$$
(32)
These equations are supplemented by the initial conditions $`C(𝐤,t=0)=\left(\mu (𝐤)D(𝐤,t=0)\right)^1`$ and $`C_\sigma (𝐤,t=0)=\left(1D_\sigma (𝐤,t=0)\right)^1`$ which are easily shown to be identical to the equilibrium, static SCSA equations first derived by BrayAJBa74 .

## IV Transition from ergodic to non-ergodic behavior

It is well known that mode-coupling and related approximations, when applied to glassforming systems, may lead to a dynamical singularitygotze91 ; BCKM96 . This latter corresponds to a transition from an ergodic to a non-ergodic behavior and is not associated with any thermodynamic equilibrium transition. For searching for such a singularity in the above equations, it is convenient to Laplace transform Eq. (III), which gives

$$\widehat{C}(𝐤,z)=\frac{C(𝐤,t=0)}{z{\displaystyle \frac{1}{C(𝐤,t=0)(i+\widehat{D}(𝐤,z))}}}.$$
(33)
An ergodicity-breaking transition is associated with the appearance of a non-zero value of the long-time limit of the correlation function, $`C(𝐤,t\mathrm{})`$; as a result, in the small-$`z`$ limit, $`\widehat{C}(𝐤,z)\frac{C(𝐤,t\mathrm{})}{z}`$ and $`\widehat{D}(𝐤,z)\frac{D(𝐤,t\mathrm{})}{z}`$, which when inserted in Eq. (33) leads to

$$D(𝐤,t\mathrm{})=\frac{C(𝐤,t\mathrm{})}{C(𝐤,t=0)(C(𝐤,t=0)C(𝐤,t\mathrm{}))}.$$
(34)
The kernel $`D(𝐤,t\mathrm{})`$ is obtained from Eq. (III) for the mode-coupling approximation and from Eqs. (30)-(32) for the dynamical SCSA.

Consider first the equation resulting from the mode-coupling approximation,

$$\frac{C(𝐤,t\mathrm{})}{C(𝐤,t=0)C(𝐤,t\mathrm{})}=6u^2C(𝐤,t=0)\frac{d^3𝐤_1}{(2\pi )^3}\frac{d^3𝐤_2}{(2\pi )^3}C(𝐤_1,t\mathrm{})C(𝐤_2,t\mathrm{})C(𝐤𝐤_1𝐤_2,t\mathrm{}).$$
(35)
Note that this actually represents a set of coupled equations for the various $`𝐤`$-modes. The necessary input for solving this equation is the knowledge of the equilibrium static correlation function $`C(𝐤,t=0)`$ (see discussion above). We consider here two standard approximations:

(i) the Hartree approximation, already presented in Eqs. (5) and (6), and

(ii) the SCSA, described in the previous section and leading to

$$C^1(𝐤,0)=\mu (𝐤)+2u\frac{d^3𝐪}{(2\pi )^3}\frac{C(𝐤𝐪,0)}{1+u\mathrm{\Pi }(𝐪)}.$$
(36)
where

$$\mu (𝐤)=r_0+k^2+\frac{Q}{k^2}+u\frac{d^3𝐪}{(2\pi )^3}C(𝐪,0)$$
(37)
and

$$\mathrm{\Pi }(𝐤)=\frac{d^3𝐪}{(2\pi )^3}C(𝐪,0)C(𝐤𝐪,0).$$
(38)
Note that if one neglects the term $`u\mathrm{\Pi }(𝐪)`$ in Eq. (36), one recovers the Hartree approximation.

In both cases, only the paramagnetic phase $`<\varphi _𝐤>=0`$, is considered. From now on, we take $`u=1`$ (recall that all momentum integrations are cut-off at $`\mathrm{\Lambda }=1`$).

Other approximations will be discussed in section VI.

We have solved the set of coupled equations, Eq. (35), by an iterative method. We find that with the above two approximations, the mode-coupling approach does lead to an ergodicity-breaking transition for the Coulomb frustrated model. When decreasing the temperature, i.e., the bare mass $`r_0`$, one reaches a point at which $`C(𝐤,t\mathrm{})`$ discontinuously jumps to a non-zero value. The transition temperature increases as frustration decreases; it seems to reach continuously, when $`Q0`$, the (equilibrium) critical temperature of the unfrustrated system, i.e., of the standard $`\varphi ^4`$ theory in either the Hartree approximation or the SCSA. This behavior is illustrated in Fig. 2. One observes a discrepancy between the results obtained with the two different approximations (i) and (ii) for $`C(𝐤,t=0)`$, but it stays within reasonable bounds: the relative difference is about $`20\%`$ or less. As can be seen from Fig. 3, the two approximations predict very similar correlation functions, both at equilibrium ($`t=0`$) and in the non-ergodic state $`(t\mathrm{})`$ when the temperature is at (or just below) the dynamical transition. Roughly speaking, this latter takes place when the maximum of the equilibrium correlation function $`C(𝐤,t=0)`$, a maximum that occurs for $`|𝐤|k_m=Q^{1/4}`$, reaches a given, $`Q`$-dependent value: this is illustrated for $`Q=0.1`$ in Fig. 4 where the dynamical transition occurs when $`\mathrm{Max}_𝐤\{C(𝐤,0)\}40`$. As the maximum increases when $`Q`$ decreases slightly more rapidly with the SCSA than with the Hartree approximation, the former predicts a somewhat higher transition temperature than the latter: see Fig. 4. Note that the fact that the ergodicity-breaking transition is driven by the maximum of the $`𝐤`$-dependent equilibrium correlation function is well established in the context of mode-coupling approachesgotze91 . As illustrated in Fig. 5, the transition takes place at the same temperature for all $`𝐤`$-modes.

The transition form ergodic to non-erogodic dynamical behavior can also been studied within the dynamical SCSA. The corresponding equations to be solved are given above and in Appendix B. A dynamical transition is indeed found, and the frustration-dependence of the transition temperature is shown in Fig. 2. The predicted transition line is not much different from those obtained with the above mode-coupling approximations. We show in Appendix B that the expressions for the correlation function $`C(𝐤,t)`$ for $`t=0`$ and $`t\mathrm{}`$ derived within the dynamical SCSA when ergodicity is broken are identical to those obtained in Refs.ScWo00 ; SW201  by using the purely thermodynamic analysis based on the replica formalism and the static SCSA. As a result, the dynamical SCSA predicts that the dynamics looses ergodicity precisely at the point at which an exponentially large number of metastable states occurs in the Schmalian-Wolynes treatment. Below this temperature, ergodicity is broken: the fluctuation-dissipation theorem and the time-translation invariance no longer apply, and Eqs. (III)-(30) should be generalized to describe the evolution of two-time correlation and response functions and the associated aging behaviorBCKM96 .

Finally, it is instructive to compare the location of the mode-coupling-like dynamical transition with that of the equilibrium, thermodynamic transition discussed in section II. The dynamical transition occurs at a temperature that is lower than the critical temperature of the unfrustrated system, a temperature that was shown in the Monte Carlo study of Ref.GrTV01  to mark the onset of fragile glassforming behavior; it seems to occur at a temperature close to that of the fluctuation-induced first-order transition from the paramagnetic to the modulated phases. This is illustrated in Fig. 6 where we display the first-order transition obtained within the Hartree approximation (see section II and Fig. 1) and the dynamical transition obtained within the mode-coupling approximation supplemented by the static Hartree approximation (as discussed above, the other predictions are quite close to this latter).

Actually at small $`Q`$’s the dynamical transition appears even below the temperature of the equilibrium first-order transition: in such a case, the dynamical transition takes place in the supercooled paramagnetic (“liquid”) regime, a regime that appears because of the first-order nature of the transition to the modulated phases and that can be described by the Hartree approximation. As discussed in section II, the paramagnetic phase is (meta)-stable at all finite temperatures within this approximation, and the equilibrium correlation length is therefore finite in all the region where the dynamics are studied.

## V Evolution with time of the correlation function

We have also solved the full set of coupled equations describing the time evolution of the equilibrium correlation function $`C(𝐤,t)`$ in the mode-coupling approximation, Eqs. (III)- (III). (The algorithm is described in RefsFGHL91 ; K00 .) For the input quantity, $`C(𝐤,t=0)`$, we have used the Hartree approximation. The results are shown in Fig. 7 for the time-dependent correlation function at a momentum $`k_{\mathrm{max}}k_m=Q^{1/4}`$ that corresponds to the maximum value of the function; curves for the frustration parameter $`Q=0.1`$ and several temperatures (i.e., several values of the bare mass $`r_0`$) are shown. One observes a behavior typical of the mode-coupling equations with a so-called $`B`$-type transitiongotze91  as those used to describe glassforming liquidsgotze91  and those describing the dynamics of a class of mean-field generalized spin glassesBCKM96 .

At high temperature, the correlation function decays in one step, but as temperature is lowered a second relaxation step appears, that becomes slower and slower so that a plateau develops between the two relaxation steps. When temperature is further decreased, one reaches a point at which the slow (“$`\alpha `$”) relaxation time diverges. The correlation function no longer decays to zero, but stays at the plateau value. Below this point, ergodicity is broken and the mode-coupling equations derived under the condition of equilibrium (with the fluctuation-dissipation theorem and the time-translation invariance) are no longer valid.

In the vicinity of the dynamical transition, various scaling laws are observed, and the slow (“$`\alpha `$”) relaxation time diverges as a power law,

$$\tau _\alpha (Q,T)(TT_c)^\gamma ,$$
(39)
where the exponent $`\gamma 1.851.89`$ is weakly dependent on the frustration parameter, provided $`Q>0`$. (For $`Q=0`$, the system shows standard critical slowing down with $`\tau (T)(TT_c)^{z\nu }`$, where $`z`$ is the dynamical exponent and $`\nu `$ the (static) correlation length exponentHH77 .) For illustration, we have plotted the logarithm of the $`\alpha `$ relaxation time versus temperature $`r_0`$ for two different frustrations in Fig. 8.

## VI Sensitivity of the results to the approximation scheme

We have already mentioned (see Ref. endnote3 ) that enough of the non-linearities of the original dynamical equation or of the equations in replica space must be kept in any approximate treatment in order to find non-trivial phenomena such as ergodicity breaking and appearance of an exponentially large number of metastable states. For this reason, the dynamical as well as the replica-space Hartree approximations are unable to generate such phenomena. One must therefore consider improved resummation schemes such as the mode-coupling approximation and the SCSAendnote4 .

The additional point we would like to make here is that even in the mode-coupling approximation, the results are somewhat dependent upon the supplementary approximation which is made to describe the static properties of the system entering as an input in the dynamical equation. It is well known, and was recalled above, that the location, or even the existence of an ergodicity-breaking transition is sensitive to the amplitude of the peak in the equilibrium (static) correlation function $`C(𝐤,t=0)`$. We have seen that the Hartree and SCSA give slightly different, but compatible results. If one uses instead a somewhat less renormalized version of the SCSA with $`\mu (𝐤)`$ in Eq. (37) now defined with the Hartree and not the full correlation function, i.e.,

$$\mu (𝐤)=r_0+k^2+\frac{Q}{k^2}+u\frac{d^3𝐪}{(2\pi )^3}\frac{1}{\mu (𝐪)},$$
(40)
a different behavior is obtained. As shown in Fig. 9, the maximum of the static correlation function appears to saturate, as one lowers the temperature, to a value that is too small to trigger a breaking of ergodicity.

Finally, considering the static analog of the mode-coupling approximation to compute $`C(𝐤,t=0)`$ (see Eq. (III) and below), i.e.,

$`C^1(𝐤,t=0)=`$ $`r_0+k^2+{\displaystyle \frac{Q}{k^2}}+3u{\displaystyle \frac{d^3𝐪}{(2\pi )^3}C(𝐪,t=0)}`$
$``$ $`D(𝐤,t=0).`$ (41)
where $`D(𝐤,t=0)`$ is given by Eq. (III) leads to a situation in which the limit of stability (spinodal) of the paramagnetic phase is reached at a finite temperature, before the occurrence of an ergodicity-breaking transition.

The validity of these various approximations should of course be checked by performing a computer simulation of the model. However, one can tentatively conclude from the above exercise that despite the formal similarity between the dynamics and the statics that comes from using the Martin-Siggia-Rose functional formalismMSR78 , different levels of approximation may be required to describe the dynamical and the static properties of the Coulomb frustrated model.

## VII Conclusion

We have studied the Langevin dynamics of the soft-spin, continuum version of the Coulomb frustrated Ising ferromagnet. By using the dynamical mode-coupling approximation, coupled with reasonable approximations for describing the equilibrium static correlation function, and the dynamical self-consistent screening approximation, we have found that the system’s dynamics display a transition from ergodic to non-ergodic behavior, similar to that obtained in the idealized mode-coupling theory of glassforming liquidsgotze91  and in the mean-field generalized spin-glasses with one-step replica symmetry breakingKT87 ; BCKM96 . This transition occurs in the paramagnetic phase, either in the stable or the metastable region. It is related to the emergence of an exponentially large number of metastable states found by a purely static replica approach: the system looses ergodicity because it gets trapped in free-energy minima separated from each other by infinite barriers. This whole description, as can be inferred from the very nature of the approximations that amount to partial resummations of perturbative expansions, has a mean-field character: thermally activated processes are completely ignored. The predicted singularity is “avoided” in the true dynamics of the system and it remains to be seen, e.g., in computer simulations of Coulomb frustrated models, what signatures may still be observed in the time evolution of the correlation function. As for describing the activated processes, other, non-perturbative approaches, such as the phenomenological frustration-limited domainKKZNT95  and entropic-droplet picturesKTW89 , must be used.

## Appendix A Fluctuation-induced first order transition

In this appendix, we calculate the temperature of the equilibrium transition between the paramagnetic phase and phases with spatially modulated order within the self-consistent Hartree approximation. The derivation given below closely follows Brazovskii’s original treatmentB75 .

The starting point is the Hamiltonian in Eq. (2) augmented by the introduction of spatially varying external fields $`h_𝐤`$ that are linearly coupled to the scalar field $`\varphi _𝐤`$. As a result, $`\varphi _𝐤`$ is now the sum of an average component, $`m_𝐤=<\varphi _𝐤>`$, and a fluctuation $`\psi _𝐤=\varphi _𝐤m_𝐤`$. The self-consistent Hartree approximation is then equivalent to a gaussian variational approximation for the fluctuations. The resulting equation of state reads

$`h_𝐤=`$ $`\left(r_0+k^2+{\displaystyle \frac{Q}{k^2}}\right)m_𝐤+u{\displaystyle \frac{d^3𝐤_1}{(2\pi )^3}\frac{d^3𝐤_2}{(2\pi )^3}\left(m_{𝐤_\mathrm{𝟏}}m_{𝐤_\mathrm{𝟐}}+3C(𝐤_\mathrm{𝟏},𝐤_\mathrm{𝟐})\right)m_{𝐤𝐤_\mathrm{𝟏}𝐤_\mathrm{𝟐}}},`$ (42)
where the connected correlation function $`C(𝐤,𝐤^{})=<\psi _𝐤\psi _𝐤^{}>`$, is obtained self-consistently by solving

$$C^1(𝐤,𝐤^{})=\left(r_0+k^2+\frac{Q}{k^2}\right)\delta (𝐤+𝐤^{})+3u\frac{d^3𝐪}{(2\pi )^3}\left(m_𝐪m_{𝐤+𝐤^{}𝐪}+C(𝐪,𝐤+𝐤^{}𝐪)\right),$$
(43)
together with the inversion formula

$$\frac{d^3𝐪}{(2\pi )^3}C^1(𝐤,𝐪)C(𝐪,𝐤^{})=\delta (𝐤^{}𝐤).$$
(44)
In the paramagnetic phase when all $`h_𝐤`$’s, and therefore $`m_𝐤`$’s, are equal to zero, Eqs. (43) and (44) reduce to Eq. (5) with $`C(𝐤,𝐤^{})=C(𝐤)\delta (𝐤+𝐤^{})`$.

In the vicinity of the transition between paramagnetic and modulated phases and for small enough frustration $`Q`$ ($`Q<<1`$), the modulated order is one-dimensional and characterized by a wave-vector $`𝐤_m`$ with $`|𝐤_m|=k_m=Q^{1/4}`$. It is then sufficient to consider

$$h_𝐤=\stackrel{~}{h}\left(\delta (𝐤𝐤_𝐦)+\delta (𝐤+𝐤_𝐦)\right)$$
(45)
and

$$m_𝐤=\stackrel{~}{m}\left(\delta (𝐤𝐤_𝐦)+\delta (𝐤+𝐤_𝐦)\right).$$
(46)
In this region, the fluctuations of wave-vector $`𝐤`$ with $`|𝐤|=k_m`$ are dominant, and BrazovskiiB75  has shown that the effect of the off-diagonal terms with $`𝐤𝐤^{}`$ could be neglected in the correlation function. As a result,

$$C(𝐤,𝐤^{})\frac{\delta (𝐤+𝐤^{})}{r+k^2+\frac{Q}{k^2}},$$
(47)
where the renormalized mass in a phase characterized by Eqs. (45,46), is given by:

$`r=`$ $`r_0+3u{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\left(\frac{1}{r+k^2+\frac{Q}{k^2}}+|m_𝐤|^2\right)}`$
$`=`$ $`r_0+3u{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\frac{1}{r+k^2+\frac{Q}{k^2}}}+6u|\stackrel{~}{m}|^2.`$ (48)
By introducing Eqs. (45)- (A) in Eq. (42) and recalling that $`k_m=Q^{1/4}`$, one obtains the following equation of state:

$`\stackrel{~}{h}=`$ $`(r_0+2Q^{1/2}`$
$`+`$ $`3u{\displaystyle }{\displaystyle \frac{d^3𝐤}{(2\pi )^3}}{\displaystyle \frac{1}{r+k^2+\frac{Q}{k^2}}}+3u|\stackrel{~}{m}|^2)\stackrel{~}{m}`$
$`=`$ $`\left(r+2Q^{1/2}3u|\stackrel{~}{m}|^2\right)\stackrel{~}{m}.`$ (49)
In the Hartree approximation, and below some temperature, there is a coexistence of the paramagnetic phase and the modulated phase. In zero field ($`\stackrel{~}{h}=0`$), the former is characterized by $`\stackrel{~}{m}=0`$ and the latter by $`\stackrel{~}{m}0`$, where $`\stackrel{~}{m}`$ is solution of $`\left(r+2Q^{1/2}3u|\stackrel{~}{m}|^2\right)=0`$. The transition point, which is then associated with a first-order transition, is obtained as the temperature at which the free-energies of the two phases are equal. Following BrazovskiiB75 , it is convenient to calculate directly the free-energy difference $`\mathrm{\Delta }F(r_0)`$ between the modulated ($`\stackrel{~}{m}0`$) and the paramagnetic ($`\stackrel{~}{m}=0`$) phases at a given temperature $`r_0`$ from the following expression:

$`\mathrm{\Delta }F=`$ $`{\displaystyle _0^{\stackrel{~}{m}}}𝑑\stackrel{~}{m}^{}{\displaystyle \frac{F}{\stackrel{~}{m}^{}}}=2{\displaystyle _0^{\stackrel{~}{m}}}𝑑\stackrel{~}{m}^{}\stackrel{~}{h}(\stackrel{~}{m}^{})`$ (50)
where $`\stackrel{~}{h}_{\stackrel{~}{m}^{}}`$ is given by Eq. (A). One can change the integration variable from $`\stackrel{~}{m}^{}`$ to $`r^{}`$ with $`r^{}(\stackrel{~}{m}^{})`$ solution of Eq. (A). After some algebra, Eq. (50) can be recast as

$`u\mathrm{\Delta }F=`$ $`{\displaystyle _{r(\stackrel{~}{m}=0)}^{r(\stackrel{~}{m})}}𝑑\stackrel{~}{r}^{}\left({\displaystyle \frac{r^{}+r_0}{2}}+2Q^{1/2}+{\displaystyle \frac{3u}{4\pi ^2}}{\displaystyle 𝑑k\frac{k^2}{r+k^2+\frac{Q}{k^2}}}\right)\left({\displaystyle \frac{1}{6}}+{\displaystyle \frac{u}{4\pi ^2}}{\displaystyle 𝑑k\frac{k^2}{\left(r^{}+k^2+\frac{Q}{k^2}\right)^2}}\right),`$ (51)
where $`r(\stackrel{~}{m}=0)`$ is solution of Eq. (A) with $`\stackrel{~}{m}=0`$, i.e., of Eq. (6), and $`r(\stackrel{~}{m})`$ and $`\stackrel{~}{m}`$ are solutions of the two coupled equations, Eq. (A) and $`\left(r+2Q^{1/2}3u|\stackrel{~}{m}|^2\right)=0`$. By solving Eq. (51) numerically for several values of $`Q<<1`$ (and for $`u=1`$), we have found that the sign of $`\mathrm{\Delta }F`$ changes at a finite value of $`r_0`$ that marks the first-order transition between paramagnetic and modulated phases. The result is shown in Fig. 1. The transition being second-order in the mean-field approximation, it is then driven first-order by the fluctuations.

## Appendix B Ergodicity breaking in the Dynamical SCSA

One can see from Eqs. (30)- (32) that ergodicity breaking requires that both $`C(𝐤,t)`$ and $`C_\sigma (𝐤,t)`$, and as a consequence $`D(𝐤,t)`$ and $`D_\sigma (𝐤,t)`$, go to non-zero values in the limit $`t\mathrm{}`$. From Eq. (34)), one obtains

$$C(𝐤,t\mathrm{})=\frac{D(𝐤,t\mathrm{})C(𝐤,0)^2}{1+D(𝐤,t\mathrm{})C(𝐤,0)}.$$
(52)
A similar expression can be derived for $`C_\sigma (𝐤,t\mathrm{})`$ by first Laplace transforming Eq. (III),

$$\widehat{C}_\sigma (𝐤,z)=\frac{C_\sigma (𝐤,0)}{z{\displaystyle \frac{1}{C_\sigma (𝐤,0)(i+\widehat{D}_\sigma (𝐤,z))}}},$$
(53)
and by looking for the dominant behavior in the small-$`z`$ limit, $`\widehat{C}_\sigma (𝐤,z)C_\sigma (𝐤,t\mathrm{})/z`$, $`\widehat{D}_\sigma (𝐤,z)D_\sigma (𝐤,t\mathrm{})/z`$; one finally gets

$$C_\sigma (𝐤,t\mathrm{})=\frac{D_\sigma (𝐤,t\mathrm{})C_\sigma (𝐤,0)^2}{1D_\sigma (𝐤,t\mathrm{})C_\sigma (𝐤,0)}.$$
(54)
By introducing the time-dependent polarization

$$\mathrm{\Pi }(𝐤,t)=\frac{d^3𝐪}{(2\pi )^3}C(𝐪,t)C(𝐤𝐪,t),$$
(55)
one can express the memory kernel $`D_\sigma (𝐤,t)`$ given in Eq. (32) as

$$D_\sigma (𝐤,t)=u\mathrm{\Pi }(𝐤,t),$$
(56)
so that the $`t=0`$ and $`t\mathrm{}`$ values of $`C_\sigma (𝐤,t)`$ (given below Eq. (32) and in Eq. (54), respectively) can be written as

$`C_\sigma (𝐤,0)=`$ $`{\displaystyle \frac{1}{1+u\mathrm{\Pi }(𝐤,0)}}`$ (57)
$`C_\sigma (𝐤,t\mathrm{})=`$ $`{\displaystyle \frac{u\mathrm{\Pi }(𝐤,t\mathrm{})C_\sigma (𝐤,0)^2}{1+u\mathrm{\Pi }(𝐤,t\mathrm{})C_\sigma (𝐤,0)}}.`$ (58)
Recalling that $`C(𝐤,0)=(\mu (𝐤)D(𝐤,0))^1`$ with $`\mu (𝐤)=r_0+k^2+Q/k^2+u\frac{d^3𝐪}{(2\pi )^3}C(𝐪,0)`$ and that $`D(𝐤,t)`$ is given by Eq. (III), one obtains with Eqs. (52), (56) and (57) a closed set of equations that determines the non-ergodicity parameter $`C(𝐤,t\mathrm{})`$. If one changes the notations from $`C(𝐤,0)`$ and $`C(𝐤,t\mathrm{})`$ to $`𝒢(𝐤)`$ and $`(𝐤)`$, from $`D(𝐤,0)`$ and $`D(𝐤,t\mathrm{})`$ to $`\mathrm{\Sigma }(𝐤)`$ and $`\mathrm{\Sigma }_\sigma (𝐤)`$, from $`uC(𝐤,0)`$ and $`uC_\sigma (𝐤,t\mathrm{})`$ to $`𝒟_𝒢(𝐤)`$ and $`𝒟_{}(𝐤)`$, one can easily check that the above equations are identical to those obtained in RefsScWo00 ; SW201  with the replica formalism and the static SCSA.


format: text

text:
# Microstructure and Velocity of Field-Driven SOS Interfaces: Analytic Approximations and Numerical Results

## I Introduction

The motion of surfaces and interfaces plays a central role in many scientific and technological disciplines. In particular, the dynamics of interfaces such as phase- and grain boundaries in solid materials MEND00  and domain walls in magnets RIKVturkey  and ferroelectrics heavily influence both dynamic and static materials properties. Among interfaces characteristic of two-dimensional systems are steps on crystal surfaces Surfstep , domain walls in thin magnetic and dielectric films RIKVturkey , and boundaries between different types of vegetation such as savanna and rainforest RAINF .

An enormous amount of work in recent years has been devoted to the dynamics and structure of moving and growing interfaces BARA95 ; MEAK98 . However, despite the fact that many important interface properties, such as mobility and catalytic and chemical activity, are largely determined by the microscopic interface structure, most of this effort has been concentrated on large-scale scaling properties. Although the detailed atomistic mechanisms by which interfaces move are often not known, useful understanding can be obtained from stochastic models in which the motion occurs through random nucleation and migration of local topological features such as kinks and steps MEND00 . It is therefore important to gain better insight into how, for different stochastic dynamics, the driving force (such as an applied magnetic or electric field, a chemical-potential gradient, or the amount of rainfall in the case of models of vegetation distribution) may alter the microscopic structure of the interface, thereby leading to a highly nonlinear velocity response.

In a previous paper RIKV00  we introduced a dynamic mean-field approximation for the microstructure of an interface in a two-dimensional kinetic Ising ferromagnet with a single-spin-flip Glauber dynamic KAWA72 ; LAND00 , driven by an applied field DEVI92 ; SPOH93 . This model is directly applicable to many magnetic and ferroelectric systems and other cases where the interface dynamics are not inhibited by coupling to a conserved field RAMO99 ; NOVO99B . Based on the resulting local interface structure, we obtained a nonlinear-response approximation for the steady-state propagation velocity, which was shown to be in good agreement with dynamic Monte Carlo (MC) simulations for a wide range of fields and temperatures. However, since the approximation was based on the Burton-Cabrera-Frank (BCF) unrestricted solid-on-solid (SOS) model BURT51 , the overhangs and bubbles in the Ising interface were handled in an uncontrolled way. Here we therefore consider the performance of our approximation for the unrestricted SOS model, so that overhangs and bubbles are absent at all times by definition of the model. In particular, we obtain the surface velocity under the Glauber dynamic as a function of applied field and temperature, as well as its anisotropy for tilt angles between $`0^{}`$ and $`45^{}`$.

In a recent paper RIKV01  we showed that the microscopic interface structure, and thus the mobility, can depend dramatically on the details of the dynamics. The most significant difference is between dynamics in which the transition probabilities of the local spins factorize into one part that depends only on the change in interaction energy due to the transition and one that depends only on the change in the field energy (soft dynamics MARR99 ), and dynamics that cannot be factorized in this way (hard dynamics MARR99 ). For soft dynamics, the interface structure for all values of the field remains the same as in equilibrium at zero field RIKV01 , so that the mobility can be calculated exactly at the level of linear response. In the present paper we concentrate on the standard Glauber dynamic KAWA72 ; LAND00 , which belongs to the class of hard dynamics and thus leads to more complicated and interesting behavior.

Both the driven Ising and SOS surfaces belong to the Kardar-Parisi-Zhang (KPZ) dynamic universality class BARA95 ; KARD86 , in which the macroscopic, stationary distribution for flat, moving interfaces is Gaussian, corresponding to a random walk with independent increments. Nevertheless, the step heights in several discrete models in this class are correlated at short distances NEER97 ; KORN00 . In the mean-field approximation developed here, these short-range correlations are ignored. The resulting discrepancies, which are minor, are elucidated by comparison with MC simulations.

The remainder of this paper is organized as follows. In Sec. II we introduce the SOS interface model and derive a linear-response approximation for its velocity in a nonzero external field. In Sec. III we develop a mean-field approximation for the time evolution of the single-step probability density function (pdf), as well as for its stationary form. The latter enables us to extend the approximation for the interface structure and velocity to a nonlinear-response level. The analytical approximations are compared with dynamic MC simulations in Sec. IV. In Sec. IV.1 we numerically solve the mean-field equation of motion for the single-step pdf and compare the resulting values of the time-dependent average interface step height with MC simulations. In Sec. IV.2 we compare the simulated stationary single-step pdfs with the theoretical predictions. In Sec. IV.3 we compare the simulated stationary interface velocity with the theoretical predictions for various values of applied field, temperature, and interface tilt angle. In Sec. IV.4 we compare simulations and analytical predictions for the detailed stationary interface structure, including the asymmetry of the simulated interface at nonzero fields. A summary and conclusions are found in Sec. V.

## II Model and Dynamics

The original BCF SOS model considers an interface in a lattice gas or $`S=1/2`$ Ising system on a square lattice of unit lattice constant as a series of integer-valued steps $`\delta (x)`$, parallel to the $`y`$ axis. The interface height is thus a single-valued integer function of the $`x`$-coordinate, with steps at integer values of $`x`$. The configuration energy is given by the nearest-neighbor Ising Hamiltonian with anisotropic, ferromagnetic interactions $`J_x`$ and $`J_y`$ in the $`x`$ and $`y`$ direction, respectively:

$$=\underset{x,y}{}s_{x,y}\left(J_xs_{x+1,y}+J_ys_{x,y+1}+H\right).$$
(1)
Here the two states at lattice site $`(x,y)`$ are denoted by $`s_{x,y}=\pm 1`$, and $`_{x,y}`$ runs over all sites. The quantity $`H`$ is the applied field, and the interface is introduced by fixing $`s_{x,y}=+1`$ and $``$1 for large negative and positive $`y`$, respectively. Without loss of generality we take $`H0`$, such that the interface on average moves in the positive $`y`$ direction.

A single-spin-flip (nonconservative) dynamic which satisfies detailed balance, such as the Metropolis or Glauber algorithms KAWA72 ; LAND00 , ensures the approach to equilibrium, which in this case is a uniformly positive phase with the interface pushed off to positive infinity. Such algorithms are defined by a single-spin transition probability, $`W(s_{x,y}s_{x,y})=W(\beta \mathrm{\Delta }E,\beta U)`$. Here $`\beta `$ is the inverse of the temperature $`T`$ (Boltzmann’s constant is taken as unity), $`\mathrm{\Delta }E`$ is the energy change corresponding to a successful spin flip, and the optional parameter $`U`$ is an energy barrier between the two states that enters into Arrhenius-type stochastic dynamics KANG89  and other dynamics that include a transition state TDA . The detailed-balance condition is expressed as $`W(\beta \mathrm{\Delta }E,\beta U)/W(\beta \mathrm{\Delta }E,\beta U)=e^{\beta \mathrm{\Delta }E}`$, where the right-hand side is independent of $`U`$. (In the case of soft dynamics, the detailed-balance condition is satisfied independently for the two parts of $`W`$.) In order to preserve the SOS configuration at all times, flips are allowed only at sites which have exactly one broken bond in the $`y`$ direction.

With the Ising Hamiltonian there are only a finite number of different values of $`\mathrm{\Delta }E`$, and the spins can therefore be divided into classes SPOH93 ; BORT75 ; NOVO95A , labeled by the spin value $`s`$ and the number of broken bonds between the spin and its nearest neighbors in the $`x`$ and $`y`$ direction, $`j`$ and $`k`$, respectively. The ten spin classes consistent with the SOS model are denoted $`jks`$ with $`j\{0,1,2\}`$ and $`k\{0,1\}`$. They are listed in Table 1.

In this paper, as in Ref. RIKV00 , we use the standard discrete-time Glauber dynamic with the transition probability KAWA72 ; LAND00

$$W_\mathrm{G}(s_{x,y}s_{x,y})=\left[1+e^{\beta \mathrm{\Delta }E}\right]^1.$$
(2)
Time is measured in MC steps per spin (MCSS).

In the BCF SOS model the heights of the individual steps are assumed to be statistically independent and identically distributed. The step-height probability density function (pdf) is given by the interaction energy corresponding to the $`|\delta (x)|`$ broken $`J_x`$-bonds between spins in the columns centered at $`(x1/2)`$ and $`(x+1/2)`$ as

$$p[\delta (x)]=Z(\varphi )^1X^{|\delta (x)|}e^{\gamma (\varphi )\delta (x)}.$$
(3)
The Boltzmann factor $`X=e^{2\beta J_x}`$ determines the width of the pdf, and $`\gamma (\varphi )`$ is a Lagrange multiplier which maintains the mean step height at an $`x`$-independent value, $`\delta (x)=\mathrm{tan}\varphi `$, where $`\varphi `$ is the average angle between the interface and the $`x`$ axis. The partition function is

$$Z(\varphi )=\underset{\delta =\mathrm{}}{\overset{+\mathrm{}}{}}X^{|\delta |}e^{\gamma (\varphi )\delta }=\frac{1X^2}{12X\mathrm{cosh}\gamma (\varphi )+X^2}$$
(4)
with $`\gamma (\varphi )`$ given by

$$e^{\gamma (\varphi )}=\frac{\left(1+X^2\right)\mathrm{tan}\varphi +\left[\left(1X^2\right)^2\mathrm{tan}^2\varphi +4X^2\right]^{1/2}}{2X\left(1+\mathrm{tan}\varphi \right)}$$
(5)
(see details in Ref. RIKV00 ). For $`\varphi =0`$, $`\gamma (\varphi )=0`$, and the partition function simplifies to

$$Z(0)=(1+X)/(1X).$$
(6)
The mean spin-class populations, $`n(jks)`$, are all obtained from the product of the independent pdfs for $`\delta (x)`$ and $`\delta (x`$+1). Symmetry of $`p[\delta (x)]`$ under the transformation $`(x,\varphi ,\delta )(x,\varphi ,\delta )`$ ensures that $`n(jk)=n(jk+)`$ for all $`j`$ and $`k`$. Numerical results illustrating the breakdown of this up/down symmetry for large $`|H|`$ are discussed in Sec. IV.4. As discussed in Ref. RIKV00 , calculation of the individual class populations is straightforward but somewhat tedious, especially for nonzero $`\varphi `$. The final results are summarized in Table 2.

Whenever a spin flips from $``$1 to +1, the corresponding column of the interface advances by one lattice constant in the $`y`$ direction. Conversely, the column recedes by one lattice constant when a spin flips from +1 to $``$1. The corresponding energy changes are given in the third column in Table 1. Since the spin-class populations on both sides of the interface are equal in this approximation, the contribution to the mean velocity in the $`y`$ direction from sites in the classes $`jk`$ and $`jk+`$ becomes

$$v_y(jk)=W(\beta \mathrm{\Delta }E(jk),\beta U)W(\beta \mathrm{\Delta }E(jk+),\beta U).$$
(7)
The results corresponding to the hard Glauber transition probabilities used here, Eq. (2), are listed in the last column of Table 2. The mean propagation velocity perpendicular to the interface becomes

$$v_{}(T,H,\varphi )=\mathrm{cos}\varphi \underset{j,k}{}n(jks)v_y(jk),$$
(8)
where the sum runs over the classes included in the tables. While the general result is cumbersome if written out in detail, the special case of $`\varphi =0`$ leads to a compact formula:

$`v_{}(T,H,0)`$ $`=`$ $`{\displaystyle \frac{\mathrm{tanh}(\beta H)}{(1+X)^2}}\left\{2X+{\displaystyle \frac{1+X^2}{1+\left[\frac{\mathrm{sinh}(2\beta J_x)}{\mathrm{cosh}(\beta H)}\right]^2}}\right\}.`$ (9)
It was shown in Ref. RIKV00  that Eq. (8) reduces to the results for the single-step DEVI92 ; SPOH93 ; MEAK86 ; PLIS87  and the polynuclear growth DEVI92 ; KRUG89 ; KERT89  models at low temperatures for strong and weak fields, respectively.

## III Nonlinear-Response

With $`X=e^{2\beta J_x}`$, the results in Table 2 correspond to a linear-response approximation. In Ref. RIKV00  we developed a mean-field approximation leading to a field-dependent $`X(T,H)`$, based on a detailed-balance argument for the stationary state. Here we show that this detailed-balance relation follows naturally from a dynamic mean-field approximation for the equation of motion for the single-step pdf during the approach to the stationary state.

We denote the total transition probability for the height of the single step at $`x`$ to change from $`\delta (x)`$ to $`\delta (x)\pm 1`$ as $`𝒲[\delta (x)\delta (x)\pm 1]`$. In terms of these transition probabilities, the equation of motion for the single-step pdf, $`p[\delta (x),t]`$ becomes

$`{\displaystyle \frac{dp[\delta (x),t]}{dt}}`$ $`=`$ $`p[\delta (x),t]\{𝒲[\delta (x)\delta (x)1]`$ (10)
$`+𝒲[\delta (x)\delta (x)+1]\}`$
$`+p[\delta (x)+1]𝒲[\delta (x)+1\delta (x)]`$
$`+p[\delta (x)1]𝒲[\delta (x)1\delta (x)],`$
where the coupling to the joint multistep probability density is hidden in the single-step transition rates $`𝒲`$.

To obtain an approximation for $`𝒲[\delta \delta \pm 1]`$, we employ the same mean-field assumption of independent steps as in equilibrium. For $`\delta (x)1`$ to increase to $`\delta (x)+1`$, either the spin in front of the interface at $`x`$+1/2 can flip from $``$1 to +1, or the spin behind the interface at $`x`$$``$1/2 can flip from +1 to $``$1. In each of these cases, $`\mathrm{\Delta }E`$ can have two different values, depending on the value of $`\delta (x`$+1) and $`\delta (x`$$``$1), respectively. The same argument holds also for the reverse transition, $`\delta (x)+1\delta (x)`$. The energy changes and corresponding conditions on $`\delta (x`$+1) and $`\delta (x`$$``$1), which are shown in Fig. 1, yield

$`𝒲[\delta \delta \pm 1,t]`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{[W(2H)+W(+2H)]\mathrm{\Pi }_\pm (t)`$
$`+[W(2H+4J_x)+W(+2H+4J_x)][1\mathrm{\Pi }_\pm (t)]\}`$
$`𝒲[\delta \pm 1\delta ,t]`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{[W(2H)+W(+2H)][1\mathrm{\Pi }_{}(t)]`$ (11)
$`+[W(2H4J_x)+W(+2H4J_x)]\mathrm{\Pi }_{}(t)\},`$
where the upper signs refer to $`\delta +1`$ and the lower signs to $`\delta 1`$, $`\mathrm{\Pi }_+(t)=_{\delta =+1}^+\mathrm{}p[\delta ,t]`$, and $`\mathrm{\Pi }_{}(t)=_{\delta =1}^{\mathrm{}}p[\delta ,t]`$. For simplicity, we here write the single-site transition rates $`W(\beta \mathrm{\Delta }E,\beta U)`$ as $`W(\mathrm{\Delta }E)`$.

In the stationary limit, Eqs. (10) and (3) lead to the detailed-balance condition

$$X(T,H)e^{\pm \gamma (\varphi )}\frac{p[\delta \pm 1]}{p[\delta ]}=\frac{𝒲[\delta \delta \pm 1]}{𝒲[\delta \pm 1\delta ]},$$
(12)
where the upper and lower signs have the same interpretation as in Eq. (11). Using Eqs. (3) and (4), we get the stationary values for $`\mathrm{\Pi }_+`$ and $`\mathrm{\Pi }_{}`$,

$`\mathrm{\Pi }_+`$ $`=`$ $`{\displaystyle \frac{Xe^{\gamma (\varphi )}\left(1Xe^{\gamma (\varphi )}\right)}{1X^2}}`$
$`\mathrm{\Pi }_{}`$ $`=`$ $`{\displaystyle \frac{Xe^{\gamma (\varphi )}\left(1Xe^{\gamma (\varphi )}\right)}{1X^2}},`$ (13)
which, when inserted together with Eq. (11) in Eq. (12), yield a self-consistency equation for $`X`$. The self-consistency equation reduces to a linear equation for $`X^2`$, and with the help of the detailed-balance condition for the spin transition rates $`W`$, the solution takes the form

$$X(T,H)=e^{2\beta J_x}\left\{\frac{e^{2\beta H}W(2H4J_x)+e^{2\beta H}W(+2H4J_x)}{W(2H4J_x)+W(+2H4J_x)}\right\}^{1/2},$$
(14)
which is independent of $`\gamma (\varphi )`$. (This solution is the same as the one obtained in Ref. RIKV00 , but the derivation given in that paper did not explicitly show that all dependence on $`\gamma (\varphi )`$ cancels out.)

Equation (14) shows that $`X(T,H)`$ depends on the specific dynamic, except for $`H=0`$, where it reduces to its equilibrium value, $`X(T,0)=e^{2\beta J_x}`$, for all dynamics that satisfy detailed balance. A situation in which the $`H`$-dependence in Eq. (14) cancels out, is that of the soft dynamics discussed in Sec. I. (The barrier energy $`U`$ can be contained in one or the other of the factors.) In Ref. RIKV01  we demonstrated that a soft dynamic yields an SOS interface that is identical to the equilibrium SOS interface at $`H=0`$ and the same temperature, regardless of the value of $`H`$. Neither the Glauber dynamic used here, nor the equally common Metropolis dynamic with transition probability $`W_\mathrm{M}(s_{x,y}s_{x,y})=\mathrm{min}[1,e^{\beta \mathrm{\Delta }E}]`$ KAWA72 ; LAND00 , satisfies this factorization condition. Such hard dynamics lead to a nontrivial field dependence in $`X`$. Inserting the Glauber dynamic defined by Eq. (2) into Eq. (14), we explicitly get

$$X_\mathrm{G}(T,H)=e^{2\beta J_x}\left\{\frac{e^{2\beta J_x}\mathrm{cosh}(2\beta H)+e^{2\beta J_x}}{e^{2\beta J_x}\mathrm{cosh}(2\beta H)+e^{2\beta J_x}}\right\}^{1/2}.$$
(15)
All the results for the spin-class populations of the zero-field equilibrium interface, which are listed in Table 2, can now be applied to obtain a nonlinear-response approximation for the steady-state propagation velocity of flat, driven interfaces with hard dynamics. This simply requires replacing the zero-field $`X=e^{2\beta J_x}`$ used in the linear-response approximation by the field-dependent $`X(T,H)`$, obtained from Eq. (14) using the transition probabilities corresponding to the particular dynamic used. For soft dynamics the linear-response result with $`X=e^{2\beta J_x}`$ is exact RIKV01 .

A physical reason for the marked difference between hard and soft dynamics is best seen by comparing concrete examples of dynamics in the two classes, such as the hard Glauber dynamics used here, Eq. (2), and the soft Glauber dynamics used in Ref. RIKV01 ,

$$W_{\mathrm{SG}}(s_{x,y}s_{x,y})=\left[1+e^{\beta \mathrm{\Delta }E_H}\right]^1\left[1+e^{\beta \mathrm{\Delta }E_J}\right]^1,$$
(16)
in the case of a very strong field. In the hard case, the effect of the field completely dominates the transition rates, such that the rate is near unity for transitions that bring a spin parallel to the applied field, and near zero for transitions in the opposite direction, irrespective of the change in interaction energy. In the soft case, the probability of bringing a spin antiparallel to the field is also near zero, but the probabilities of different transitions bringing a spin parallel to the field differ according to the corresponding change in the interaction energy, as given by the second factor in Eq. (16).

In the next Section we show that this nonlinear-response approximation gives very good agreement with MC simulations of driven, flat SOS interfaces evolving under the hard Glauber dynamic for a wide range of fields and temperatures.

## IV Comparison with Monte Carlo Simulations

We have compared the analytical estimates of step-height distributions, propagation velocities, and spin-class populations developed above with MC simulations of the same model for $`J_x=J_y=J`$. The details of our particular implementation of the discrete-time $`n`$-fold way rejection-free MC algorithm BORT75  are the same as described in Ref. RIKV00 , except that only transitions from the classes with one broken $`y`$-bond ($`k=1`$) are allowed. By keeping only the interface sites in memory, the algorithm is not subject to any size restriction in the $`y`$ direction, and simulations can be carried out for arbitrarily long times.

The numerical results presented here are based on MC simulations at the two temperatures, $`T=0.2T_c`$ and 0.6$`T_c`$ ($`T_c=2J/\mathrm{ln}(\sqrt{2}1)`$ is the critical temperature for the isotropic, square-lattice Ising model ONSA44 ), with $`L_x=\mathrm{10\hspace{0.17em}000}`$ and fixed $`\varphi `$ between 0 and $`45^{}`$. In order to ensure stationarity we ran the simulation for 5 000 $`n`$-fold way updates per updatable spin (UPS) before taking any measurements (100 000 UPS for some of the strongest fields and largest values of $`\varphi `$ at $`T=0.2T_c`$). Exploratory simulations with both larger and smaller $`L_x`$ (up to 100 000) and “warm-up” times (see Sec. IV.1) showed that the values used in the production runs were sufficient to ensure a stationary interface. Stationary class populations and interface velocities were averaged over 50 000 UPS. In the stationary limit 1 UPS corresponds to between 2 MCSS for strong fields at both temperatures, and about 75 MCSS for $`H=0`$ at $`T=0.2T_c`$. Adequate statistics for one- and two-step pdfs was ensured by the large $`L_x`$, ten times the value used in Ref. RIKV00 . Each data point took approximately nine hours on a DEC Alpha Unix workstation or four hours on a Pentium II PC running Linux.

### IV.1 Approach to stationarity

In order both to check the applicability of the mean-field approximation at early times, and to decide the approximate time needed to reach the stationary state, we first studied the transient behavior of the average step height $`|\delta |`$ for $`\varphi =0`$ at $`T=0.2J`$ and 0.6$`J`$.

In the dynamic MC simulation $`|\delta |`$ was measured during a “time window” which was opened after a specified number of UPS, corresponding to a given approximate average evolution time $`t`$ in MCSS. To obtain an optimum balance between time resolution and accuracy, the width of the window was varied from approximately one MCSS at early times, to about ten MCSS at late times, and five independent runs were performed for each value of $`t`$. Standard errors for $`t`$ and $`|\delta |`$ were estimated in the usual way from the spread in their measured values over the five realizations.

For comparison with the MC simulations, we solved the mean-field equation of motion for the single-step pdf, Eq. (10), numerically by a first-order iterative scheme with a time step of 10<sup>-4</sup> MCSS (shorter time steps made no discernible difference). Both the simulations and the solution of the equation of motion were started from a sharp interface at $`t=0`$. Results for fields between 0 and 10$`J`$ are shown in Fig. 2.

In general, we find overall qualitative agreement between the simulations and the equation of motion. For $`T=0.2T_c`$ and $`|H|2J`$, both methods have reached a common stationary value by $`t=\mathrm{10\hspace{0.17em}000}`$ MCSS, while for $`T=0.6T_c`$ and $`|H|3J`$, stationarity is reached by $`t=\mathrm{1\hspace{0.17em}000}`$ MCSS. Our choice of 5 000 UPS as “warm-up time” in our studies of the stationary properties, corresponding to at least $`\mathrm{10\hspace{0.17em}000}`$ MCSS, is thus well justified. However, there are significant quantitative differences between the simulation results and the solution of Eq. (10) for early times. We believe this indicates that the mean-field assumption of statistically independent step heights is not well justified until the interface structure has been “randomized” through a sufficient number of updates. For the extremely strong field, $`H=10J`$, $`𝒲(\delta \delta \pm 1)𝒲(\delta \pm 1\delta )1/2`$, so that the evolution of the interface width is essentially diffusional and $`|\delta |t^{1/2}`$. Both methods agree with this result, although the amplitude for the MC simulation is larger than predicted by the mean-field equation of motion.

### IV.2 Single-step probability densities

Single-step pdfs were obtained by MC simulation at $`T=0.2T_c`$ and $`0.6T_c`$ for $`\varphi =0`$ and several values of $`H`$ between 0 and 3.0$`J`$. The simulation results for $`p[\delta ]`$ are shown in Fig. 3, together with the theoretical result, Eq. (3) with $`X(T,H)`$ from Eq. (15). For both temperatures, the agreement is excellent in the whole range of $`\delta `$ and $`H`$ shown.

A simple comparison between the analytical and simulation results is given in Fig. 4(a), which shows $`|\delta |`$ vs $`H`$ for $`\varphi =0`$ at $`T=0.2T_c`$ and 0.6$`T_c`$. The solid curves represent the theoretical result obtained by summation of Eq. (3), $`|\delta |=2X/\left(1X^2\right)`$, with $`X`$ from Eq. (15). There is excellent agreement between the theoretical field dependence and the MC data. Additional confirmation of the form of the single-step pdf, Eq. (3), is obtained from the simulation results by calculating $`|\delta |`$, both directly by summation over the numerically obtained pdf and from the probability of zero step height as $`|\delta |=\left\{p[0]^1p[0]\right\}/2`$.

A slightly different way to check the agreement between the analytical predictions and the simulation results for the single-step pdf, is to compare $`X(T,H)`$ as given by Eq. (15) with the same quantity obtained from the simulations under the assumption that Eq. (3) holds. From this equation for the pdf, using $`Z(0)`$ from Eq. (6), it follows that $`X`$ is given in terms of $`p[0]`$ as $`X=\{1p[0]\}/\{1+p[0]\}`$, and in terms of $`|\delta |`$ as $`X=\sqrt{1+|\delta |^2}|\delta |^1`$. Both these MC estimates for $`X(T,H)`$ are shown in Fig. 4(b). Again, the agreement is excellent. The slight deviations of the estimate based on $`|\delta |`$ for large $`H`$ are probably due to the fact that data were only recorded in separate bins for $`|\delta |<64`$, so that the calculated average becomes inaccurate whenever higher steps cannot be ignored. The estimate based on $`p[0]`$ does not suffer from this problem. However, the slight discrepancy between both MC estimates on the one hand and Eq. (15) on the other, which is seen between $`H/J=0.5`$ and 2 for $`T=0.2T_c`$, is probably a real effect.

### IV.3 Interface velocities

In this section we compare the simulated interface velocities with the analytical approximation, Eq. (8). Figure 5 shows the normal velocity vs $`H`$ for $`\varphi =0`$. Included are both the linear-response approximation (i.e., $`X=e^{2\beta J_x}`$) and the nonlinear-response result with $`X(T,H)`$ from Eq. (15). Overall, there is excellent agreement between the MC results and the nonlinear-response theory, while the linear-response approximation seriously underestimates the velocity, especially at the lower temperature. As for the Ising model with hard Glauber dynamics studied in Ref. RIKV00 , these results show that the latter approximation is clearly inadequate, and we include no further linear-response results in this paper. A slight disagreement between the simulations and the analytical predictions is seen at 0.2$`T_c`$ in the same range of field values as for $`X(T,H)`$ obtained from the single-step pdf in Sec. IV.2.

The dependence of the normal velocity on the tilt angle $`\varphi `$ is shown in Fig. 6 for several values of $`H/J`$ between 0.1 and 3.0. At $`T=0.2T_c`$ the anisotropy undergoes a gradual change from increasing with $`\varphi `$ in agreement with the polynuclear growth model at small angles and the single-step model for larger angles at weak fields, to Eden-type inverse anisotropy DHAR86 ; MEAK86B ; HIRS86  at strong fields \[Fig. 6(a)\]. At $`T=0.6T_c`$, on the other hand, inverse anisotropy is found for the stronger fields, growing gradually more pronounced with increasing $`H`$ \[Fig. 6(b)\], while for the weakest fields studied the velocity is nonmonotonic in $`\mathrm{tan}\varphi `$ \[Fig. 6(c)\]. The agreement between the simulations and the analytical results is excellent everywhere.

The temperature dependence of the normal interface velocity is shown in Fig. 7 for several values of $`H/J`$ between 0.2 and 3.0. The agreement between the simulations and the analytical results is excellent everywhere.

### IV.4 Spin-class populations and skewness

A closer look at the performance of the mean-field approximation for the interface structure is provided by the mean spin-class populations. The analytical predictions for the class populations are based on the assumption that different steps are statistically independent. A comparison of the simulation results with the analytical predictions therefore gives a way of testing this assumption.

The six mean class populations, $`n(01s)`$, $`n(11s)`$, and $`n(21s)`$ with $`s=\pm 1`$ are shown vs $`H`$ in Fig. 8 for $`\varphi =0`$ and $`T=0.2T_c`$ and $`0.6T_c`$. At both temperatures the analytical approximations follow the average of the populations for $`s=+1`$ and $`s=1`$ quite well, but at intermediate fields in particular, the populations in front of the surface ($`s=1`$) and behind it ($`s=+1`$) are distinctly different.

The skewness between the spin populations on the leading and trailing edges of the interface are a consequence of short-range correlations between neighboring steps, and it is quite commonly observed in driven interfaces. This is the case, even when the long-range correlations vanish as they do for interfaces in the KPZ dynamic universality class, to which the present model belongs. Such skewness was also observed in our study of the Ising model in Ref. RIKV00 , but in that case it was difficult to separate it from the effects of bubbles and overhangs. Skewness has also been observed in other SOS-type models, such as the body-centered SOS model studied by Neergaard and den Nijs NEER97  and a model for the local time horizon in parallel MC simulations studied by Korniss et al. KORN00 . (However, no skewness is observed for the soft Glauber dynamic, a result which may be general for soft dynamics.) The correlations associated with the skewness generally lead to a broadening of protrusions on the leading edge (“hilltops”), while those on the trailing edge (“valley bottoms”) are sharpened NEER97 , or the other way around KORN00 . In terms of spin-class populations, the former corresponds to $`n(21)>n(21+)`$ and $`n(11+)>n(11)`$. The relative skewness can therefore be quantified by the two functions,

$$\rho =\frac{n(21)n(21+)}{n(21)+n(21+)},$$
(17)
introduced in Ref. NEER97 , and

$$ϵ=\frac{n(11+)n(11)}{n(11+)+n(11)}.$$
(18)
These two skewness parameters are shown together in Fig. 9. The relative skewness is seen to be considerably stronger at the lower temperature. This temperature dependence is especially pronounced for $`\rho `$.

Yet another way to visualize the skewness is to consider the joint two-step pdf, $`p[\delta (x),\delta (x+1)]`$. Logarithmic contour plots of this quantity for different values of $`H`$ are shown in Fig. 10 for $`\varphi =0`$ at $`T=0.6T_c`$. It is clearly seen how the contours change with $`H`$. For $`H`$=0 a symmetric diamond shape with equidistant contours indicates statistical independence with single-step pdfs given by Eq. (3). For stronger fields we find shapes that are elongated in the second quadrant \[$`\delta (x)<0`$, $`\delta (x+1)>0`$\] and foreshortened in the fourth quadrant \[$`\delta (x)>0`$, $`\delta (x+1)<0`$\]. This shape indicates that large negative $`\delta (x)`$ tend to be followed by large positive $`\delta (x+1)`$ (sharp valleys), while positive $`\delta (x)`$ tend to be followed by smaller negative $`\delta (x+1)`$ (rounded hilltops). We have not been successful in attempts to construct an analytical approximation which describes this evolution of the joint two-step pdf with increasing field.

## V Conclusion

In this paper we have considered in detail the microstructure of an unrestricted solid-on-solid (SOS) interface with Glauber dynamics, which is driven far from equilibrium by an applied field. The microstructure is of interest because it determines a number of interface properties, such as mobility and chemical reactivity. We adapted to this model a mean-field, nonlinear-response approximation previously developed for driven Ising interfaces without the SOS restriction RIKV00 . In comparison to the Ising driven interface, which leaves bubbles of the unstable phase in its wake and exhibits “overhangs,” the SOS interface is a relatively simple object. The absence of overhangs and of fluctuations in the stable and unstable phases (bubbles behind and in front of the interface) makes the SOS interface more suitable for description in terms of a mean-field type model. Moreover, unlike the Ising model, in which there are several effects that simultaneously contribute to the inaccuracy of the approximate treatment, the simpler SOS structure makes it possible to identify the short-range correlations as the only significant factor causing deviations between the true interface behavior and the mean-field theory.

To study the microstructure of the interface in detail, we investigated the interface velocity as a function of driving field, temperature, and angle relative to the lattice axes. We also studied the local shape of the interface in terms of the spin-class populations and the probability density for individual steps in the interface. In essentially all cases we found excellent agreement between our theoretical description of the stationary moving interface and the results of our dynamic MC simulations.

The microstructure of the moving interface depends crucially on the details of the stochastic dynamics, and for the Glauber dynamic used here, the average height of a step in the interface was found to increase strongly with the applied field. Our theory predicts that this should be the case (with quantitative variations depending on the particular dynamic) for any dynamic in which the single-spin transition rates cannot be factorized into a part that depends only on interaction energies and another that depends only on the applied field (hard dynamics MARR99 ). In contrast, for factorizing (soft) dynamics the interface structure should remain independent of the field. This was recently confirmed for the unrestricted SOS interface by MC simulations RIKV01 . It is therefore important to use great caution in drawing conclusions about the microstructure of driven interfaces, based on dynamic MC simulations. For such conclusions to be valid, the dynamic must be chosen appropriately to the physical system of interest. The hard type of dynamics would appear to be particularly suited for certain interfaces in magnetic or dielectric systems, where the local order parameter is not conserved.

One aspect of the interface dynamics not completely captured by our model is, of course, the short-range correlations. Namely, within the mean-field approximation used here, individual steps of the interface are assumed to be statistically independent. However, for increasing fields the interface undergoes a gradual breakdown of up/down symmetry. This is clearly seen in our simulations here, as well as in other examples of driven interfaces NEER97 ; KORN00 . It would seem likely that one could construct a mean-field approximation at the two-spin level, which might be able to predict this skewness for hard dynamics, as well as its absence for soft dynamics. However, such a theory has not yet been developed.

## Acknowledgments

We acknowledge useful conversations with B. Schmittmann and K. Kawasaki. P. A. R. appreciates the hospitality of the Faculty of Integrated Human Studies at Kyoto University. The research was supported in part by National Science Foundation Grant Nos. DMR-9981815 and DMR-0120310, and by Florida State University through the Center for Materials Research and Technology and the School of Computational Science and Information Technology.