Influence of Disorder Strength on Phase Field Models of Interfacial Growth

We study the influence of disorder strength on the interface roughening process in a phase-field model with locally conserved dynamics. We consider two cases where the mobility coefficient multiplying the locally conserved current is either constant throughout the system (the two-sided model) or becomes zero in the phase into which the interface advances (one-sided model). In the limit of weak disorder, both models are completely equivalent and can reproduce the physical process of a fluid diffusively invading a porous media, where super-rough scaling of the interface fluctuations occurs. On the other hand, increasing disorder causes the scaling properties to change to intrinsic anomalous scaling. In the limit of strong disorder this behavior prevails for the one-sided model, whereas for the two-sided case, nucleation of domains in front of the invading front are observed.


I. INTRODUCTION
Interface growth in disordered systems has been a subject of great interest in non-equilibrium statistical physics for more than a decade. Many different examples of interfaces growing in heterogeneous systems have been found in phase separation by nucleation and growth [1,2], solidification [3,4], and fluid invasion in porous media [5,6], among others.
Phase field models of increasing complexity have in recent years been extensively used in studying interface roughening as well as microstructure formation [5,7,8,9,10,11,12,13]. A particularly interesting problem of interface roughening is that associated with a fluid invasion front moving into a disorder medium [5], which can be experimentally studied with the Hele-Shaw cell set-up [14,15,16,17,18]. In modeling such a fluid invasion experiment, two different ways to consider the mobility parameter in a Model B type of phase field model, called the one-sided and symmetric models, were used by Hernández-Machado et al. [19] and Dubé et al. [20,21], respectively. The difference between these two cases is that in the two-sided model the mobility is constant throughout the system, while in the one-sided model it becomes zero in the phase which is being invaded.
In this paper our aim is to carry out a detailed analysis about the influence of the strength of the disorder in the two types of models described above. Both cases can be described by a generalized Cahn-Hilliard equation or Model B with quenched disorder in the background medium. The boundary conditions are used to couple the system to a reservoir of the invading phase with a constant mass flux. At the linear level of small front fluctuations, both phase field models can be analyzed through linearized interface equations, that is, the evolution of the interface can be described in terms of the interface profile alone. Moreover, the bulk diffusion fields implicit in this description cause the interface equation to become spatially non-local [19,20,21]. It is thus of interest to examine how the models are influenced by varying disorder strength, which can be easily realized in the experiments, too. To this end, we first define a critical value of the disorder strength σ c , above which the disorder becomes strong. We find that at weak disorder strengths (σ < σ c ) both models have the same interface scaling behavior given by super-rough anomalous scaling, which is also predicted by the linear theory. We observe clear deviation from this scaling when the disorder strength comes close to σ c , where super-roughness disappears and weak intrinsic anomalous scaling arises. Furthermore, in the limit of strong disorder (σ > σ c ), the two models are no longer equivalent. The one-sided model can still be applied to describe diffusive liquid invasion in good agreement with experimental results of Ref. [16]. However, in this limit the symmetric model exhibits nucleation of the invading phase in front of the advancing interface.
The structure of the paper is as follows. In Section II we introduce the two versions of the phase field model, consider the linearized interface equation (LIE) valid in the weak disorder limit, and give an estimate for the strong disorder limit, thus introducing a scale for the disorder strength. Section III defines the concepts of scaling in interface roughening, including super-rough and intrinsic anomalous scaling. In Setion IV we present our numerical results, and finally give our conclusions in Section V.

