Morphology transition at depinning in a solvable model of interface growth in a random medium

The motion of a driven interface in a random medium keeps attracting wide attention in condensed-matter and statistical physics, as it occurs in many different physical processes such as fluid invasion in porous media, domain wall motion in magnetic systems, contact line motion in wetting, etc. [1,2]. In such systems, the driven interface displays a transition (at zero temperature) from a pinned phase to a moving phase as the driving force surpasses some finite threshold value. This depinning transition has been extensively studied numerically and theoretically, in particular for continuum models of elastic interfaces in the presence of quenched disorder [3,4].

Driven interfaces also display remarkable changes in the domain growth morphology as a function of the disorder strength, an issue that has been explored in detail within the random-field Ising model (RFIM) for various dimensions, coordination numbers, and distributions of the random fields [5–10]. Specifically, in three dimensions, for an unbounded distribution of the random fields (e.g. Gaussian), one observes a transition from a compact-growth regime with a self-affine interface at low disorder to a self-similar, percolation-like growth regime at high disorder [9]. These two regimes are separated by a multicritical point and described by different sets of critical exponents.

The RFIM on a three-dimensional lattice, however, can only be studied numerically, and it would be useful to have a simple and exactly solvable model of interface growth that exhibits, at least qualitatively, the same type of morphology changes. This goal is achieved in the present letter where we introduce an interface-growth dynamics for the RFIM on a Cayley tree for which evolution equations can be written exactly, and solved, in the thermodynamic limit¹ . We show that the phase diagram of the model (critical depinning field vs. disorder strength) is qualitatively similar to that obtained on the cubic lattice [9]. As there is no clear-cut definition of the height and width of the interface on a Cayley tree, we then introduce a specifically tailored random graph (randomly connected chains) with a local tree structure that still makes the problem analytically tractable. The phase diagram is the same as for the Cayley tree and is shown to be indeed associated with the two distinct growth morphologies observed on the cubic lattice [9].

We consider the RFIM on a Cayley tree of degree (coordination number) c + 1. The Hamiltonian is given by

$$\begin{equation} {\mathcal H}(\boldsymbol{\sigma}) = -\frac{1}{2} \sum_{i}\sum_{j\in \partial i} \sigma_i\sigma_j -\sum_i (H + h_i)\sigma_i,\end{equation} \tag{ 1 }$$

where i runs over all sites (nodes) of the graph, $\boldsymbol{\sigma}\equiv\{\sigma_i\}$ with $\sigma_i=\pm 1$, $\partial i$ is the set of sites that are directly connected to i, H is the external field, and the random fields $\{h_i\}$ are independently taken from a Gaussian distribution $\rho(h)$ with zero mean and variance R². The strength of the exchange interaction is taken as the energy unit.

The Cayley tree can be visualized as a root (origin) from which different generations emanate. At the n-th generation level then, there are cⁿ sites. To describe interfaces, however, it is more convenient to start from an outer boundary at a very large generation level L and to proceed inward towards the origin: we therefore label a site i by two indices, p and z, with z denoting the number of generations counted from the outer boundary (z = 0 at the boundary and z = L at the origin) and p labeling the $c^{L-z}$ sites at the level z.

We study an interface-growth dynamics that is analogous to the zero-temperature front-propagation dynamics considered in [5–10]. The interface is defined as the collection of spins that are −1 ("down") and have at least one direct neighbor on the graph that is  +1 ("up"). (Note that this includes spins on the foremost "front" as well as at the interfaces and isolated down spins left behind.) The initial condition at time t = 0 corresponds to all spins at the boundary z = 0 in the up state and all the other spins in the down state. The configuration of spins is then evolved at a fixed external field H with the following rule: only spins at the interface are allowed to flip and a spin at site i on the interface at time t − 1 flips to +1 at time t when the local field $H+h_i +\sum_{j\in \partial i}\sigma_j(t-1)$ is positive (which then lowers the total energy of the system); all unstable spins at the interface are flipped simultaneously at each time step (parallel dynamics). It is worth stressing that this interface growth dynamics for the RFIM at zero temperature is different from the dynamics used to describe the hysteresis loop of the magnetization in the same model [14].

