Phase separation and critical percolation in bidimensional spin-exchange models

Phase separation is the process whereby a binary mixture of components A and B, initially in a homogeneous phase, demix. This process leads to the coexistence of two phases: one rich in A and the other in B [1–6]. The system, initially in an unstable spatially uniform state, progressively coarsens to approach its thermodynamically stable phase-separated state. Such phenomena arise in binary alloys, fluid mixtures, and polymer blends. Recently, the dynamics of phase separation have seen a revival of interest in the context of experimental [7,8] and numerical [9–12] studies of binary mixtures of Bose gases.

The late time dynamics are well understood. In the absence of driving forces, a dynamic scaling regime with statistically self-similar domain morphology sets in. This regime is well described by an extension of the Lifshitz-Slyozov-Wagner (LSW) theory [13,14], in which the typical domain radius grows as [15]

$$\begin{equation} \ell_d(t) \simeq t^{1/z_d}\qquad \text{with} \qquad z_d=3 \end{equation} \tag{ 1 }$$

(whereas for scalar non-conserved order parameter dynamics the growing length is also given by a power law but the exponent is $z_d=2$). Numerical results in favour of this law were published in [15–17] for spin-exchange models although the growth law can be more complex in particle or polymer phase-separating systems, see, e.g., [18] and references therein. The pre-asymptotic dynamics leading to this regime have not been discussed in detail in the literature.

It was noticed in [19] that the low-temperature evolution of a bidimensional 50:50 binary mixture after a quench from infinite temperature shares many points in common with the one generated by Glauber single spin-flip stochastic dynamics satisfying detailed balance [20,21]. On the one hand, an early approach to critical percolation was noticed, although the time needed to reach this state was not studied in detail. On the other hand, a separation of length scales in the statistics and morphology of finite-size cluster areas and domain wall lengths was observed. Linear or planar objects that are smaller than the typical ones, $\ell_d(t)$ or $\ell_d^2(t)$, satisfy dynamic scaling with respect to $\ell_d(t)$, while larger objects were found to be very close to the ones of critical percolation.

In this letter we characterise the early stages of the dynamical process. More precisely, we analyse the way in which the system approaches a state with a stable pattern of critical percolating domains. We monitor a number of observables (to be defined in the main part of the text) and we explain how their behaviour constitutes evidence for this claim. We prove the approach to critical percolation for balanced mixtures whereas different behaviour is found for asymmetric ones [11].

One of the quantitative goals of this work is to measure the growing length, $\ell_p(t)$, that characterises the approach to critical percolation for phase-separating systems on the lattice. Accordingly, we define the time needed to reach critical percolation, t_p, from $\ell_p(t_p) = L$ with L the linear size of the sample. For Glauber dynamics the approach to the percolating state is characterised by $\ell_p(t) \simeq t^{1/z_p}$, with z_p < z_d, and the numerical analysis suggested [22]

$$\begin{equation} \ell_p(t) \simeq \ell^n_d(t) \end{equation} \tag{ 2 }$$

with n the coordination of the lattice. The results that we will present here indicate that such an algebraic relation is also satisfied for the local spin-exchange (Kawasaki) dynamics of binary mixtures with equal concentration of the two species. For Glauber dynamics, reasoning in terms of z_p and z_d gives very good results, since $z_d=2$ over a long time interval. However, for locally conserved order parameter dynamics the power law is hard to establish and the representation of $\ell_p$ in terms of $\ell_d$, as in eq. (2), is better. As for the exponent n, we argued that, in Glauber dynamics, it is the coordination of the lattice that we can associate to half the maximal possible energy change induced by a single spin flip [22]. For local Kawasaki dynamics neighbouring pairs of spins are updated simultaneously. Such a change can induce, at most, an energy change equal to twice the number of nearest neighbours of a neighbouring pair of sites. If we follow the argument for the Glauber case, n = 6 for the triangular lattice (we built this lattice in such a way that there are 8 nearest neighbours of a pair of sites but two of them are shared by the two spins in the pair), n = 6 for the square lattice, and n = 4 for the honeycomb lattice. Note, however, that the triangular lattice is special as the 50:50 initial conditions are right at the critical percolation point. We will discuss this guess in the body of this letter.