II. THE PHASE FIELD MODEL
The phase field model we are using describes a system of two phases separated by an interface. We introduce a locally conserved field φ to represent the two phases of the problem with the equilibrium values φ e = ±1 in such a way that the interface position is at φ(x, h) = 0. The dynamics of the field is then assumed to follow a conserved equation based on a time-dependent Ginzburg-Landau Hamiltonian (model B [7]): where the current is proportional to the gradient of the chemical potential j = −M (φ)∇µ. Here, the chemical potential is µ = δF /δφ and the free energy is taken to be of the form F [φ] = dr r r(V (φ) + (ǫ∇φ) 2 /2). A simple double-well potential is chosen with a linear random term: The variable η(r) is taken to be stochastic and it plays the role of spatially quenched disorder in the system. The disorder is characterized by its average η , and its standard deviation σ, which characterizes the disorder strength. The disorder also has a correlation length l corr , which in a numerical scheme is most conveniently set to the lattice spacing. Note that by considering η as local average over area of size l corr , the choice of the length l corr also enters the disorder strength. That is, the standard deviation of η is σ when observed at scale l corr . The surface tension of this model can be calculated with the disorder-free "kink" solution φ 0 , or Goldstone mode, given by the lowest energy solution of the boundary conditions φ(−∞) = −1 and φ(+∞) = 1. The resulting dimensionless surface tension is γ = √ 2/3 ≃ 0.47 [9]. The equation for the dynamics of the phase field reads where M (φ) is a mobility parameter which may depend on the phase field. In the sharp interface limit ǫ → 0, the normal velocity of the interface can be obtained by integrating equation (3) in a region around the interface, where the subscripts + and − correspond to the two phases of the system. In our study, two different functional forms for the mobility M (φ) will be considered. First, we assume a symmetric parameter M = M 0 constant for the whole system independent of the field φ.
In this case, the velocity of the interface is controlled by the difference between incoming and outgoing currents. Second, we will consider a one-sided parameter M = M 0 θ(φ), θ(φ) being the Heaviside step function, which is zero in the phase where φ < 0. Note that the normal velocity of the interface then reduces to v n ≃ −M 0 ∂ n µ| − , that is, it is only proportional to the outgoing current. As we will see in our numerical results, the two models can give different results and describe different physical situations depending on the strength of the disorder.
Here we consider the so-called driven case, where the mean velocity of the interface is fixed to a constant value. The relevant boundary condition is to impose a fixed gradient of the chemical potential at the bottom of the system, ∇µ = −Vŷ. Combined with the initial condition of the boundary condition leads to phase φ = +1 invading phase φ = −1 and the interface moving with a constant average velocity. In the general context of the phase field model we will refer to these phases as A and B, respectively. In terms of liquid front invasion into a Hele-Shaw cell, these phases would be liquid and air, respectively, whereas in terms of phase separation they would be the phases rich in components A and B.

A. The linearized interface equation
The dynamics of a front in the phase field model described above can also be considered in terms of an integro-differential equation for the interface, which reduces the dimensionality of the problem by one, obviously a desirable property. However, this description can only be obtained as a perturbation expansion around a flat front with small fluctuations, and since the fluctuations are caused by the disorder, this also means weak disorder. Because the front dynamics are only obtained to first order in small fluctuations we call it the linearized interface equation (LIE). It is noteworthy that the longranged effects of mass conservation in the phase field results in the LIE being spatially non-local, even though it is linear. This means that the LIE takes the form of convolutions in position, which can be made local in Fourier space.
The LIE can be obtained using the Green's function G to express the phase field equation, Eq. (3), in terms of an integro-differential equation [20], and then projecting it to an interfacial description with the operator . The LIE takes the form of dispersion relation for Fourier components of small front fluctuationsĥ(k, t) around the average interface height where γ is surface tension,Ḣ 0 = V is the (constant) average front propagation velocity, and the disorder term is the Fourier transform of the two-dimensional (2D) disorder along the front, orη(k, t) = dxe −ikx η(x, H(x, t)).
From the dispersion relation (6), one immediately obtains the crossover length scale when the two dispersion terms are equally strong. Physically the dominant dispersion mechanism then changes from surface tension to mass transport. In addition, because the conservative character of the capillary disorder term |k|η(k, t), it turns out that the crossover length ξ × acts as an upper cutoff for the correlation length of fluctuations. This means that in the long wavelength region, where mass transport controls the dissipation of front fluctuations, the interface is asymptotically flat. This effect has been numerically shown in Refs. [5,20,21,22], and also in a general context of kinetic roughening [23].