For a sufficiently large negative value of H, the interface remains pinned whereas for a sufficiently large positive value, it keeps moving. In between a depinning transition takes place at a critical value H_c that depends on the disorder strength R. This phenomenon can be captured by studying the probability P_t that a randomly chosen spin at time t and at level z = t is up: in the large t limit, $P_t=0$ in the pinned phase and P_t > 0 in the moving one. (By construction, the interface cannot be higher than level t at time t.) Thanks to the tree structure of the graph, one can easily write down a recursion equation for P_t:

$$\begin{equation} P_t=\sum_{k=1}^{c} {\binom ck} p_k (H) P_{t-1}^k[1-P_{t-1}]^{c-k},\end{equation} \tag{ 2 }$$

where the k-th term of the above sum is the probability that the chosen site at level t, which has a down neighbor at level t + 1, has $k\geq 1$ up and c − k down neighbors at level t − 1 and has a positive local field; $p_k(H)$ is then the probability that the local field is positive when k neighbors are up and $c+1-k$ are down, i.e. $p_k(H)= \int_{- H -2k + (c+1)}^{\infty} \text{d}h \rho(h)$.

For large t, P_t tends to a fixed point $P^*$ given by the self-consistent equation obtained from eq. (2). This is a polynomial equation that has generically c solutions, including $P^*=0$. For concreteness we discuss more specifically the case c = 3, but similar results are obtained for c > 3. The two nontrivial solutions are then given by $P_{\pm}^*(H,R)=-(v \pm \sqrt{v^2-4uw})/2u$, where $u\equiv (p_3-3p_2+3p_1)$, $v\equiv 3(p_2-2p_1)$, and $w\equiv (3p_1-1)$.

As one increases H from a large negative value, the only physical fixed point is at first $P^*=0$ (the two other fixed points are initially real but unphysical, and then become complex). Below some critical disorder $R_{\text{c}}$, the two complex solutions merge ($v^2=4uw$) at a critical field $H=H_{\text{sn}}(R)$ and become real: $P_+^*=P_-^*=-v/(2u)$, with $0<-v/(2u)\leq 1$. $P_+^*$ is then the stable fixed point for $H\ge H_{\rm sn}(R)$ while $P_-^*$ remains unstable. This corresponds to a saddle-node bifurcation. On the other hand, for a strong disorder $R>R_{\text{c}}$, $P_+^*$ becomes the physical, stable fixed point (while $P_-^*$ is unphysical) by passing through zero from the negative side $(w=0)$ at a critical field $H=H_{\text{tc}}( R)$. This is a transcritical bifurcation. The two distinct bifurcations meet at a multicritical point $(R_{\text{c}},H_{\text{c}}(R_{\text{c}}))$ where $P_+^*=P_-^*=0$ (implying $v=w=0$, i.e., $p_2=2p_1=2/3$): this yields $R_{\text{c}}= 2.32165\cdots$ and $H_{\text{c}}(R_{\text{c}})=1$.

One can also describe the long-time dynamical behavior of P_t in the thermodynamic limit L → ∞. Below the line formed by $H_{\text{sn}}(R)$ and $H_{\text{tc}}(R)$, P_t approaches zero exponentially fast as $t\rightarrow \infty$ and one is in the pinned phase. Above the line, P_t approaches the nonzero fixed-point value $P_+^*(H,R)$ exponentially fast and one is in the moving phase. The line corresponds to a depinning transition whose nature is different above and below $R_{\text{c}}$:

• i)  
  Above $R_{\text{c}}$, the transition is continuous and P_t approaches zero with a power law $t^{-1}$. When H is close to $H_{\text{tc}}(R)$, P_t can be cast in a scaling form

  $$\begin{equation} P_t\simeq \vert H- H_{\text{tc}}\vert f_{\pm}(\vert H- H_{\text{tc}}\vert t)\end{equation} \tag{ 3 }$$

  with $f_-(-\infty)=0$, $f_+(\infty)>0$, and $f_{\pm}(x\rightarrow 0)\sim x^{-1}$. (More precisely, $P_t\simeq (-1/v)t^{-1}$ with v < 0 for $H=H_{\rm tc}$.) The transition then belongs to the same universality class as the mean-field percolation.