Concretely, we study a lattice gas model with local Kawasaki rules [23–26] used to mimic phase separation in systems in which hydrodynamic effects can be neglected [2]. We use the spin language, in which up and down Ising variables correspond to the presence of the A and B species on a given site. The energy function is the familiar ferromagnetic Ising one,

$$\begin{equation} H = -J \sum_{\langle ij\rangle} s_i s_j, \end{equation} \tag{ 3 }$$

with the sum running over pairs of nearest-neighbour sites on the lattice and J > 0. We set J = 1 so that the critical temperature on a square lattice is $T_c = 2 \,[\ln{(1+ \sqrt{2})}]^{-1}$ in units of k_B. The spin-exchange (Kawasaki) dynamic rules are defined as follows. At each time step, a pair of nearest-neighbour sites is chosen at random. If the spins sitting on these sites are antiparallel they are exchanged with the heat bath Monte Carlo rule. If the two sites are occupied by the same kind of spin their state remains unchanged. The control parameters are temperature and the relative concentration of the two species. All data shown are for the Ising model on a square lattice with periodic boundary conditions after quenches from infinite temperature to $T_c/4$ and $T_c/2$. Results on other lattice geometries will be mentioned but not shown.

In fig. 1 we show the evolution of the characteristic domain pattern in a system with linear size L = 128. The concentration of up (red) and down (white) spins is a half. The clusters that percolate in at least one direction are highlighted in green and blue, respectively. Very early, already at t = 4, a percolating cluster appears, then breaks (not shown) and rebuilds again $(t=128)$ until in the late-time snapshot, at $t \geq 512$, two large clusters of opposite orientation are interlaced and percolate horizontally. This configuration belongs to the first class sketched in fig. 2, named "h" for horizontal along the torus. Note that other runs can lead to configurations of the "v" and "diag" type with also two percolating clusters or "hv" with only one percolating cluster. We will show the probability of reaching each of these in fig. 6 and discuss them in the text below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A way to estimate the time dependence of the dynamic growing length involved in the dynamic scaling is to use the time-dependent excess energy, $E(t) = \langle H(t)\rangle$, with respect to the ground-state energy E₀

$$\begin{equation} \ell_{G}(t) = E_{0} /[E_{0} - E(t)]. \end{equation} \tag{ 4 }$$

(e.g., $E_{0}= - {2}L^2 {J}$ for the square lattice). The angular brackets indicate an average over different realisations of the dynamics (initial conditions and/or thermal noise). This quantity is shown in fig. 3 for Kawasaki and Glauber dynamics at two sub-critical temperatures, $T_c/2$ and $T_c/4$. Contrary to what happens for Glauber dynamics, nowhere in the time span shown in the figure a stable algebraic increase of $\ell_G$ established. The evolution of the effective exponent in $\ell_G(t) \simeq t^{1/z_{\rm eff}(t)}$ is followed in the inset. In this time window $z_{\text{eff}} (t)$ goes roughly from 10 to 4 for $T_c/2$, and the measurement slowly approaches the expected value for the dynamic exponent $z_d =3$ [15] but time scales of two orders of magnitude longer are needed to reach convergence [16,17]. For other geometries, such as triangular or honeycomb lattices, the values of $z_{\text{eff}}$ are similar to the ones for the square lattice [27]. The fact that the effective exponent varies so much in time and depends on temperature suggests to use the time-dependent growing length itself to analyse the pre-asymptotic regime with the eventual approach to percolation.

Zoom In Zoom Out Reset image size

The first observable that we study to determine whether the coarsening dynamics drives the system towards a critical percolating state is the winding angle, $\theta(x)$, of an interface. It is defined as the angle between the local tangents at two points separated by a curvilinear distance x along the curve (for a spin configuration on a lattice the local tangent to an interface is a perpendicular vector to a broken bond and its direction takes values on a discrete set). Bidimensional critical interfaces are, in the continuum limit, conformally invariant curves described by the stochastic Loewner evolution $\text{SLE}_\kappa$. The average squared $\theta(x)$ can then be exactly computed [28,29]:

$$\begin{equation} \langle \theta^2(x) \rangle = \text{cst} +{4 \kappa \over 8+\kappa } \ln{x}. \end{equation} \tag{ 5 }$$