B. Definition of disorder strength
In order to study the influence of disorder strength in the phase-field model of Eq. (3), we need to define a measure for the relative strength of the disorder. This can be achieved by comparing the disorder contribution in the dimensionless bulk free energy of Eq (2) to the surface energy. We do this by considering a domain of linear size r, where the local disorder average η r is a stochastic variable with standard deviation σ r ∝ r −1 σ, where σ is the standard deviation of a single disorder site, which is of linear size l corr . The underlying disorder then has a correlation length l corr . Considering a circular domain of radius r, it is energetically favorable for this domain to be of the opposite phase than its surroundings if where the miscibility gap in our model is ∆φ = 2. The LHS is the energy cost of the interface, whereas the RHS is the energy gain due to disorder. We consider the local disorder average η r on the fluctuation site to be as large as its standard deviation σ r = σ l corr / √ πr. Then Eq. (8) gives the condition for when the disorder can locally dominate the bulk energy and thus is defined to be strong. The order of magnitude estimate for strong disorder in our dimensionless units (l corr = 1, ǫ = 1) is thus obtained as the variance being of the same order as the surface tension σ = σ c ≈ γ. Note in particular that no r dependence remains in the estimate, and thus the relative disorder strength will be the same at all length scales (larger than the interface width ǫ and disorder site size l corr ).

III. SCALING OF ROUGH INTERFACES
The statistical treatment of a 1D interface H(x, t) is usually done by studying the scaling properties of its fluctuations over the whole system of total size L [24]. For scale-invariant growth the lateral correlation length of the fluctuations is expected to grow in time following a power law ℓ c ∼ t 1/z until it reaches the system size L, defining a saturation time t s ∼ L z . Alternatively, the global width of the interface W (L, t) = [H(x, t) − H] 2 1/2 increases as W (L, t) ∼ t β for t < t s and becomes constant W (L, t) ∼ L χ for t t s . Here, .. denotes average over different noise realizations and the overbar is an spatial average in the x direction. The quantities χ, β and z are the roughness, growth and dynamical exponent, respectively, and they are related through the scaling relation χ = βz. In the standard Family-Vicsek scaling framework [25], this set of scaling exponents fully describe the scaling properties of the interface fluctuations.
However, experimental results and several growth models have appeared in the last decade showing that global and local scales are not always equivalent. This is called anomalous scaling [26,27]. Therefore, one has to compute the interface width at different window sizes, .. ℓ denotes an average over x in windows of size ℓ < L. For scaleinvariant interfaces local fluctuations are expected to increase as with the corresponding scaling function where χ loc is the local rough exponent and it characterizes the roughness at small scales. One of the implications of anomalous scaling is that the local width saturates when the correlation length reaches the system size, i.e. at the time t s and not at the local time t ℓ ∼ ℓ z as occurs in the Family-Vicsek scaling. There is an intermediate regime between t ℓ and t s where the local width grows as w(ℓ, t) ∼ t β * with β * = β − χ loc /z. Another useful technique to determine the whole set of scaling exponents is to study the power spectrum of the interface S(k, t) = h (k, t)h(−k, t) . In the presence of anomalous scaling it is expected to show the following scaling where the scaling function has the general form and χ s defines the spectral roughness exponent. Different scaling arises from this scaling function [26]. For χ s < 1 it is always valid that χ loc = χ s and the Family-Vicsek scaling is recovered whenever χ loc = χ. In contrast, the so-called intrinsic anomalous scaling is observed when χ s = χ. Note that for this kind of anomalous scaling a temporal shift in the power spectrum is observed. The quantification of this shift based on the above scaling form is unreliable and inaccurate, however, and thus we refrain from giving a measure for χ − χ s in the one case we observe scaling of this type. On the other hand, for χ s > 1 χ loc = 1 and the super-rough scaling appears when χ = χ s . It is worth to remember here that the crossover length ξ × , present in the imbibition phenomenon becomes a very important scale in the kinetic roughening process. More precisely, extensive numerical studies [5,20,21,22] have shown that fluctuations of the interface do not evolve in time beyond this crossover length. Therefore, the interface fluctuations reach the steady-state at t s ≃ ξ z × instead of saturating at L z . Additionally, since the interface is asymptotically flat at scales larger than the crossover length, the crossover shows up as a maximum at the corresponding wavevector in the saturated structure factor S(k, t → ∞).

