The role of magnetization in phase-ordering kinetics of the short-range and long-range Ising model

The kinetics of phase ordering has been investigated for numerous systems via the growth of the characteristic length scale ℓ(t)∼tα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell (t) \sim t^{\alpha }$$\end{document} quantifying the size of ordered domains as a function of time t, where α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} is the growth exponent. The behavior of the squared magnetization ⟨m(t)2⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle m(t)^{2}\rangle$$\end{document} has mostly been ignored, even though it is one of the most fundamental observables for spin systems. This is most likely due to its vanishing for quenches in the thermodynamic limit. For finite systems, on the other hand, we show that the squared magnetization does not vanish and may be used as an easier to extract alternative to the characteristic length. In particular, using analytical arguments and numerical simulations, we show that for quenches into the ordered phase, one finds ⟨m(t)2⟩∼m02tdα/V,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle m(t)^{2} \rangle \sim m_0^2 t^{d\alpha }/V,$$\end{document} where m0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_0$$\end{document} is the equilibrium magnetization, d the spatial dimension, and V the volume of the system.


Introduction
Phase-ordering kinetics is omnipresent in nature and a very general phenomenon [1][2][3].It occurs in many physical situations, e.g., when quenching a system from a disordered state at high temperature to a low temperature where in equilibrium it would be in an ordered state.In case of spin models such as the Ising model, this corresponds to a nonequilibrium relaxation from the paramagnetic to a ferromagnetic state.During this process, one observes the growth of ordered regions (or domains) with time t which is driven by the reduction of the total area of domain walls, as every domain wall contributes energetically.The goal is to quantify this phase-ordering process.Usually this is achieved by investigating the characteristic length scale (t) as a measure for the average linear extent of the domains, which typically grows like a power-law (t) ∼ t α .We show that instead of measuring (t), one may alternatively investigate the conceptually much easier squared magnetization m(t) 2 during this process; a quantity apparently overlooked in the past as being interesting.
Henrik Christiansen and Suman Majumder contributed equally to this work.a e-mail: wolfhard.janke@itp.uni-leipzig.de(corresponding author) b e-mail: henrik.christiansen@itp.uni-leipzig.dec e-mail: smajumder@amity.edu We find m(t) 2 ∼ (t) d ∼ t dα , both from our analytical arguments as well as from simulations of the Ising model with nearest-neighbor interactions and the long-range Ising model with power-law interactions in one and two dimensions.This should not be confused with the zerofield-cooled or thermoremanent magnetization in the response of the system after quenching with an applied (small) magnetic field [4][5][6][7].While the equilibrium and nonequilibrium properties of the nearest-neighbor Ising model are, a priori, extremely well studied, both properties are not studied to the same extent in the long-range model.The Ising model with long-range interactions is of special interest, as it allows one to tune the interaction range, thereby serving as a generic model for long-range interactions that are omnipresent in natural and other sciences, ranging from electrostatic forces over neuroscience to economical phenomena [8][9][10][11][12].