The parameter κ is related to the interface fractal dimension and determines the universality class. For critical percolation interfaces, $\kappa = 6$, while for critical Ising ones $\kappa=3$.

After the quench the averaged variance of the winding angle approaches the form (5) with $\kappa\simeq 6$ at a time t_p, and the later evolution conforms with dynamic scaling with the dynamic growing length. As a proof of this claim, we show in fig. 4 $\langle \theta^2(x,t) \rangle$ for the largest cluster interface in a system with L = 640. The dotted straight line is a fit to the data at t = 8192 and yields $\kappa \simeq 6.831$. We then used the rescaled curvilinear coordinate $x/\ell_G(t)$ to successfully scale the data at subsequent times, as shown in the inset. We also used the data for $\langle \theta^2(x,t) \rangle$ to obtain an alternative estimate of $1/z_{\mathrm{eff}}(t)$. We looked for pair-wise collapse of datasets $\{\langle \theta^2(x,t_i) \rangle\}$ at consecutive times with $1/z_ {\mathrm{eff}}(t_i)$. The outcome is shown with open (green) points in the inset of fig. 3.

Zoom In Zoom Out Reset image size

We now turn to the analysis of the time scale t_p and the growing length towards critical percolation $\ell_p(t)$. In [30] we used the asymptotic value of the overlap $Q(t,t_w) = N^{-1} \sum_{i=1}^N s_i(t) \sigma_i(t_w)$ between two copies of the system, $\{ s_i(t), \sigma_i(t)\}$, created at a waiting time $t_w, \, \{s_i(t_w) = \sigma_i(t_w) \}$, that later evolve independently to measure $t_p(L)$ in an Ising model evolved with Glauber dynamics. The idea was to scale t_w with L and search for the weakest L-dependence in the form $t_w = L^{z_p}$ such that $\lim_{t\gg t_w} Q(t,t_w(L)) >0$ (while it vanishes otherwise). This yields $t_p = t_w(L) =L^{z_p}$. In the conserved order parameter case, the excess energy growing length does not reach a stable power law regime and, consequently, it is hard to use a fixed power of L to quantify t_p. For this reason, we used other observables to estimate $t_p(L)$.

In fig. 5 we study the geometric and scaling properties of the largest cluster. With the use of the area and interface fractal dimensions of critical percolation clusters [31–33] we find good scaling of its area and interface length using $A_c(t)/L^{D_A}$ and $l_c(t)/L^{D_\ell}$ against $\ell^n_{\text{G}} (t)/L$ with n an adjustable parameter. The best collapse was found using n = 5, a value that is not far from the guess n = 6 on the square lattice. We performed the same analysis on other lattices and we report here the results without showing the scaling plots. On the honeycomb lattice we found $n_{\text{honey}}=6$ from the scaling of both A_c and l_c. Instead, on the special triangular case, we could not scale the data for A_c while we found acceptable scaling of the data for l_c using $n_{\text{triang}}= 6$.

Zoom In Zoom Out Reset image size

These results are confirmed by the analysis of $\pi_{\text{hv}}$, $\pi_{\text{h}} + \pi_{\text{v}}$ and $\pi_{\text{diag}}$ [34] defined as the probability of having a cluster wrapping in both directions of the lattice, in one Cartesian direction only (horizontally or vertically) and diagonally (as in the fourth sketch of fig. 2), respectively. These wrapping probabilities have been calculated for critical percolation on a torus in [35] and are shown with dotted lines in the upper panel in fig. 6. The dynamic data, also shown in this figure against the scaling variable $\ell_G^n(t)/L$, approach these values asymptotically for sufficiently large system sizes. Consistently with what we found in the analysis of A_c and l_c, the best scaling is obtained with n = 5. A similar analysis on the triangular and honeycomb lattices yield $n_{\text{triang}}=n_{\text{honey}}=6$. In all cases other choices of the parameter n give scaling plots that are of considerable less quality.

Zoom In Zoom Out Reset image size