IV. NUMERICAL RESULTS
Our study will be focused on the influence of the disorder on the scaling behavior of the fluctuating interface for the two cases of the mobility M in the phase field model, that is, the one-sided model and the symmetric case. In our simulations we have used gaussian distributed disorder with zero mean η = 0, and different disorder strengths (standard deviation of η) σ. In the numerical scheme the disorder will have a correlation length as long as the lattice spacing, meaning that the lattice spacing normalizes the standard deviation when dimensionless numbers in the numerical scheme are turned to physical units.
For the driven boundary condition, the methods used to obtain the interface description, and ultimately the LIE, predict that the mean value of the disorder will be irrelevant, as it has no contribution for the interface propagation. This has been numerically verified by our simulations of the full phase field models using different values for the disorder average. Here we only report results with η = 0.
As numerical method, we chose the simple explicit Euler scheme, where the disorder can be straightforwardly added. In two dimensions the timestep requirement of this method is not too restrictive for systems large enough for present consideration.

A. Weak disorder
For weak disorder, both models of the phase field (onesided and symmetric) are expected to be equivalent, since the same LIE describes both cases in this limit. The weak disorder regime corresponds to σ ≪ γ, where γ is the surface tension. In our dimensionless units (ǫ = 1 in Eq. (3)) γ ≃ 0.47. Numerically the roughness and growth exponents of χ s ≃ 1.3 and β ≃ 0.4, within the super-rough scaling description (χ loc = 1, χ = χ s ), were observed by integrating the LIE [22]. This is in agreement with our results for both phase field models at weak disorder, which are shown in Figs. 1 and 2 on the left. From Fig. 2(a) and (d) we observe χ = χ s , since no temporal shift in the structure factor is present (unlike Fig. 2(c)). Moreover, in Fig. 3 we observe that at small disorder strengths σ < 0.2 in accordance with σ ≪ γ, the spectral roughness exponent saturates to the LIE value, and is independent of the disorder strength.
We also find that numerical artifacts in the interpretation of the interface from the phase field model appear when the disorder is so small that the global interface width W (L, t) is much less than the numerical lattice constant. This is due to the interpolation required to obtain the location of the interface between two sites of the numerical grid, and shows up as prominent periodic oscillation in interface width with oscillation time ∆x/Ḣ 0 . While this numerical artefact can be removed by decreasing ∆x, it means that in the convenient and typically used [5,19,21,22] dimensionless units, which lead to Eq. (3) with ǫ = 1, only a very limited disorder strength range leads to the universal weak disorder limit. This range is roughly 0.1 < σ < 0.2, when ∆x = 1.

B. Strong disorder
In the regime of strong disorder, σ > 0.5, different scenarios occur depending on which model is used. Below, we discuss the cases of the symmetric and two-sided models separately.

Symmetric model
According to the analysis of the disorder strength above, in this limit it is favorable for the phase field model to spontaneously create disperse domains of one phase (A) within the region that is initially of the other phase (B). Droplets of phase A will form in phase B, and this mixture will initially cover most of the system. However, local mass conservation must still be valid regardless of any nucleation events. This means that mass must be diffusively transported from the A phase to the location where it nucleates within the B phase to facilitate a growing droplet.
This is exactly what happens in the symmetric model (see Fig. 4), when one imposes the initial condition of Eq. (5). Since there is no characteristic scale for the domain creation, the droplets are not restricted by the crossover length ξ × , which acts as a cutoff for the interface fluctuations and therefore, the surface roughening at large scales is different respect to a weaker disorder strength. This is observed in Fig. 2(c), where the interface power spectrum is plotted for a higher disorder strength at different times using the symmetric model. We can see that the fluctuations are not saturated, indicating that the crossover length ξ × (represented by the dashed line) is not acting anymore as an upper cutoff. In addition, our numerical results for the symmetric model at strong disorder show three differences to the weak disorder case. First, the local growth exponent β * approaches the global exponent β (see Fig. 1(c)). Second, the spectral roughness exponent decreases drastically to the range of χ s ≃ 0.5 (see Fig. 3) and third, a temporal shift appears in the power spectrum (see Fig. 2(c)). We can thus conclude that the scaling picture of interface fluctuations changes from superrough to intrinsic anomalous scaling, where χ s = χ.
The picture becomes problematic for strong disorder, however, because the interface becomes less and less representable by a single valued function H(x, t). This is due to the interfacial area becoming more and more tattered by overhangs, droplets and bubbles. This also means that some numerical tricks are needed to distinguish the interface from these bubbles and droplets. This distinction is essentially made by finding a path for the phase boundary across the system that locally has as small height jumps as possible. This works relatively well as long as the disorder is not much stronger than σ = 1 in our dimensionless units. However, note anomalous fluctuations in Fig. 1(c) around value 2.6 on the vertical axis, that are due to the abovementioned reason.