• ii)  
  Below $R_{\text{c}}$, the transition takes place with a jump in the asymptotic value of P_t. On the moving side, near the transition, there is a slow approach to the asymptotic value $P_+^*(H_{\text{sn}},R)=-v/(2u)$ and P_t can be cast in a scaling form

  $$\begin{equation} P_t-P_+^*\simeq (H- H_{\text{sn}})^{1/2}f_{+}( (H- H_{\text{sn}})^{1/2} t)\end{equation} \tag{ 4 }$$

  with $f_+(\infty)>0$ and $f_{+}(x\rightarrow 0)\sim x^{-1}$. (More precisely, $P_t-P_+^*\simeq (2/v)t^{-1}$ with v > 0 for $H=H_{\text{sn}}$ (see footnote ² ).) On the pinned side, close to the transition, the trajectory of P_t is first attracted to the fixed point $P_+^*(H_{\text{sn}},R)$ as t increases, but eventually, at a rather well-defined crossover time $t_{\text{co}}$, P_t goes exponentially to the stable fixed point 0. The crossover time $t_{\text{co}}$ diverges as $(H_{\text{sn}}-H)^{-1/2}$ as $H \nearrow H_{\text{sn}}$. Thus, the transition at the saddle-node bifurcation has a mixed, first- and second-order, character with a jump in the order parameter and a diverging characteristic time scale. A similar situation is found in the zero-temperature Glauber dynamics of the random-field Ising model [15], the for k-core (or bootstrap) percolation [16–18] and in kinetically constrained models [19].
• iii)  
  Finally, near the critical point $R_{\text{c}}$ separating the two types of transitions, one finds that along the saddle-node transition line, $P^{*}_+(H_{\text{sn}},R) \sim (R_{\text{c}}-R)$ when $R \nearrow R_{\text{c}}$. Exactly at R_c, P_t can also be put in a scaling form

  $$\begin{equation} P_t\simeq \vert H-H_{\text{c}}\vert ^{1/3}f_{\pm}(\vert H-H_{\text{c}}\vert^{2/3}t)\end{equation} \tag{ 5 }$$

  with $f_-(-\infty)=0$, $f_+(\infty)>0$, $f_{\pm}(x\rightarrow 0)\sim x^{-1/2}$. (More precisely, $P_t\simeq(t^{-1/2}+t^{-3/2}/4)/\sqrt{-2u}$ with u < 0 at $(R_{\text{c}},H_{\text{c}})$.)

The phase diagram corresponding to the interface growth model for the RFIM on a Cayley tree with c = 3 is shown in fig. 1. There is a striking resemblance with the phase diagram obtained for the same model on the cubic lattice: see fig. 3 of [9]. In both cases, the nature of the depinning transition changes at a critical value of the disorder close to the maximum of the transition line $H_{\text{c}}(R)$. On a cubic lattice the growth proceeds with a self-affine interface at low disorder and with a self-similar one at high disorder. On the Cayley tree, a similar pattern is found with a continuous, percolation-like, transition at high disorder and a mixed transition with features akin to those of the self-affine case at low disorder (see also below). In addition, just as observed on Euclidean lattices, only the continuous percolation-like transition survives when the coordination number is small enough (here, c = 2, which can be compared to the case of the honeycomb lattice studied in [6] for bounded distributions of random fields). Furthermore, all transitions disappear when c = 1, which corresponds to d = 1.

Zoom In Zoom Out Reset image size

The distinct signature of the changes of growth mechanism at the depinning transition on a cubic lattice lies in the morphology of the interface. A central quantity to characterize the latter is the height of the "front" (the foremost part of the interface) at time t, $l_p(t)$, for each position $p=1,\ldots,c^L$ on the initial boundary; l_p can be defined as the maximum z at which one finds an up spin in the shortest path connecting the tree origin $(z=L)$ to the position p at the boundary $(z=0)$ (see footnote ³ ). However, there are some difficulties to consider this quantity on the Cayley tree, due to the peculiar, hyperbolic-like, geometry of the latter, with an exponentially decreasing number of sites with increasing height z. It is hard to visualize and assess the interface morphology in such a case, and numerical experiments are extremely demanding.