Until now we used an equal concentration of up and down spins. We now search to determine whether the same critical percolation phenomenon exists for different values of, say, the up spin concentration p. In the lower panel in fig. 6 we display the various wrapping clusters probabilities for $p=0.4$, 0.42, 0.45, 0.47, 0.5. The probability of having a wrapping cluster along one principal direction of the lattice, $\pi_v + \pi_h$ (green curves), does not reach the dotted horizontal line for $p = 0.4 < 1-0.5927$ since one of the two species has a percolating cluster initially. For $p=0.42$, 0.45, 0.47, none of the two species percolates initially, and the (green) curves increase and approach the critical percolation value but they rapidly detach from it and decrease to zero. The complementary (red) curves simultaneously approach one. The curves for $p = 1/2$ are the only ones that approach the non-trivial asymptotes shown with horizontal dotted lines. This result is consistent with the observation in [11], where the segregating dynamics of a mixture of Bose-Einstein condensates was studied.

We stress the fact that t_p is just a characteristic time scale associated to the approach to the critical percolation state, mostly used in the form $t_p \simeq L^{z_p}$ to scale time as done quite successfully in [30]. However it can be useful to have an idea of the order of magnitude for the system sizes used. One criterium to get a numerical estimate is to look for the time t at which the wrapping probabilities shown in the top panel of fig. 6 reach approximately the asymptotic values corresponding to critical percolation. This condition corresponds roughly to $\ell^n_G(t_p) / L = 10$ and yields $t_p = 1100$, 2400, 4900, 9400, 17200 MCs for $L= 40$, 80, 160, 320, 640, respectively. This relation is approximately linear suggesting $z_p =1$. This value is in acceptable agreement with the guess $z_p = z_d/n$, if we take into account the fact that the effective z_d is $z_{\text{eff}}\simeq 5$, see fig. 3, and the measured n is close to 5.

In fig. 7, we study the time-dependent pair connectedness correlation function $g(r,t)$ against distance r at different times t for a system with L = 640. This function measures the probability that two occupied sites (say, with spins up) separated at distance r belong to the same cluster [31,33]. We compare the dynamic results to the ones of random site percolation at the site occupation probability $p=0.5927$ that is approximately the threshold value p_c (red dashed curve), a case for which $g(r) \simeq r^{-2 \Delta}$ where $\Delta = 2 -D_A$ for $a\ll r \ll L$ with a a microscopic length scale. Data are for square lattices with the same linear sizes. In the inset of the same figure we show $g(r,t) \cdot [r/t^{1/ z_{\text{eff}}}]^{2\Delta}$ against r, where we take $1/z_{\text{eff}}$ as an adjustable parameter. Using $1/z_{\text{eff}}\approx 0.27$ the dynamic data collapse onto the curve $g_{\text{perc}}(r) \cdot r^{2\Delta}$ of critical percolation (similarly, for L = 160 we find $1/z_{\text{eff}}\approx 0.28$). The value of $1/z_{\text{eff}}$ is consistent with the effective growth exponent measured from the excess energy, see the inset in fig. 3. At short length scales, the clusters emerging from the dynamical process are very different from the ones of critical percolation. The percolation and dynamic curves are not flat at large r due to the periodic boundary conditions.

Zoom In Zoom Out Reset image size

The time-dependent number density of geometric domain areas A, that we call ${\mathcal N}(A,t, L)$, should scale as

$$\begin{equation} {\mathcal N}(A,t, L) \simeq N(A,t) + N_p(A,t,L) \end{equation} \tag{ 6 }$$

with the first term describing the weight of the finite areas and the second one the weight of the areas that span the sample. At t_p the last term should become independent of time, and scale with $A/L^{D_A}$. A discussion of the scaling properties of the first term in terms of $A/\ell_d^2(t)$, and the crossover of the scaling function from the regime in which the scaling variable varies from being much smaller to being much larger than one was presented in [19]. In the second regime the first term approaches the algebraic finite-size clusters distribution

$$\begin{equation} N(A)=2 c_d \ A^{-\tau_A} \end{equation} \tag{ 7 }$$

with $\tau_A = 1 + d/D_A$ [31–33]. There is no analytic prediction for c_d but a rough estimate suggests $c_d \simeq 0.029$ [19,21,36].