One-sided model
Using the one-sided model allows us to suppress the domain creation in phase B, where the mobility parameter is zero. Then, the position of the interface H(x, t) can be found by taking the largest height where the phase A has advanced to at time t, coming from the phase B when the phase field is above zero. In Fig. 4 we show an example of the interface profile at different times for a strong disorder, σ = 1.0.
The growth exponent β ≃ 0.5 measured for strong disorder strength (see Fig. 1(e) and (f)) agrees with the experimental value of β = 0.50 ± 0.02 reported in Ref. [16] for liquid front dynamics into a Hele-Shaw cell. Likewise, a similar variation in the spectral roughness exponent χ s , which changes from χ s ≃ 1.23 to χ s ≃ 0.91 when the disorder strength is increased (see Fig. 3), was also experimentally observed in the same reference, whith a variation from χ s = 1.1 ± 0.1 to χ s = 0.9 ± 0.1 when the capillary forces of the Hele-Shaw cell were increased (see Fig. 15 in Ref. [16]). On the other hand, we numerically observe that the crossover length ξ × still acts as a cutoff length for the interface fluctuations at strong disorder (see Fig. 2). These results indicate that the model can still describe the imbibition phenomenon at strong disorder.

V. CONCLUSIONS AND DISCUSSION
In this work, we have studied two different ways of considering the influence of the mobility parameter in a Model B type of phase field model with a Ginzburg-Landau type free energy. The main experimental context  considered here is liquid front invasion into a Hele-Shaw cell with quenched disorder [16]. We have focused on the case of driven front invasion, where there is a forced constant mass flux into the system that follows locally conserved dynamics. The symmetric model, studied for example in Refs. [5,21,22], uses a constant mobility factor, whereas the one-sided model, studied for example in Refs. [19,28], uses a mobility that is zero in the receding phase, which we call phase B.
We note that both models have previously been turned into non-local linear interface equations (LIEs) in the limit of small front fluctuations, which is equivalent to weak disorder [19,21,22]. These LIEs are identical for both models, and therefore both models are expected to have identical scaling behavior at the weak disorder limit. This is verified by direct comparison of numerical simulations. Furthermore, these results also agree with the relevant Hele-Shaw experiments [16].
We give an estimate for the strong disorder limit by comparing the disorder contribution to bulk energy to the surface tension. We find that the linear weak disorder limit is found well below this disorder value, and that only in this limit does the roughness exponent not continuously depend on the disorder strength in either model. This means that a well-defined region of universality only exists at the weak disorder limit.
Numerically we study the dependence of the growth and roughness exponents of the invasion fronts as a function of the disorder strength. Our results are consistent with a (continuous) change of scaling behavior from superrough to intrinsic anomalous scaling, when the disorder strength is increased from weak to strong.
At strong disorder, the symmetric model is no longer found to correspond to the Hele-Shaw experiment due to domain creation of the invading phase in front of the propagating interface. As our analysis shows, the domain growth can occur at all length scales (larger than the disorder site size l corr of the system) without any characteristic radius, which is a phenomenon observed in other experimental situations such as nucleation on dislocations [29] or binary mixtures [2]. In contrast, the results for the one-sided model do agree well with the Hele-Shaw experiments [16], even as the disorder strength is increased.
We hypothesize that the change in scaling behavior is due to decreased effect of surface tension and mass transport in interface roughening as the disorder becomes strong. Since the argument of strong disorder is the dominance of disorder over surface tension, the scaling in the regime dominated by the surface tension (characterized by having a correlation length ℓ c < ξ × ) should change for both models at strong disorder. In the one-sided model, mass transport still restricts interface roughening, and thus the crossover scale ξ × persists. Conversely in the symmetric model the nucleated domains create roughening by avalanches that are not controlled by mass transport from the reservoir, and thus the crossover scale ξ × becomes irrelevant.