Model and methods
The nearest-neighbor Ising model (NNIM) has the ubiquitous Hamiltonian where s i = ±1 are Ising spins, J 0 is the coupling strength, and i, j symbolizes summation over only the nearest neighbors of a d -dimensional (hypercubic) lattice.We use periodic boundary conditions to minimize boundary effects so that spins at opposing boundaries interact with each other as if they were nearest neighbours.As an extension to this model, we consider the long-range Ising model (LRIM) with Hamiltonian with J chosen to decay according to a power law, that is where d is the spatial dimension and σ > 0 can be tuned to capture different long-range behavior.For σ ≤ 0, the potential (without further normalization) diverges.This case will not be considered here.To mitigate the strong finite-size effects of this long-range model as much as possible, we will use Ewald summation to calculate effective interaction strengths J [13][14][15] over the minimum image convention often associated with periodic boundary conditions.Both models are investigated in d = 1 and d = 2 (on square L × L lattices).In our simulations J 0 sets the unit of energy and J 0 /k B the unit of temperature where k B is the Boltzmann constant (i.e., we set J 0 = k B = 1 for convenience).
The system is simulated using Monte Carlo (MC) simulations.One of the big advantages of MC simulations in equilibrium is that one may use collective updates in which many degrees of freedom are updated at once.Especially successful for spin models with both nearest-neighbor and long-range interactions are cluster algorithms [14,[16][17][18][19][20].For nonequilibrium studies, it is, however, important to only perform physical and thereby local moves, so that one is restricted to singlespin flip dynamics.The computational complexity for a single spin flip scales as O(1) in the NNIM, but as O(N ) in the LRIM since all N ≡ V = L d spins interact with each other.Exploiting the observation that during a sweep (which consists of N single-spin flip attempts) only a small fraction of spins is updated, we were able to significantly lower the prefactor of this scaling by a factor ∼ 10 3 [15,21].This allows us to simulate rather big systems of size N up to 4096 2 , which was sufficient for the present purpose.For a recently developed much faster but more involved hierarchical and adaptive tree-like single-spin update algorithm that is compatible with O(log N ) complexity, see Refs.[22,23].
The simulation protocol in our phase-ordering kinetics study is to prepare an initial disordered configuration with magnetization m = i s i /V ≈ 0. Subsequently, we quench the system below the critical temperature T c , where the system is ferromagnetic.The process occurring is driven by the reduction of energetically unfavorable domain boundaries, where eventually at late times, one of the two magnetic branches with m fluctuating around ±m 0 is "randomly" selected with m 0 denoting the equilibrium magnetization at the quench temperature.

Results
We start our analysis by presenting in Fig. 1 snapshots during the phase-ordering process in d = 2 of the NNIM compared to the LRIM with σ = 0.6 for a quench from a randomly drawn disordered configuration at T = ∞ to T = 0.1T c , where T c (σ) for the LRIM was extracted from Ref. [14].Phenomenologically, the events appear similar and in both systems, the domains grow with time.While in the NNIM, the domains form with very smooth boundaries, for the LRIM, one observes more often relatively sharp corners, especially at late times.Apparent are of course the different time scales involved.For the NNIM, the lattice is close to be aligned at t ≈ 10 5 (in units of sweeps), whereas in the LRIM, the same is observed already at t ≈ 10 3 .There could be different reasons for this: (i) a larger amplitude in the LRIM, (ii) a different growth exponent, or (iii) a combination of both.

Characteristic length
Therefore, in phase-ordering kinetics, one is mainly interested in the quantification of the growth of domains consisting of ordered regions with time.This can be quantified by measuring the characteristic length scale (t) of the average domain size.It is known that for the NNIM one has a power-law growth with α = 1/z = 1/2, where z = 2 is the usual equilibrium dynamical exponent [24,25].For the nonequilibrium behavior of the LRIM, there exists a long-standing prediction for the growth of the characteristic length based on continuum descriptions of this model [26][27][28]: where for σ > 1, one has the same growth exponent α = 1/2 as in the NNIM.This scaling law does not depend on the spatial dimension d and has been numerically confirmed by us in d = 2 [15,21] and subsequently also settled to be valid in d = 1 [29].
In our studies, we extracted the values of the characteristic length from the decay of the equal-time twopoint correlation function where . . .stands for an average over different initial (disordered) configurations and time evolutions and r is the distance between the locations of the spins s i and s j .To get an appropriate estimate of (t), we take the intersection of C(r, t) (= C( r, t) averaged over all directions) with a constant value c ∈ (0, 1) (we choose c = 0.5).This method is rather involved and has the drawback that for an efficient calculation of the correlation function, one has to make use of a fast Fourier transform [30,31] to calculate s i s j .In Fig. 2, we show the characteristic length for the d = 2 LRIM with σ = 0.6, 0.7, 0.8, and 1.5.The system was prepared in a disordered configuration and quenched to T = 0.1T c (σ) as explained above.Here, we present new data for large systems of linear size up to L = 4096 and average over 50 independent time evolutions.For the values of σ that are in the longrange regime, a growth with α = 1/(1 + σ) is expected, whereas for σ > 1, one expects α = 1/2.As can be seen from the log-log plots of this figure, the data are consistent with this prediction depicted by the black solid lines.This is a reconfirmation of our results previously reported for up to L = 2048 [15,21].

Magnetization
In the paramagnetic phase in equilibrium, one trivially has m(T ) = 0 since the entropic contributions dominate, where m = i s i /V is the usual magnetization.Below the critical temperature, in the ferromagnetic phase, the spontaneous magnetization m 0 (T ) = 0, i.e., one in practice observes |m(T )| ≈ m 0 .However, due to the Z 2 symmetry, for finite systems in equilibrium, one still in principle has m(T ) = 0.In practice, this is often not observed in equilibrium simulations for large systems since there is a first-order transition between the magnetic branches below T c suppressing the crossover.During a quench from the paramagnetic phase to the ferromagnetic phase, m of finite systems goes for each individual run from the initial magnetization m ≈ 0 in the disordered state1 into one of the two ordered branches with m = ±m 0 (T ) (≈ ±1 for low T ), where the sign is realized randomly for each run.This means that for the time evolution of the expectation value of the magnetization, one has m(t) = 0. What is nonzero, however, is the expectation value of the squared magnetization m(t)2 as a function of time since each individual run still evolves into one of the magnetization branches.
In Fig. 3, we present these individual time evolutions of the magnetization following a quench in the ordered region.In total, we have run 500 independent time evolutions for the d = 1 NNIM, 150 for the d = 1 LRIM, and 50 for the d = 2 NNIM and LRIM, respectively.The d = 1 model is shown in (a) for the NNIM and in (b) for the LRIM with σ = 0.6.Plots (c) and (d) are the analog plots for d = 2.The quench temperature was generally chosen to be T = 0.1T c with the "exception" of the d = 1 NNIM in (a) where this choice degenerates to T = 0. 2 It is apparent that in most cases, roughly 50% of the runs choose a positive magnetization, whereas the other 50% are negatively magnetized.The trajectories for all four different systems look very different.While the time evolutions of the magnetization for the d = 1 NNIM more or less resemble a "random walk", the curves for the corresponding LRIM look much more directed.In some of the simulations, configurations occur that get "stuck" into a (strictly speaking meta-stable) striped state configuration as can be appreciated from the curves that are roughly constant at m = ±1.This is a well-known phenomenon for the d = 2 NNIM (especially at zero temperature) [34] and has recently also been investigated in the LRIM [35].We now present an analytic argument for the growth of the expectation value of the squared magnetization m(t) 2 , which is clearly nonzero from above data.Since we know that the average linear domain size scales with time t as (t) ∼ t α , the average volume of a domain scales in d dimensions as so that the average number of domains in a total volume V = L d is given by where the scaling function f (x ) satisfies ∞ 0 dxf (x) = 1 and x ≡ ∞ 0 dxxf (x) = 1.The first condition implies that P ( d , t) is normalized to unity and the second condition leads to the desired property that Higher-order moments scale similarly with (t) as The total magnetization M (t) = V m(t) results by summing up the N d (t)contributions where the sign ±1is random and it is assumed that each domain has already equilibrated locally to the spontaneous magnetization density m 0 (T )at T < T c .Due to the random sign, M (t) = N d i=1 M d of each time evolution is equally likely positive or negative so that M (t) = 0.For the squared magnetization, all contributions add up and by a random-walk argument (with M (t) playing the role of the end-to-end distance), one obtains This result confirms that the squared magnetization is zero in the limit V → ∞, as one would expect for a global quantity.It is the amplitude of the leading finite-size scaling correction that carries the signature of the growing length scale.Note that a similar power-law scaling behavior is known from the critical short-time dynamics following a quench from a disordered initial state with vanishing magnetization to the critical point of a system: m(t) 2 ∼ t (d−2β/ν)/z /V where β and ν are the critical exponents of the magnetization and correlation length, respectively, and z is the dynamical critical exponent [36,37].Putting formally β/ν = 0 for a noncritical quench to temperatures below the critical temperature and identifying 1/z = α with its non-critical value, Eq. ( 14) is recovered.Even though the physical situations are quite different, the origin of the similar scaling behavior can be traced to the applicability of the random-walk argument (or central limit theorem) in both cases.
The simplest but unrealistic realization of the dynamic scaling function in (9) would be f (x) = δ(x−1) (with x n = 1) which corresponds to assuming that all N d domain sizes are equal, d = (t), so that the random-walk picture assumes its simplest form.A pure exponential, f (x) = e −x (with x n = n!), is clearly a more realistic ansatz.Empirically an even better choice [38] is the double exponential f (x) = a(e − 2ax a+1 −e − 2ax a−1 ), parameterized by the constant a.For a = 2, it fits numerical data for the two-dimensional nearestneighbor Blume-Capel model very well [38].For this parameter choice, f (x ) steeply increases for small x as 16x /3, develops a peak at x max = (3/8) ln 3 of height f max = 4/(3 √ 3), and for large x decays exponentially as 2e −4x/3 .Recall that the Blume-Capel model is a spin-1 model with s i = ±1 and 0, where the vacancies s i = 0 can be suppressed by a large crystal-field coupling so that in this limit the model becomes equivalent to the Ising model.
A simpler, less microscopic argument goes as follows.The dynamical scaling properties of the coarsening process imply in general that C( r, t) depends only on the scaling variable r/ (t) but not on r and t separately.This in turn implies that the structure factor, its Fourier transform S( k, t) = d rC( r, t)e i k r , exhibits the dynamic scaling behavior S( k, t)/ (t) d = f ( k (t)).Since this does not depend on t as k → 0 and a direct calculation yields the usual result S( k = 0, t) = V m(t) 2 (recall that m = s i = 0), one concludes that m(t) 2 ∼ (t) d /V , in agreement with the prediction (14) based on microscopic arguments.
We now check the prediction in Eq. ( 14) numerically: The upper two subplots of Fig. 4 show V m(t) 2 versus t in d = 1 after a quench to (a) T = 0 for the NNIM and (b) T = 0.1T c for the LRIM with σ = 0.6.The other subplots (c) and (d) show the analogous data for d = 2. Here, instead of m(t) 2 , we have plotted V m(t) 2  to check for the predicted prefactor in (14).The factor m 2 0 is non-universal and dependent on temperature.We find in all cases that the data for different system sizes collapse reasonably well onto each other, confirming the prefactor of the prediction.In particular they agree well with the predicted time dependence drawn as a solid line with the theoretical slope (and adjusted amplitude).Furthermore, adapting the recipe of Ref. [39], we have verified that least-square fits using the ansatz V m(t) 2 = at dα provide estimates for the growth exponent α that are compatible with the theoretical expectations within a few percent accuracy.Details of this quite elaborate fitting exercise are discussed in the Appendix.
The big advantage of analyzing m(t) 2 is of course that this is extremely easy to calculate compared to the rather involved estimation of (t) via the correlation function.When compared to the estimation of (t) via counting the spin sign-changes it is not (overly) sensitive to single thermally excited spins, either.One of the drawbacks is, however, that this quantity has significantly larger error bars.

Conclusion
We have investigated the behavior of the (squared) magnetization following a quench from a disordered state into the ordered phase, which is an often overlooked quantity during coarsening.While for finite systems of size V = L d , the average magnetization is zero, the average squared magnetization carries the signature of growing domains.Using analytical arguments, it is shown that m(t) 2 ∼ (t) d /V , i.e., the squared magnetization may be used as an alternative observable to quantify the coarsening of domains.The predicted time dependence is verified numerically in d = 1 and d = 2 both for the NNIM and the LRIM with σ = 0.6.We also confirm the prefactor 1/V of the growth, validating our prediction.Additionally, we have also presented new data for the growth of the characteristic length (t) itself for the d = 2 LRIM with σ = 0.6, 0.7, 0.8, and 1.5 for system sizes up to L = 4096 using for (t) the standard definition from the intercept of the equal-time two-point correlation function (with c = 0.5).It would be interesting to investigate the time dependence of the squared magnetization in other coarsening systems such as the XY model both in the short-range and long-range variant.For both considered quantities, one expects corrections to asymptotic scaling for small times t and finitesize corrections for large times.Since these corrections are not taken into account in the fit ansatz, we have estimated their influence on the parameter estimates by systematically varying both ends of the fit interval t ∈ [t min , t max ] over a rather wide range.
In the following, we go through the data presented in Fig. 4a-d point by point:

(b) 1D LRIM (σ=0.6):
A glance on the data for the 1D LRIM in Fig. 4b shows that here the situation is more difficult.In particular, the finite-size corrections creep in much more gradually than for the 1D NNIM and even for the large lattice size L = 32768 they obviously start to become important already for t ≈ 5 × 10 3 .This limits the available asymptotic fitting range substantially.A quite conservative fit interval is t ∈ [60, 3000], still covering more than 1.5 decades in t.This gives α = 0.62(3), ln(a) = 3.7(2), a = 41 (7) with an acceptable χ 2 /dof = 0.12 for 14 dof.

(d) 2D LRIM (σ=0.6):
For the 2D LRIM, the data for the largest lattice size L = 4096 shown in Fig. 4d looks somewhat off the trend for the smaller lattices.We have checked that this behavior is reproduced when averaging over different start configurations and time evolutions (and also when employing a completely independent computer code).

Fig. 3 1 −
Fig. 3 Magnetization m versus time t in d = 1 for a the NNIM and b the LRIM with σ = 0.6, whereas c and d show the analogous data in d = 2.For the d = 1 NNIM we used T = 0 and for all other cases T = 0.1Tc as our quench temperature.For details, see text.Note the (very) different time scales of the plots

Fig. 4 5 10 1 10 2 10 3 10 4 10 5 10 6 10 7 (
Fig. 4 Shown is V m(t) 2 versus t after a quench into the ordered phase in d = 1 for a the NNIM and b the LRIM with σ = 0.6.The other panels c and d show analogous data for d = 2

Table 1
Summary of fitted growth exponents α using the ansatz V m(t)2= at dα respectively (t) = at α (last line) The numbers in square brackets[ . ..]  indicate the systematic error due to the choice of the fit interval (read off from the last 3 columns).χ 2 r stands short for χ 2 /dof interval, we finally opted for this latter choice using 10 t values per decade (we have explicitly checked that using 20 or more t values per decade does not significantly alter the parameter estimates).All fit results reported below refer to this case if not otherwise indicated.Another advantage of the logarithmic spacing is that the used t values are naturally widely spaced (in particular for larger t), so that correlations of the data for closeby times do not play such a prominent role (of course, for a high-precision study, one could deal with such correlations by jackknife blocking techniques which are, however, rather laborious and hence left for future work).