We first focus on the large area regime of the number density of non-wrapping domain areas, that is to say, the contribution $N(A,t)$ to the complete number density ${\mathcal N}(A,t,L)$, see eq. (6). The data for a system with L = 640 are presented in fig. 8, where $A^{\tau_A} \, {\mathcal N}(A,t,L)$ is plotted against A in the upper panel and against $A/\ell^{n{D_A}}_G(t)$ in the lower panel with n = 5. In both cases the critical percolation value $\tau_A=187/91$ was used.

Zoom In Zoom Out Reset image size

In this representation the flat part is close to 2c_d (see [27] for more details on the normalisation), the vertical looking parts of the curves belong to N_p, and the remaining bump corresponds to a power law decay of N with an exponent smaller than $\tau_A$. This regime is well described by the scaling function

$$\begin{equation} A^{\tau_A} \, {\cal N}(A,t,L)\sim \Phi ( A/\ell^{n{D_A}}_G(t)), \end{equation} \tag{ 8 }$$

where $\ell_G^{n{D_A}}(t)$ has the meaning of a characteristic domain area at time t, and Φ is a scaling function such that $\Phi(u) \rightarrow 2c_d$, as $u \rightarrow 0$, and $\Phi(u) = C \, u^{a}$, with a > 0. By fitting this function to the data at early times (up to t = 32), in the region [1, 10²] of the scaling variable $A/\ell^{n{D_A}} _G(t)$ with n = 5, we find $C\simeq 0.101$ and $a \simeq 0.304$. The same scaling behaviour with, interestingly enough, the same scaling function is found for Glauber dynamics [27]. The collapse of the data for the number density of domain areas is less sensitive to the choice of the value of n than the other observables studied. For this and the other lattices we found acceptable scaling for n between 4 and 6, approximately.

It is also instructive to study the shape of the probability distribution of the largest clusters. In percolation, the distribution of the largest cluster area is a Gumbel distribution for p < p_c [37–39] while it is approximately Gaussian for p > p_c [40]. Little is known about its distribution at p_c, except for a remarkable exact result in the mean-field case [41]. Numerical studies [42] suggest that there is a smooth crossover between the subcritical and the supercritical phase, and that the probability distribution of the order parameter can be approximated by a weighted sum of a Gumbel and a Gaussian distribution. The next size in the hierarchy is much smaller than the largest one.

We reckoned in fig. 1 that two very large clusters with size comparable to the system size develop dynamically. In fig. 9 we follow the time evolution of the size distribution of the two largest clusters, $\mathcal{N}_ {\mathrm{(LC+SLC)}}$, and we compare it to the one of the largest cluster in random site percolation at $p=0.5927$. At the late stages of the approach to percolation, $t \lesssim t_p$, the distribution of the two largest cluster sizes (after proper rescaling) resembles more and more a Gumbel distribution. Beyond t_p the distribution approaches a Gaussian distribution (not shown).

Zoom In Zoom Out Reset image size

In conclusion, the ordering dynamics of a phase-separating 50:50 mixture starts with an approach to critical percolation that lasts for a macroscopic time t_p determined by $\ell_p(t_p) \simeq L$, when a stable pattern of percolating domains establishes. After this time, the percolating cluster(s) become fatter. This second regime is characterized by the expected growing length $\ell_d(t) \simeq t^{1/z_d}$ with $z_d=3$. The growing length of critical percolation clusters $\ell_p(t)$ is related to the dynamic growing length as determined from the excess energy or other observables via eq. (2) with the exponent n being numerically close, though not identical, to the number of neighbours of a pair of nearest neighbours on the lattice. The approach to percolation and the proper dynamic scaling regime with $z_d=3$ are well separated, similarly to what happens in Glauber dynamics [30] (where, say on the square lattice, $z_p=1/2$ and $z_d=2$) but differently from what was found for the voter model [43] (where $z_p \simeq 1.67$ and $z_d=2$). For unbalanced mixtures in which one species is more present than the other, the dynamics does not reach a long-lasting critical percolation pattern.

We have also studied non-local Kawasaki rules in which the pairs of spins updated are not necessarily nearest neighbours finding similar results to the ones discussed here. A detailed comparison with this and other (Glauber, voter) dynamic rules will be presented in a long publication [27].

We thank F. Corberi and F. Insalata for very useful discussions. LFC is a member of Institut Universitaire de France.