To sidestep this issue and provide some insight on possible changes in the interface morphology while retaining the analytically solvable character of tree-like graphs, we introduce a specifically tailored random graph, which can be described as randomly connected chains (see fig. 2). It consists of L_x parallel vertical chains of connected sites of height $L \geq L_x$. A site at height z is (deterministically) connected to its two nearest neighbors on the same chain, as well as randomly connected to b additional sites at level z − 1 and to b additional sites at level z + 1 on the other chains. The total coordination number of each site is therefore $2(b+1)$, except for the bottom sites forming the initial interface (at z = 0) and the top sites at z = L, for which it is b + 1.

Zoom In Zoom Out Reset image size

In spite of the presence of chains, the graph has a local tree structure. From heuristic arguments one expects that the typical loop length essentially behaves as $\log L_x$ when $L_{x}\rightarrow \infty$, so that loops become irrelevant in the thermodynamic limit. A rough argument goes indeed along the usual lines [20,21] as follows. Consider a random site, say at height z, and assume that in a ball of large radius ℓ the graph is a tree: the number of sites in the ball then goes essentially as $[2(b+1)]^{\ell}$. On the other hand, the total number of sites in the strip of the graph including the ball is $2\ell\, L_x$. The tree-like assumption then breaks down and loops must appear when $[2(b+1)]^{\ell}\sim 2\ell\, L_x$, i.e. when $\ell \sim \log L_x$ in the large L_x limit. A somewhat refined argument accounting more precisely for the presence of (determistic) chains of links is given in the online supplementary material [22], but it leads to the same conclusion, up to subdominant terms. (Although the height L does not explicitly appear in the reasoning, it should of course be large, which is satisfied by requiring that L is of the order of L_x.) Note finally that the tree-like character of the graph in the limit $L_x \rightarrow \infty$ is also confirmed by the excellent agreement between numerical simulations and theoretical expressions (see below and online supplementary information [22]).

The main advantage of the randomly connected chains is that there is now a natural definition of the front height l_p by considering the situation on the pth chain with $p=1,\ldots, L_{x}$: l_p corresponds to the highest up spin on the chain. The interface growth proceeds in a similar way as on the Cayley tree (see above) and one follows the motion of the interface from the initial boundary at z = 0. However, to reproduce the results previously obtained on the Cayley tree, it turns out to be necessary to introduce a modification in the dynamical rule. Otherwise, the peculiar structure of the randomly connected chains allows for additional back-reaction effect compared to the Cayley tree: a spin on the front that flipped at time t can influence the front on another chain at a later time by triggering the flip of spins at lower heights. To avoid this effect, spins at height z < t are not allowed to flip at time t. This modification is the price to pay for mimicking real-space interfaces while retaining a tree-like structure. (Note that this does not suppress the fact that there may be several interfacial spins on the same chain. On the other hand, it may change the roughness of the front.)

With this additional constraint, the probability P_t that a randomly chosen spin flips up at z = t can still be described by an exact recursion equation (in the thermodynamic limit). At time t all spins at $z=t+1$ are down by construction, and one only needs to consider the various possibilities for the b + 1 spins at level $z=t-1$. One can then use the statistical independence of the branches reaching the chosen site at z = t to derive

$$\begin{eqnarray} &&P_t = \sum_{k=1}^{b+1} {\binom {b+1}{k}} \tilde p_k (H)\, P_{t-1}^k[1-P_{t-1}]^{b+1-k}, \end{eqnarray} \tag{ 6 }$$

where $\tilde{p}_k(H)= \int_{- H -2k + 2(b+1)}^{\infty} \text{d}h \rho(h)$. By comparison with eq. (2), one concludes that the interface-growth dynamics on the randomly connected chains exactly reproduces that on the Cayley tree provided one chooses $b=c-1$ and shifts the external field from H to $H+c-2b-1=H-b$. For b = 2, the phase diagram of the model is thus identical to the one shown in fig. 1, with the y-axis shifted by 2. To illustrate the change of growth morphology as a function of disorder, we show in figs. 3(a) and (b) the results of numerical simulations performed near the depinning transition in the two different regimes. It is manifest that the interface roughness is much larger in the percolation-like regime at high disorder.

Zoom In Zoom Out Reset image size

The interest of the model is that it remains analytically tractable. To provide further information we consider the probability $Q_t(l)$ that the front height $l_p(t)$ on a randomly chosen chain p is equal to l. (Note that due to the modification of the dynamics, the only meaningful interface is now the foremost front as spins behind are no longer allowed to evolve.) A closed-form equation for $Q_t(l)$ can be derived in the thermodynamic limit thanks to the local tree structure of the randomly connected chains. It is useful to first introduce the conditional probability $q^{\pm}_t$ that a spin at a site $\{p,z=t\}$ flips at time t provided that the spin just below at $\{p,z=t-1\}$ is equal to $\pm 1$; by using the statistical independence of the branches arriving at the site under consideration, it can be expressed as

$$\begin{equation} q^{\pm}_t=\sum_{k=1-\delta_{\sigma,1}}^{b} {\binom bk} \tilde{p}_{k+\delta_{\sigma,1}}(H) P_{t-1}^k[1-P_{t-1}]^{b-k},\end{equation} \tag{ 7 }$$

where $\delta_{\sigma,1}$ is the Kronecker symbol and $\sigma=\pm 1$. Recalling that at time t spins are allowed to flip only at height z = t and not below and that a spin at the front is the highest up spin along its chain, one immediately finds

$$\begin{equation} Q_t(l=t)=P_t.\end{equation} \tag{ 8 }$$

The probability that the front height is at $l=t-1$ at time t requires that a spin flipped at $z=t-1$ at time t − 1 and that the spin directly above it does not flip at time t, hence,

$$\begin{equation} Q_{t}(l=t-1)=P_{t-1}[1-q^+_t].\end{equation} \tag{ 9 }$$

Finally, along the same lines, one easily derives that for lower heights $0\leq l \leq t-2$,

$$\begin{equation} Q_t(l)=Q_{t-1}(l)[1-q^-_t].\end{equation} \tag{ 10 }$$

From $Q_t(l)$ one can compute the moments $\langle l^n\rangle_t\equiv \sum_{l=0}^t l^n Q_t(l)$ and derive the mean height $l_0(t)\equiv\langle l\rangle_t$ and the mean width $W_0(t)\equiv[\langle l^2\rangle_t-l_0(t)^2]^{1/2}$. To study the asymptotic behavior at large time of l₀ and W₀, it is convenient to use the above equations to derive

$$\begin{equation} l_{0}(t)=tP_{t}+(t-1)P_{t-1}[1-q^+_{t}]+X_{t}^{(1)},\end{equation} \tag{ 11 }$$

$$\begin{equation} W_{0}(t)^2+l_0(t)^2=t^2P_{t}+(t-1)^2P_{t-1}[1-q^+_{t}]+X_{t}^{(2)},\end{equation} \tag{ 12 }$$

where $X_t^{(n)}$ for $n=1,2$ satisfies the following equation:

$$\begin{align} X_{t}^{(n)}=[1-q^-_{t}]\left(X_{t-1}^{(n)}+(t-2)^nP_{t-2}[1-q^+_{t-1}]\right). \end{align} \tag{ 13 }$$

The asymptotic analysis for $t\rightarrow \infty$ then leads to the following results:

• i)  
  In the pinned phase, both $l_0(t)$ and $W_0(t)$ are of $\mathcal{O}(1)$, while in the moving phase, $l_0(t)\simeq t$ and $W_0(t)=\mathcal{O}(1)$.
• ii)  
  At the discontinuous (saddle-node) transition for ${R<R_{\text{c}}}$, $l_0(t)\simeq t$ and $W_0(t)=\mathcal{O}(1)$.
• iii)  
  At the continuous percolation-like transition (transcritical bifurcation) for $R>R_{\text{c}}$, $l_0(t)\simeq a t$ with $a=2/(8-9\tilde{p}_2)<1$ and $W_0(t )\simeq b t$ with $b=3(2 -3 \tilde{p}_2)/[(8-9\tilde{p}_2)(7-9\tilde{p}_2)^{1/2}]<1$.
• iv)  
  Finally, at the multicritical point for $(R_{\text{c}},H_{\text{c}})$, $l_0(t)\simeq t$ and $W_0(t)\simeq c t^{1/2}$ with $c=3/\sqrt{-2u}\simeq 0.664$. The full solution is illustrated in fig. 4.

Thus, as anticipated by looking at fig. 3, the behavior of $l_0(t)$ and $W_0(t)$ provides a clear signature of the change in the morphology of the front at depinning as a function of disorder, from rather compact for $R<R_{\text{c}}$ to rough and percolation-like for $R>R_{\text{c}}$.

Zoom In Zoom Out Reset image size

In addition to the analytical study we have also carried out numerical simulations of the interface-growth process. It allows us, on the one hand, to confirm that the randomly connected chains is equivalent to a tree in the thermodynamic limit and, on the other hand, to assess finite-size characteristics of the interface (front) that are hard to compute analytically. We focus on the mean width at time t for graphs of lateral extension $L_{x} \gg 1$ (typically from 50 to 6400):

$$\begin{equation} W(t,L_{x})\equiv\overline{[\langle l_p(t)^2\rangle - \langle l_p(t)\rangle^2]}^{1/2},\end{equation} \tag{ 14 }$$

where the overline denotes an average over the random fields and the random graphs (the number of samples varies from 4096 for $L_x=50$ to 32 for $L_x=6400$), and $\langle \; \rangle$ an average over the L_x chains, $(1/L_{x})\sum_{p=1}^{L_{x}}$. When $L_{x}\rightarrow \infty$, $W(t,L_{x})$ goes to $W_0(t)$ already computed (see also the supplementary information [22]). As can be seen in fig. 5, $W(t=L_{x},L_{x})$ at the depinning transition appears to follow a power law $(L_{x})^{\alpha}$ at large L_x, where α can then be taken (with a pinch of salt due to the nature of the graph) as a "roughness" exponent. The numerical data leads to $\alpha=1/2$ for $R>R_{\text{c}}$ and $\alpha=0$ for $R<R_{\text{c}}$. (The data for $R=R_{\text{c}}$, not shown here, seem compatible with $\alpha \simeq1/2$ but the convergence to the asymptotic limit is much slower so that the estimate is not as reliable as for the other cases.) From all the above results, the mean interface width at the depinning transition is then expected to follow a scaling form

$$\begin{equation} W(t,L_{x})\simeq (L_{x})^{\alpha} w[t/(L_{x})^{\alpha/\beta}]\end{equation} \tag{ 15 }$$

with $w[\infty]=\mathcal{O}(1)$, $w[y\rightarrow 0]\sim y^\beta$, $\beta=0$ for the saddle-node bifurcation (on the moving side), $\beta=1$ for the transcritical bifurcation, and $\beta=1/2$ for the multicritical point.

Zoom In Zoom Out Reset image size

In summary, we have proposed a solvable model of interface growth in a random medium on a Cayley tree that displays a morphology transition at depinning as a function of disorder strength. Along the depinning line, a multicritical point separates a low-disorder regime where the transition has a mixed continuous-discontinuous character and is described by a saddle-node bifurcation and a high-disorder regime where the transition is continuous and characterized by a transcritical bifurcation.

To analyze the interface (front) morphology in more detail, we have introduced a trick: we have replaced the problem on a Cayley tree by one on a graph formed by randomly connected chains. With some adjustments in the dynamics, the latter is described by the same recursion equations as the former, but it allows a clear definition and visualization of the front. The results then illustrate that the low-disorder regime corresponds to a compact growth while the high-disorder one corresponds to percolation-like growth, showing that the multicritical point can be associated with a morphology transition. Remarkably, the phase diagram is similar to that of the same model on the cubic lattice [9]. In the latter case the multicritical point was conjectured to coincide with the maximum observed in the depinning transition line $H_{\text{c}}(R)$ (error bars in the simulations make it difficult to settle this issue) [9]. We show here analytically that on tree-like graphs this maximum unambiguously takes place in the low-disorder region (see fig. 1), below the multicritical disorder strength.

We are grateful to A. Rosso, T. Sasamoto and G. Semerjian for useful discussions and HO thanks the Yukawa Institute for Theoretical Physics where this work was initiated. This work is also supported by EC-Grant "STAMINA" No. 265496.