- Email: [email protected]
Quantifying nonequilibrium behavior with varying quenching rates
Physica D 181 (2003) 121–131
Quantifying nonequilibrium behavior with varying quenching rates Carmen J. Gagne∗ , Marcelo Gleiser Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA Received 7 May 2002; received in revised form 29 January 2003; accepted 27 February 2003 Communicated by L. Kramer
Abstract We investigate the approach to thermal equilibrium of (1 + 1)-dimensional stochastic Ginzburg–Landau models at varying cooling rates. The nonequilibrium dynamics is modeled by coupling the field to an external heat bath with damping rate η. We argue that the departure from thermal equilibrium can be measured from the absolute value of the rate of change of the momentum-integrated structure function, (t). If the field equilibration time-scale is faster than the cooling (or quench) time-scale τq , then (t) → 0. Otherwise, (t) displays a peak which scales as τqn , with 1/3 ≤ n ≤ 1/2 as η is varied from the underdamped to the overdamped limit. Furthermore, we show that the amplitude of the peak in (t) measures the departure −6/5 from equilibrium and scales as τq , independent of η. © 2003 Elsevier Science B.V. All rights reserved. PACS: 05.70.Ln; 11.10.Wx; 98.80.Cq Keywords: Nonequilibrium dynamics; Stochastic PDEs; Ginzburg–Landau models; Finite quench rate
1. Introduction Given the widespread role of phase transitions in condensed matter physics, cosmology, and high energy physics, it is imperative that the nonequilibrium dynamics of these phase transitions be better understood. In particular, the effect of varying cooling rates is still a widely open topic of investigation. Because all but the simplest problems are notoriously difficult to handle analytically, most approaches assume one of the two extreme limits: either one assumes a quasi-adiabatic (infinitely slow) cooling, or an instantaneous (infinitely fast) quench, often even in numerical studies. However, in realistic applications all cooling occurs with a finite rate and, therefore, lies somewhere between these two extremes. Some notable exceptions are Wong and Knobler’s [1] and Binder’s [2] work on quenches within the one-phase region, Wong and Knobler’s [3] and Ruiz’s [4] work on double quenches, and Onuki’s [5] work on periodic quenches. For reviews, see [6,7,9]. More recently, Laguna and Zurek [10], and Bettencourt et al. [11] have investigated how different cooling rates affect the formation of topological defects in classical fields, while Bowick and Momen [12] have examined this issue for quantum fields. ∗
Corresponding author. Present address: Department of Physics, Clark University, Worcester, MA 01610, USA. E-mail addresses: [email protected] (C.J. Gagne), [email protected] (M. Gleiser). 0167-2789/03/$ – see front matter © 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0167-2789(03)00092-7
122
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
Although these works have done much to advance our understanding of this issue, the effect of different cooling rates on the dynamics of phase transitions, as well as other nonequilibrium situations that commonly arise in nonlinear field theories, requires much further study. Even if advances in computer technology in recent years have greatly facilitated numerical studies of phase transitions, some aspects of nonequilibrium dynamics of nonlinear field theories remain difficult to study numerically, as they necessarily require large lattices and/or long simulations. In particular, dynamical studies near the critical region of continuous phase transitions, or of nucleation in discontinuous phase transitions, require both long runs and large lattices. Once we focus on dynamical issues, it is often desirable to observe how a system evolves towards its final equilibrium state. It is thus important to develop tools designed to quantify the nonequilibrium behavior of a system being cooled (or warmed up) in a way that is numerically efficient. With this goal in mind, in this paper we propose a possible measure to quantify the departure from equilibrium of a system coupled to a heat bath, which clearly correlates the approach to equilibrium with the relevant control parameters of the system, namely, the absolute value of the rate of change of the momentum-integrated structure function, (t), which will be defined below. In what follows, we will apply this measure to two relevant situations: a system undergoing a temperature quench, and a system undergoing a pressure quench, while the temperature is held constant. The paper is organized as follows. In Section 2 we introduce the dynamical equation to be used throughout this work, a generalized Langevin equation, and define the measure of the departure from equilibrium, (t). In Section 3, we apply this measure to a system with a quadratic time-independent potential, undergoing a temperature quench. In Section 4, we investigate a system with a time-dependent Ginzburg–Landau potential undergoing a pressure quench, with fixed temperature. In Section 5 we describe our simulation results for this system, illustrating the use of our measure to quantify the departure from equilibrium. We conclude in Section 6 summarizing our results, and pointing out possible future avenues of research.
2. Defining (t) In general, the model is characterized by a potential V(φ, t) (or, equivalently, the homogeneous part of the free-energy density for the order parameter φ) coupled in the simulations to a thermal bath using a generalized Langevin equation ∂2 φ ∂2 φ ∂φ ∂V − + ξ(x, t), = −η 2 2 ∂t ∂φ ∂t ∂x
(1)
where the temperature of the bath T is related to the damping rate η and the stochastic force of zero mean ξ(x, t) by the fluctuation–dissipation relation ξ(x, t)ξ(x , t ) = 2ηTδ(x − x )δ(t − t ).
(2)
We use a single classical real scalar field φ(x, t) in (1 + 1) dimensions. Although there is no phase transition in one dimension, we can still examine (local) symmetry breaking and phase ordering through the formation of kink–antikink pairs. The present work should be considered a further step in the quantification of nonequilibrium behavior, which can be extended to higher-dimensional systems. All quantities have been scaled appropriately to be dimensionless and periodic boundary conditions are used. At any given time during the evolution of the system, there will be local fluctuations around the spatially-averaged ¯ ¯ order parameter, φ(t). Writing these fluctuations as u(x, t) = φ(x, t) − φ(t), and their Fourier transform for a given wave number k as uk (t), we define the structure function as Sk (t) = |uk (t)|2 . Clearly, Sk (t) describes how individual Fourier modes evolve in time [7]. The derivative ∂Sk (t)/∂t gives information about the rate of change of the
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
123
individual modes. (Note that Sk (t) is not the dynamical structure function.) The integral of ∂Sk (t)/∂t over wave number k gives the net change of the fluctuations in the order parameter at a given time. Of course, one could equally well investigate the evolution of individual momentum modes. However, we are interested in obtaining a global property of the system, independent of the particular equilibration scale of a given mode. We thus define the quantity (t) as 2π/δx ∂Sk (t) (t) ≡ dk (3) . 2π/L ∂t δx is the lattice spacing and L the lattice size, defining the ultraviolet and infrared cutoffs, respectively. The integral is over the first Brillouin zone, and is logarithmically divergent with the UV cutoff, of O ∼ 1/δx. If needed (not the case here, as discussed further below), one could control this divergence by adding a lattice counterterm to the potential used in the numerical simulation, as was done in Ref. [8]. It should be noted that the time derivative makes (t) odd in time. Therefore, ignoring the absolute value for a moment, (t) seems a good candidate for a quantity that will tend to zero in equilibrium. Since we are only interested in departure, or “distance”, essentially, from equilibrium, the overall sign is not important and we have taken the absolute value for convenience and clarity of presentation. As we will see in Sections 3 and 4, (t) will display a peak at a certain t which can be used to characterize the nonequilibrium dynamics of the system. In equilibrium, (t) 0. Nonvanishing values of (t) can be used to measure the departure from equilibrium of the system. If needed, (t) may be smoothed by sampling every few time steps, or, equivalently, by averaging over the same number of time steps. The analytical calculation in Section 3, and the corresponding numerical simulations, confirm the sensitivity of (t) to the departure from equilibrium, as expected. In Section 3, we investigate a simple model, where the system, prepared initially in equilibrium at temperature T0 , is instantaneously quenched to a final temperature T = 0. If, for example, the initial temperature is doubled, (t) should reflect that. We will show, both analytically and numerically, that this is indeed the case. Notice that this is not a sensitivity to initial conditions, since the system starts in thermal equilibrium (only the equilibration temperature changes) but rather a sensitivity that we want, measuring the departure from equilibrium, which in this example is induced by an instantaneous quench. Also, as noted above, due to the logarithmic UV dependence, the absolute height of the peak is sensitive to the choice of lattice spacing. However, we have confirmed that the qualitative behavior of the peak height and its location in time as the damping coefficient or cooling rates are varied, do not depend on lattice spacing. Since our goal here is simply to demonstrate the usefulness of our measure, we did not add the lattice counterterm to our simulations. This could easily be done, following Ref. [8]. Finally, it should be noted that stationary states (such as metastable states) would also have a (t) ≈ 0, for the lifetime of that state. Thus, if one computes (t) for time-scales longer than the stationary state’s lifetime, our method can actually be used to identify the presence of these states in the system, while also providing a measure of their lifetimes.
3. (t) for an instantaneous temperature quench in a time-independent quadratic potential As a first illustration of the usefulness of (t) to quantify the departure from equilibrium, we will investigate an instantaneous temperature quench in a simple quadratic potential. We chose this situation not for its experimental usefulness, but because it shows departure and return to equilibrium, while allowing us to calculate (t) analytically and compare the results to the numerical simulation.
124
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
3.1. (t) for an instantaneous quench to T = 0: analytical computation The field φ(x, t) is initially thermalized in a time-independent potential V(φ) = (1/2)m2 φ2 at a temperature T0 . At t = 0 the system is instantaneously quenched to the temperature T = 0. (Any choice of final temperature here would do.) At t ≥ 0, the equation of motion is (in units where c = kB = 1) ∂2 φ(x, t) ∂φ(x, t) ∂2 φ(x, t) = −η − m2 φ(x, t). ∂t ∂t 2 ∂x2
(4)
Since the average value of φ(x, t) is zero, this is also the equation for the field fluctuations. Fourier transforming φ(x, t) gives an equation of motion for φ(k, t), ∂φ(k, t) ∂2 φ(k, t) = −(k2 + m2 )φ(k, t) − η , 2 ∂t ∂t
(5)
which has the solution (6) φ(k, t) = e−ηt/2 [A(k) cos (ωt) + B(k) sin (ωt)], where ω = k2 + m2 − η2 /4 is the dispersion relation in the presence of viscosity, and A(k) and B(k) are to be determined. The exponentially growing solution has been dropped because of the boundary condition for φ(x, t) and for φ(k, t) → 0 at t → ∞. The structure function is defined as S(k, t) = |φ(k, t)|2 ,
(7)
where the angular brackets refer to a statistical average. S(k, t) can be written as S(k, t) = e−ηt |A(k)|2 cos 2 (ωt) + |B(k)|2 sin 2 (ωt) + [A(k)B (k) + B(k)A (k)] cos (ωt) sin (ωt) .
(8)
|A(k)|2 is obtained by the condition that at t = 0 the structure function can be approximated by the Ornstein–Zernicke form, so that S(k, t = 0) = |A(k)|2 =
k2
T0 . + m2
(9)
Taking the time derivative of S(k, t), ∂S(k, t) = e−ηt {−η|A(k)|2 + ω[A(k)B (k) + B(k)A (k)]} cos 2 (ωt) ∂t + {−η|B(k)|2 − ω[A(k)B (k) + B(k)A (k)]} sin 2 (ωt) + {−η[A(k)B (k) + B(k)A (k)] + 2ω[−|A(k)|2 + |B(k)|2 ]} cos (ωt) sin (ωt) .
(10)
The structure function should be stationary at t = 0, before the quench, so that ∂S(k, t = 0) = −η|A(k)|2 + ωA(k)B (k) + B(k)A (k) = C, ∂t
(11)
where C is a constant. C is added for convenience, so that we can match the analytical results with the numerical simulation. This follows from the use of a discretized time in the simulation, which precludes a precise match at t = 0. Clearly, the exact analytical results would have C = 0. And, in fact, the values obtained for the numerical match have C 1. (As δt → 0, we expect C → 0 as well.) This gives A(k)B (k) + B(k)A (k) =
C + η|A(k)|2 , ω
(12)
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
125
leaving only |B(k)|2 undetermined. If we assume A(k) and B(k) are real and statistically independent, |B(k)|2 = and
C2 (k2 + m2 ) η2 T0 ηC + , + 2 2 4ω T0 2ω (k2 + m2 ) 2ω2
2 2 ηC (k + m2 ) ηT0 η2 C ∂S(k, t) −ηt 2 =e C cos (ωt) − + 2 + + C sin 2 (ωt) ∂t 4ω2 T0 ω 2ω2 2 2 C (k + m2 ) T0 + − sin (2ωt) . 4ωT0 ω
The quantity we are after is 2π/δx ∂S(k, t) (t) = dk , 2π/L ∂t
(13)
(14)
(15)
where δx is the lattice spacing and L the lattice size. The integral is over the first Brillouin zone, omitting k = 0, because it is zero by convention, since we areconcerned with fluctuations from the average. If we change the integration variable to y = ω/ω0 , where ω0 = m2 − η2 /4 is the k = 0 dispersion relation, (t) becomes y cos 2 (ω0 ty) ηC y sin 2 (ω0 ty) −ηt (t) = e ω0 C dy − ω0 C 1 + dy 4T0 y2 − 1 y2 − 1
ηT0 sin 2 (ω0 ty) η2 C 2 sin (2ω0 ty) ηC 2 − 1+ dy − T0 1 − dy 2 2 ω0 4T0 16T0 y y −1 y2 − 1 2 ω02 C2 y sin 2 (ω0 ty) + dy . (16) 4T0 y2 − 1 The integrals can be evaluated analytically for η ≤ 2, if we take the ideal limits δx → 0 and L → ∞, so that the integration limits become 1 ≤ y ≤ ∞: 3 ηC 2 π 1 π ηC π ηT0 −ηt 2 2 (t) e ω0 C J1 (2ω0 t) − ω0 C 1 + J1 (2ω0 t) − ω0 t1 F 2 ; 1, ; −a t 1+ 4 4T0 4 ω0 4T0 2 2 2
ω02 C2 π J1 (2ω0 t) ω02 C2 π η2 C 2 π − T0 1 − J (2ω Im(ω0 t) = 0. t) + − J (2ω t) 0 0 2 0 sign(ω0 t), 4T0 2 2ω0 t 4T0 2 16T02 2 (17) Since numerically we sum over k, rather than integrate, (t) above must be multiplied by an additional k−1 = L/2π. After adding a normalization, N, defining the combination D ≡ C/4T0 , and scaling out m to make the equation dimensionless, (t) becomes sign(ω0 t)πNT0 L e−ηt 1 − D2 D J (t) (2ω t) + ω D 2 + Dη − J1 (2ω0 t) + 2ω02 D2 J2 (2ω0 t) 0 0 0 2 2 t ηt 1 3 + (1 + D)2 1 F2 (18) ; 1, ; −ω02 t 2 , 2 2 2 where J0 , J1 and J2 are Bessel functions of the first kind and 1 F2 is a hypergeometric function.
126
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
Fig. 1. (t) and the average kinetic energy K(t) for an instantaneous heat bath temperature quench from T0 = 0 to T = 0 within a single phase in the time-independent potential.
3.2. (t) for an instantaneous quench to T = 0: numerical results Fig. 1 shows (t) and the average kinetic energy density, K(t), for a field in a time-independent potential, V(φ) = (1/2)φ2 . The field φ is initially thermalized to a temperature T0 with fixed damping rate η = 1. We considered two different values for T0 , T0 = 0.01 and T0 = 0.005. At t = 0, the bath is instantaneously quenched to zero temperature. Notice how initially (t) is approximately zero because the field is thermalized (K(t) = T0 /2), but begins immediately to rise after the quench, reaching a maximum and returning to zero after the field thermalizes (K(t) → 0). Since the field is thermalized at t = 0, it is possible to compute the initial value of the structure function, Sk (0). As the system evolves deterministically after the quench, we can use this initial value to compute (t) for t > 0 using Eq. (18). N and D are parameters that were varied to fit the data. In Fig. 1 we show numerical (circles) and analytical (lines) results for both initial temperatures, T0 = 0.01 and 0.005. The simulation data were fit using N = 0.82 and D = 0.65 for both T0 = 0.01 (filled circles) and T0 = 0.005 (open circles). Note that this value of D is consistent with small C, as it should. We have verified (not shown) that both the time, tpeak , and amplitude of the peak, Apeak , depend on the damping rate. Also, the peak amplitude (Apeak ) increases with temperature change, as can be seen in Fig. 1; thus, the location of the peak gives a measure of the equilibration time-scale of the system, while its amplitude provides a measure of the departure from equilibrium. These properties will be further explored below.
4. (t) for a self-interacting potential for varying quenching rates at constant temperature For simplicity of comparison, we chose to use Laguna and Zurek’s model of a linear pressure quench for a single classical real scalar field φ(x, t) in (1 + 1) dimensions [10]. Although there is no phase transition in one dimension, we can still examine (local) symmetry breaking and phase ordering through the formation of kink–antikink pairs. The present work should be considered a further step in the quantification of nonequilibrium behavior, which can be extended to higher-dimensional systems. All quantities have been scaled appropriately to be dimensionless.
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
127
The linear cooling model is characterized by the following Ginzburg–Landau type potential (or, equivalently, the homogeneous part of the free-energy density for the order parameter φ), V(φ) = 18 (1 − 2'(t)φ2 + φ4 ),
(19)
where '(t) ≡ min(t/τq , 1), and the time t is set to zero at the “phase transition point” (when V (φ = 0) = 0). When t < 0, '(t) < 0 and V(φ) is a single well potential. When t > 0, '(t) > 0 and V(φ) is a double well potential. τq is the quenching time-scale. Note that the potential stops changing when '(t) = 1, which, in a Ginzburg–Landau model, is the zero-temperature limit. The terminology “pressure quench” is justified by the fact that the temperature of the bath T is kept constant, while the quadratic coefficient of the potential '(t) changes linearly in time. As argued in Ref. [10], this is equivalent to a pressure quench at constant temperature in the laboratory. An obvious extension of this work is to implement a true linear cooling, where the temperature of the bath changes, and '(t) = (Tc − T)/Tc . We will leave this case for a future investigation (see, however, Ref. [11]). Our main interest here is in proposing a measure of departure from equilibrium of classical fields, which can be adapted for several different situations, including those involving cooling through the bath temperature, as we show below. Laguna and Zurek [10] have studied the effect of different quenching rates on the density of zero crossings of the field φ. (The density is given by the number of times φ goes through zero divided by the lattice length and provides the approximate kink density [10].) In general, φ will attempt to keep up with the changing potential as best it can, that is, its Fourier modes will try to remain thermalized. As the quenching rate is increased, we can envisage a situation where the various modes of φ will not be able to maintain thermal equilibrium with the bath and we say that φ becomes “frozen”. According to Zurek’s conjecture [10], this freeze-out occurs when the ˙ for overdamped systems, and by τφ¨ |φ/φ| ¨ 1/2 for underdamped dynamical relaxation time—given by τφ˙ |φ/φ| systems—is comparable with the time to (from) the phase transition. From this freeze-out condition, one can derive 1/3 1/2 scaling relationships for the freeze-out time, τˆφ˙ ∝ τq and τˆφ¨ ∝ τq . Using these results, Zurek and Laguna find scaling laws for kink density, which they confirm with simulations and contrast with experimental results for pressure quenches. Here, we are mostly concerned with how to extend our knowledge of nonequilibrium properties of field theories. We will thus be using Zurek and Laguna’s model as a testing ground for our methods, comparing some of our results to theirs. Even though their freeze-out time is clearly not the same as our peak time, since the former measures when the system first falls out of equilibrium while the latter measures when the system is maximally out of equilibrium, both time variables scale similarly with cooling time-scale.
5. Simulation results for time-dependent self-interacting potential at constant temperature In addition to testing our second-order staggered leapfrog Langevin code for stability with the same parameters as [10] (L = 2048, δx = 0.125, δt = 0.025 and T = 0.01), we tested the reliability of our code by reproducing their defect density results. For a more thorough description of the algorithm used, see [8]. The calculation of the structure function uses an FFT routine from Numerical Recipes [13], while the change in the structure function uses a simple finite difference method, with the ability to average over a few time steps to smooth the data. We can then extract the amplitude and time of the peak. Fig. 2 shows the variation in (t/τq ) with damping rate and cooling time-scale. t/τq is used on the time axis for easier comparison between different cooling time-scales. The amplitude of the peak is higher the lower the damping rate and the shorter the cooling time-scale (faster cooling). Also, (t/τq ) peaks earlier for larger cooling time-scales (slower cooling). The simulations start with the field φ thermalized at t/τq = −2, so that (t) is initially zero. As time progresses the potential changes and, if the thermalization rate of the field is not faster than the potential’s quenching rate, (t)
128
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
Fig. 2. (t) for different damping rates η and cooling time-scales τq with the potential in Eq. (1). The amplitude of the peak, Apeak , increases with decreasing time-scale and decreases with increasing damping rate. tpeak /τq increases with decreasing cooling time-scale and increasing damping rate.
begins to change, so that, even before the phase transition, it is nonzero and increasing. The most drastic changes in (t) happen as the potential goes from a single well to a double well, as we would expect: (t) peaks some time after the critical point. The peaks occurring during the cooling show a plateau beginning after the peak and lasting until the cooling ends at ' = 1. This can be seen in Fig. 2 for τq = 48, represented by the dashed line. For η = 1,
Fig. 3. Peak time vs. cooling time-scale τq for time-dependent potential for several damping rates. The solid straight lines above and below the data show the scaling predicted in Ref. [9] for overdamped and underdamped freeze-out times as a function of cooling time-scale. The simulation data interpolates between these two extremes appropriately.
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
129
Fig. 4. Amplitude of peak vs. cooling time-scale for time-dependent potential with varying damping rate. For all but the fastest quench times, −6/5 all damping rates studied satisfy the scaling Apeak ∝ τq . All three fit lines have a slope of −6/5.
we were able to identify a peak up to τq ∼ 128 at which point the data become too noisy to distinguish the peak from the plateau. For the fastest coolings (which, with η = 1, happen for τq 16), (t/τq ) actually peaks after the potential stops changing (' = t/τq = 1): the system remains out of equilibrium during the whole cooling process, which is thus equivalent to an instantaneous quench. For τq 256, (t) 0, and the system remains thermalized during the cooling; this is equivalent to an adiabatic cooling. It is between these two regimes, which, for η = 1, are given by 16 τq 256, that the cooling most dramatically affects the dynamics of the system. The location (tpeak ) and amplitude (Apeak ) of the peak in (t) are plotted in Figs. 3 and 4, respectively, as a function of τq for several damping rates. From Fig. 3 we note that the location of the peak scales with τq , in agreement with the scaling obtained in Ref. [10] for the freeze-out time as a function of τq . This is consistent with the fact that tpeak is the time after which the field is able to begin relaxing to its equilibrium state. −p We found that the amplitude of the peak, Apeak , scales as τq , where p = 1.2 ± 0.05 for all but the fastest coolings. We conclude that there are essentially three regimes: (i) the field remains thermalized all through the cooling process, thus approximating an “adiabatic” cooling (for τq 256, with η = 1); (ii) the field remains out of equilibrium throughout the duration of the cooling, thus approximating an instantaneous quench (for τq 16 and η = 1); (iii) the intermediate regime, where the field reaches thermalization during the cooling process, signaled by the appearance of a plateau in (t), which decreases in amplitude as τq is increased, until it cannot be differentiated from noise (for τq 256).
6. Summary We have presented a novel method to measure the departure from thermal equilibrium in nonlinear field theories, based on the computation of the momentum-integrated time-derivative of the structure function, which we called (t). In order to illustrate our method, we have discussed two particular cases: in the first one, we followed the evolution of a field in a quadratic potential prepared initially in a thermal state at T0 and quenched instantaneously
130
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
to T = 0. Due to the simplicity of the model, we were able to produce an analytical solution which was compared with the results from the numerical simulation. In the second case, we studied the Ginzburg–Landau potential with a linear pressure quench as suggested by Laguna and Zurek [10]. After reproducing their results, we showed how nonequilibrium is characterized by the presence of a peak in (t). The amplitude of the peak is a measure of the departure from equilibrium, while its location is related to the equilibration time-scale of the system as compared to the quenching time-scale. The general picture can be described qualitatively as follows: the modes are initially thermalized, and their rate of change increases as they try to “catch up” with the changing potential (or, when appropriate, with the changing heat bath). This may or may not happen during the cooling, determined by the value of τq . As the quench proceeds at the steady rate of τq−1 , (t) increases until the field begins to probe the bottom of the potential (the free-energy minima), when it has essentially “caught up” with the changing potential. If this happens while the potential is still changing, the field then changes at the same rate as the potential ((t) begins to plateau), until the potential stops changing and the field can fully thermalize ((t) → 0). For the fastest coolings (τq 16 for η = 1), the field is not able to “catch up” with the changing potential before the end of the cooling, and so (t) never plateaus; the amplitude of the peak does not obey the scaling in Fig. 4 because the field is still changing after the cooling has ended, which, as we remarked above, approximates an instantaneous quench. Finally, we note that our approach is applicable to any scalar field theory in the Ising universality class, with or without topological defects. It would be interesting to apply it to systems with metastable (stationary) states, as a tool to measure the lifetime and to quantify the approach to equilibrium in those systems. Also, it should be straightforward to obtain an expression of (t) for systems described by a complex scalar order parameter, such as 4 He. The extension of our approach to systems with local gauge invariance, applicable to the description of superconductors, remains open for investigation.
Acknowledgements CG was supported in part by a National Science Foundation grant no. PHY-9453431 and Dartmouth College. MG was supported in part by NSF grants PHY-0070554 and PHY-9453431, and by the Mr. Tomkins Fund for Cosmology and Field Theory at Dartmouth College. CG thanks the Theoretical Condensed Matter Group at Boston University, where parts of this work were developed, for their hospitality. We both thank Sidney Redner, Paul Krapivsky, Wonmuk Hwang, Prashant Sharma, Victor Spirin, Harvey Gould and Bill Klein for many helpful discussions. References [1] N.C. Wong, C.M. Knobler, Phys. Rev. Lett. 43 (1979) 1733; N.C. Wong, C.M. Knobler, Phys. Rev. Lett. 45 (1979) 498. [2] K. Binder, Phys. Rev. B 15 (1977) 4425. [3] N.C. Wong, C.M. Knobler, J. Chem. Phys. 69 (1978) 725. [4] R. Ruiz, Phys. Rev. A 26 (1982) 2227. [5] A. Onuki, Prog. Theoret. Phys. 66 (1981) 1230; A. Onuki, Prog. Theoret. Phys. 67 (1982) 768; A. Onuki, Prog. Theoret. Phys. 67 (1982) 787; A. Onuki, Phys. Rev. Lett. 48 (1982) 753. [6] J.D. Gunton, M. San Miguel, P.S. Sahni, in: C. Domb, J.L. Lebowitz (Eds.), Phase Transitions and Critical Phenomena, vol. 8, Academic Press, London, 1983. [7] J.S. Langer, in: C. Godreche (Ed.), Solids Far from Equilibrium, Cambridge University Press, Cambridge, 1992. [8] C.J. Gagne, M. Gleiser, Phys. Rev. E 61 (2000) 3483; M. Gleiser, H. Müller, Phys. Lett. B 422 (1998) 69.
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
131
[9] A.J. Bray, Adv. Phys. 43 (1994) 357. [10] P. Laguna, W.H. Zurek, Phys. Rev. D 58 (1998) 085021; P. Laguna, W.H. Zurek, Phys. Rev. Lett. 78 (1997) 2519; W.H. Zurek, Phys. Rep. 276 (1996) 177; W.H. Zurek, Acta Pol. B 24 (1993) 1301; W.H. Zurek, Nature (London) 317 (1985) 505; Los Alamos Report No. LAUR 84-3818, Unpublished, 1984. [11] L.M.A. Bettencourt, N.D. Antunes, W.H. Zurek, Phys. Rev. D 62 (2000) 065005. [12] M. Bowick, A. Momen, Phys. Rev. D 58 (1998) 085014. [13] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed., Cambridge University Press, Cambridge, 1992.
Quantifying nonequilibrium behavior with varying quenching rates Carmen J. Gagne∗ , Marcelo Gleiser Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA Received 7 May 2002; received in revised form 29 January 2003; accepted 27 February 2003 Communicated by L. Kramer
Abstract We investigate the approach to thermal equilibrium of (1 + 1)-dimensional stochastic Ginzburg–Landau models at varying cooling rates. The nonequilibrium dynamics is modeled by coupling the field to an external heat bath with damping rate η. We argue that the departure from thermal equilibrium can be measured from the absolute value of the rate of change of the momentum-integrated structure function, (t). If the field equilibration time-scale is faster than the cooling (or quench) time-scale τq , then (t) → 0. Otherwise, (t) displays a peak which scales as τqn , with 1/3 ≤ n ≤ 1/2 as η is varied from the underdamped to the overdamped limit. Furthermore, we show that the amplitude of the peak in (t) measures the departure −6/5 from equilibrium and scales as τq , independent of η. © 2003 Elsevier Science B.V. All rights reserved. PACS: 05.70.Ln; 11.10.Wx; 98.80.Cq Keywords: Nonequilibrium dynamics; Stochastic PDEs; Ginzburg–Landau models; Finite quench rate
1. Introduction Given the widespread role of phase transitions in condensed matter physics, cosmology, and high energy physics, it is imperative that the nonequilibrium dynamics of these phase transitions be better understood. In particular, the effect of varying cooling rates is still a widely open topic of investigation. Because all but the simplest problems are notoriously difficult to handle analytically, most approaches assume one of the two extreme limits: either one assumes a quasi-adiabatic (infinitely slow) cooling, or an instantaneous (infinitely fast) quench, often even in numerical studies. However, in realistic applications all cooling occurs with a finite rate and, therefore, lies somewhere between these two extremes. Some notable exceptions are Wong and Knobler’s [1] and Binder’s [2] work on quenches within the one-phase region, Wong and Knobler’s [3] and Ruiz’s [4] work on double quenches, and Onuki’s [5] work on periodic quenches. For reviews, see [6,7,9]. More recently, Laguna and Zurek [10], and Bettencourt et al. [11] have investigated how different cooling rates affect the formation of topological defects in classical fields, while Bowick and Momen [12] have examined this issue for quantum fields. ∗
Corresponding author. Present address: Department of Physics, Clark University, Worcester, MA 01610, USA. E-mail addresses: [email protected] (C.J. Gagne), [email protected] (M. Gleiser). 0167-2789/03/$ – see front matter © 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0167-2789(03)00092-7
122
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
Although these works have done much to advance our understanding of this issue, the effect of different cooling rates on the dynamics of phase transitions, as well as other nonequilibrium situations that commonly arise in nonlinear field theories, requires much further study. Even if advances in computer technology in recent years have greatly facilitated numerical studies of phase transitions, some aspects of nonequilibrium dynamics of nonlinear field theories remain difficult to study numerically, as they necessarily require large lattices and/or long simulations. In particular, dynamical studies near the critical region of continuous phase transitions, or of nucleation in discontinuous phase transitions, require both long runs and large lattices. Once we focus on dynamical issues, it is often desirable to observe how a system evolves towards its final equilibrium state. It is thus important to develop tools designed to quantify the nonequilibrium behavior of a system being cooled (or warmed up) in a way that is numerically efficient. With this goal in mind, in this paper we propose a possible measure to quantify the departure from equilibrium of a system coupled to a heat bath, which clearly correlates the approach to equilibrium with the relevant control parameters of the system, namely, the absolute value of the rate of change of the momentum-integrated structure function, (t), which will be defined below. In what follows, we will apply this measure to two relevant situations: a system undergoing a temperature quench, and a system undergoing a pressure quench, while the temperature is held constant. The paper is organized as follows. In Section 2 we introduce the dynamical equation to be used throughout this work, a generalized Langevin equation, and define the measure of the departure from equilibrium, (t). In Section 3, we apply this measure to a system with a quadratic time-independent potential, undergoing a temperature quench. In Section 4, we investigate a system with a time-dependent Ginzburg–Landau potential undergoing a pressure quench, with fixed temperature. In Section 5 we describe our simulation results for this system, illustrating the use of our measure to quantify the departure from equilibrium. We conclude in Section 6 summarizing our results, and pointing out possible future avenues of research.
2. Defining (t) In general, the model is characterized by a potential V(φ, t) (or, equivalently, the homogeneous part of the free-energy density for the order parameter φ) coupled in the simulations to a thermal bath using a generalized Langevin equation ∂2 φ ∂2 φ ∂φ ∂V − + ξ(x, t), = −η 2 2 ∂t ∂φ ∂t ∂x
(1)
where the temperature of the bath T is related to the damping rate η and the stochastic force of zero mean ξ(x, t) by the fluctuation–dissipation relation ξ(x, t)ξ(x , t ) = 2ηTδ(x − x )δ(t − t ).
(2)
We use a single classical real scalar field φ(x, t) in (1 + 1) dimensions. Although there is no phase transition in one dimension, we can still examine (local) symmetry breaking and phase ordering through the formation of kink–antikink pairs. The present work should be considered a further step in the quantification of nonequilibrium behavior, which can be extended to higher-dimensional systems. All quantities have been scaled appropriately to be dimensionless and periodic boundary conditions are used. At any given time during the evolution of the system, there will be local fluctuations around the spatially-averaged ¯ ¯ order parameter, φ(t). Writing these fluctuations as u(x, t) = φ(x, t) − φ(t), and their Fourier transform for a given wave number k as uk (t), we define the structure function as Sk (t) = |uk (t)|2 . Clearly, Sk (t) describes how individual Fourier modes evolve in time [7]. The derivative ∂Sk (t)/∂t gives information about the rate of change of the
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
123
individual modes. (Note that Sk (t) is not the dynamical structure function.) The integral of ∂Sk (t)/∂t over wave number k gives the net change of the fluctuations in the order parameter at a given time. Of course, one could equally well investigate the evolution of individual momentum modes. However, we are interested in obtaining a global property of the system, independent of the particular equilibration scale of a given mode. We thus define the quantity (t) as 2π/δx ∂Sk (t) (t) ≡ dk (3) . 2π/L ∂t δx is the lattice spacing and L the lattice size, defining the ultraviolet and infrared cutoffs, respectively. The integral is over the first Brillouin zone, and is logarithmically divergent with the UV cutoff, of O ∼ 1/δx. If needed (not the case here, as discussed further below), one could control this divergence by adding a lattice counterterm to the potential used in the numerical simulation, as was done in Ref. [8]. It should be noted that the time derivative makes (t) odd in time. Therefore, ignoring the absolute value for a moment, (t) seems a good candidate for a quantity that will tend to zero in equilibrium. Since we are only interested in departure, or “distance”, essentially, from equilibrium, the overall sign is not important and we have taken the absolute value for convenience and clarity of presentation. As we will see in Sections 3 and 4, (t) will display a peak at a certain t which can be used to characterize the nonequilibrium dynamics of the system. In equilibrium, (t) 0. Nonvanishing values of (t) can be used to measure the departure from equilibrium of the system. If needed, (t) may be smoothed by sampling every few time steps, or, equivalently, by averaging over the same number of time steps. The analytical calculation in Section 3, and the corresponding numerical simulations, confirm the sensitivity of (t) to the departure from equilibrium, as expected. In Section 3, we investigate a simple model, where the system, prepared initially in equilibrium at temperature T0 , is instantaneously quenched to a final temperature T = 0. If, for example, the initial temperature is doubled, (t) should reflect that. We will show, both analytically and numerically, that this is indeed the case. Notice that this is not a sensitivity to initial conditions, since the system starts in thermal equilibrium (only the equilibration temperature changes) but rather a sensitivity that we want, measuring the departure from equilibrium, which in this example is induced by an instantaneous quench. Also, as noted above, due to the logarithmic UV dependence, the absolute height of the peak is sensitive to the choice of lattice spacing. However, we have confirmed that the qualitative behavior of the peak height and its location in time as the damping coefficient or cooling rates are varied, do not depend on lattice spacing. Since our goal here is simply to demonstrate the usefulness of our measure, we did not add the lattice counterterm to our simulations. This could easily be done, following Ref. [8]. Finally, it should be noted that stationary states (such as metastable states) would also have a (t) ≈ 0, for the lifetime of that state. Thus, if one computes (t) for time-scales longer than the stationary state’s lifetime, our method can actually be used to identify the presence of these states in the system, while also providing a measure of their lifetimes.
3. (t) for an instantaneous temperature quench in a time-independent quadratic potential As a first illustration of the usefulness of (t) to quantify the departure from equilibrium, we will investigate an instantaneous temperature quench in a simple quadratic potential. We chose this situation not for its experimental usefulness, but because it shows departure and return to equilibrium, while allowing us to calculate (t) analytically and compare the results to the numerical simulation.
124
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
3.1. (t) for an instantaneous quench to T = 0: analytical computation The field φ(x, t) is initially thermalized in a time-independent potential V(φ) = (1/2)m2 φ2 at a temperature T0 . At t = 0 the system is instantaneously quenched to the temperature T = 0. (Any choice of final temperature here would do.) At t ≥ 0, the equation of motion is (in units where c = kB = 1) ∂2 φ(x, t) ∂φ(x, t) ∂2 φ(x, t) = −η − m2 φ(x, t). ∂t ∂t 2 ∂x2
(4)
Since the average value of φ(x, t) is zero, this is also the equation for the field fluctuations. Fourier transforming φ(x, t) gives an equation of motion for φ(k, t), ∂φ(k, t) ∂2 φ(k, t) = −(k2 + m2 )φ(k, t) − η , 2 ∂t ∂t
(5)
which has the solution (6) φ(k, t) = e−ηt/2 [A(k) cos (ωt) + B(k) sin (ωt)], where ω = k2 + m2 − η2 /4 is the dispersion relation in the presence of viscosity, and A(k) and B(k) are to be determined. The exponentially growing solution has been dropped because of the boundary condition for φ(x, t) and for φ(k, t) → 0 at t → ∞. The structure function is defined as S(k, t) = |φ(k, t)|2 ,
(7)
where the angular brackets refer to a statistical average. S(k, t) can be written as S(k, t) = e−ηt |A(k)|2 cos 2 (ωt) + |B(k)|2 sin 2 (ωt) + [A(k)B (k) + B(k)A (k)] cos (ωt) sin (ωt) .
(8)
|A(k)|2 is obtained by the condition that at t = 0 the structure function can be approximated by the Ornstein–Zernicke form, so that S(k, t = 0) = |A(k)|2 =
k2
T0 . + m2
(9)
Taking the time derivative of S(k, t), ∂S(k, t) = e−ηt {−η|A(k)|2 + ω[A(k)B (k) + B(k)A (k)]} cos 2 (ωt) ∂t + {−η|B(k)|2 − ω[A(k)B (k) + B(k)A (k)]} sin 2 (ωt) + {−η[A(k)B (k) + B(k)A (k)] + 2ω[−|A(k)|2 + |B(k)|2 ]} cos (ωt) sin (ωt) .
(10)
The structure function should be stationary at t = 0, before the quench, so that ∂S(k, t = 0) = −η|A(k)|2 + ωA(k)B (k) + B(k)A (k) = C, ∂t
(11)
where C is a constant. C is added for convenience, so that we can match the analytical results with the numerical simulation. This follows from the use of a discretized time in the simulation, which precludes a precise match at t = 0. Clearly, the exact analytical results would have C = 0. And, in fact, the values obtained for the numerical match have C 1. (As δt → 0, we expect C → 0 as well.) This gives A(k)B (k) + B(k)A (k) =
C + η|A(k)|2 , ω
(12)
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
125
leaving only |B(k)|2 undetermined. If we assume A(k) and B(k) are real and statistically independent, |B(k)|2 = and
C2 (k2 + m2 ) η2 T0 ηC + , + 2 2 4ω T0 2ω (k2 + m2 ) 2ω2
2 2 ηC (k + m2 ) ηT0 η2 C ∂S(k, t) −ηt 2 =e C cos (ωt) − + 2 + + C sin 2 (ωt) ∂t 4ω2 T0 ω 2ω2 2 2 C (k + m2 ) T0 + − sin (2ωt) . 4ωT0 ω
The quantity we are after is 2π/δx ∂S(k, t) (t) = dk , 2π/L ∂t
(13)
(14)
(15)
where δx is the lattice spacing and L the lattice size. The integral is over the first Brillouin zone, omitting k = 0, because it is zero by convention, since we areconcerned with fluctuations from the average. If we change the integration variable to y = ω/ω0 , where ω0 = m2 − η2 /4 is the k = 0 dispersion relation, (t) becomes y cos 2 (ω0 ty) ηC y sin 2 (ω0 ty) −ηt (t) = e ω0 C dy − ω0 C 1 + dy 4T0 y2 − 1 y2 − 1
ηT0 sin 2 (ω0 ty) η2 C 2 sin (2ω0 ty) ηC 2 − 1+ dy − T0 1 − dy 2 2 ω0 4T0 16T0 y y −1 y2 − 1 2 ω02 C2 y sin 2 (ω0 ty) + dy . (16) 4T0 y2 − 1 The integrals can be evaluated analytically for η ≤ 2, if we take the ideal limits δx → 0 and L → ∞, so that the integration limits become 1 ≤ y ≤ ∞: 3 ηC 2 π 1 π ηC π ηT0 −ηt 2 2 (t) e ω0 C J1 (2ω0 t) − ω0 C 1 + J1 (2ω0 t) − ω0 t1 F 2 ; 1, ; −a t 1+ 4 4T0 4 ω0 4T0 2 2 2
ω02 C2 π J1 (2ω0 t) ω02 C2 π η2 C 2 π − T0 1 − J (2ω Im(ω0 t) = 0. t) + − J (2ω t) 0 0 2 0 sign(ω0 t), 4T0 2 2ω0 t 4T0 2 16T02 2 (17) Since numerically we sum over k, rather than integrate, (t) above must be multiplied by an additional k−1 = L/2π. After adding a normalization, N, defining the combination D ≡ C/4T0 , and scaling out m to make the equation dimensionless, (t) becomes sign(ω0 t)πNT0 L e−ηt 1 − D2 D J (t) (2ω t) + ω D 2 + Dη − J1 (2ω0 t) + 2ω02 D2 J2 (2ω0 t) 0 0 0 2 2 t ηt 1 3 + (1 + D)2 1 F2 (18) ; 1, ; −ω02 t 2 , 2 2 2 where J0 , J1 and J2 are Bessel functions of the first kind and 1 F2 is a hypergeometric function.
126
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
Fig. 1. (t) and the average kinetic energy K(t) for an instantaneous heat bath temperature quench from T0 = 0 to T = 0 within a single phase in the time-independent potential.
3.2. (t) for an instantaneous quench to T = 0: numerical results Fig. 1 shows (t) and the average kinetic energy density, K(t), for a field in a time-independent potential, V(φ) = (1/2)φ2 . The field φ is initially thermalized to a temperature T0 with fixed damping rate η = 1. We considered two different values for T0 , T0 = 0.01 and T0 = 0.005. At t = 0, the bath is instantaneously quenched to zero temperature. Notice how initially (t) is approximately zero because the field is thermalized (K(t) = T0 /2), but begins immediately to rise after the quench, reaching a maximum and returning to zero after the field thermalizes (K(t) → 0). Since the field is thermalized at t = 0, it is possible to compute the initial value of the structure function, Sk (0). As the system evolves deterministically after the quench, we can use this initial value to compute (t) for t > 0 using Eq. (18). N and D are parameters that were varied to fit the data. In Fig. 1 we show numerical (circles) and analytical (lines) results for both initial temperatures, T0 = 0.01 and 0.005. The simulation data were fit using N = 0.82 and D = 0.65 for both T0 = 0.01 (filled circles) and T0 = 0.005 (open circles). Note that this value of D is consistent with small C, as it should. We have verified (not shown) that both the time, tpeak , and amplitude of the peak, Apeak , depend on the damping rate. Also, the peak amplitude (Apeak ) increases with temperature change, as can be seen in Fig. 1; thus, the location of the peak gives a measure of the equilibration time-scale of the system, while its amplitude provides a measure of the departure from equilibrium. These properties will be further explored below.
4. (t) for a self-interacting potential for varying quenching rates at constant temperature For simplicity of comparison, we chose to use Laguna and Zurek’s model of a linear pressure quench for a single classical real scalar field φ(x, t) in (1 + 1) dimensions [10]. Although there is no phase transition in one dimension, we can still examine (local) symmetry breaking and phase ordering through the formation of kink–antikink pairs. The present work should be considered a further step in the quantification of nonequilibrium behavior, which can be extended to higher-dimensional systems. All quantities have been scaled appropriately to be dimensionless.
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
127
The linear cooling model is characterized by the following Ginzburg–Landau type potential (or, equivalently, the homogeneous part of the free-energy density for the order parameter φ), V(φ) = 18 (1 − 2'(t)φ2 + φ4 ),
(19)
where '(t) ≡ min(t/τq , 1), and the time t is set to zero at the “phase transition point” (when V (φ = 0) = 0). When t < 0, '(t) < 0 and V(φ) is a single well potential. When t > 0, '(t) > 0 and V(φ) is a double well potential. τq is the quenching time-scale. Note that the potential stops changing when '(t) = 1, which, in a Ginzburg–Landau model, is the zero-temperature limit. The terminology “pressure quench” is justified by the fact that the temperature of the bath T is kept constant, while the quadratic coefficient of the potential '(t) changes linearly in time. As argued in Ref. [10], this is equivalent to a pressure quench at constant temperature in the laboratory. An obvious extension of this work is to implement a true linear cooling, where the temperature of the bath changes, and '(t) = (Tc − T)/Tc . We will leave this case for a future investigation (see, however, Ref. [11]). Our main interest here is in proposing a measure of departure from equilibrium of classical fields, which can be adapted for several different situations, including those involving cooling through the bath temperature, as we show below. Laguna and Zurek [10] have studied the effect of different quenching rates on the density of zero crossings of the field φ. (The density is given by the number of times φ goes through zero divided by the lattice length and provides the approximate kink density [10].) In general, φ will attempt to keep up with the changing potential as best it can, that is, its Fourier modes will try to remain thermalized. As the quenching rate is increased, we can envisage a situation where the various modes of φ will not be able to maintain thermal equilibrium with the bath and we say that φ becomes “frozen”. According to Zurek’s conjecture [10], this freeze-out occurs when the ˙ for overdamped systems, and by τφ¨ |φ/φ| ¨ 1/2 for underdamped dynamical relaxation time—given by τφ˙ |φ/φ| systems—is comparable with the time to (from) the phase transition. From this freeze-out condition, one can derive 1/3 1/2 scaling relationships for the freeze-out time, τˆφ˙ ∝ τq and τˆφ¨ ∝ τq . Using these results, Zurek and Laguna find scaling laws for kink density, which they confirm with simulations and contrast with experimental results for pressure quenches. Here, we are mostly concerned with how to extend our knowledge of nonequilibrium properties of field theories. We will thus be using Zurek and Laguna’s model as a testing ground for our methods, comparing some of our results to theirs. Even though their freeze-out time is clearly not the same as our peak time, since the former measures when the system first falls out of equilibrium while the latter measures when the system is maximally out of equilibrium, both time variables scale similarly with cooling time-scale.
5. Simulation results for time-dependent self-interacting potential at constant temperature In addition to testing our second-order staggered leapfrog Langevin code for stability with the same parameters as [10] (L = 2048, δx = 0.125, δt = 0.025 and T = 0.01), we tested the reliability of our code by reproducing their defect density results. For a more thorough description of the algorithm used, see [8]. The calculation of the structure function uses an FFT routine from Numerical Recipes [13], while the change in the structure function uses a simple finite difference method, with the ability to average over a few time steps to smooth the data. We can then extract the amplitude and time of the peak. Fig. 2 shows the variation in (t/τq ) with damping rate and cooling time-scale. t/τq is used on the time axis for easier comparison between different cooling time-scales. The amplitude of the peak is higher the lower the damping rate and the shorter the cooling time-scale (faster cooling). Also, (t/τq ) peaks earlier for larger cooling time-scales (slower cooling). The simulations start with the field φ thermalized at t/τq = −2, so that (t) is initially zero. As time progresses the potential changes and, if the thermalization rate of the field is not faster than the potential’s quenching rate, (t)
128
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
Fig. 2. (t) for different damping rates η and cooling time-scales τq with the potential in Eq. (1). The amplitude of the peak, Apeak , increases with decreasing time-scale and decreases with increasing damping rate. tpeak /τq increases with decreasing cooling time-scale and increasing damping rate.
begins to change, so that, even before the phase transition, it is nonzero and increasing. The most drastic changes in (t) happen as the potential goes from a single well to a double well, as we would expect: (t) peaks some time after the critical point. The peaks occurring during the cooling show a plateau beginning after the peak and lasting until the cooling ends at ' = 1. This can be seen in Fig. 2 for τq = 48, represented by the dashed line. For η = 1,
Fig. 3. Peak time vs. cooling time-scale τq for time-dependent potential for several damping rates. The solid straight lines above and below the data show the scaling predicted in Ref. [9] for overdamped and underdamped freeze-out times as a function of cooling time-scale. The simulation data interpolates between these two extremes appropriately.
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
129
Fig. 4. Amplitude of peak vs. cooling time-scale for time-dependent potential with varying damping rate. For all but the fastest quench times, −6/5 all damping rates studied satisfy the scaling Apeak ∝ τq . All three fit lines have a slope of −6/5.
we were able to identify a peak up to τq ∼ 128 at which point the data become too noisy to distinguish the peak from the plateau. For the fastest coolings (which, with η = 1, happen for τq 16), (t/τq ) actually peaks after the potential stops changing (' = t/τq = 1): the system remains out of equilibrium during the whole cooling process, which is thus equivalent to an instantaneous quench. For τq 256, (t) 0, and the system remains thermalized during the cooling; this is equivalent to an adiabatic cooling. It is between these two regimes, which, for η = 1, are given by 16 τq 256, that the cooling most dramatically affects the dynamics of the system. The location (tpeak ) and amplitude (Apeak ) of the peak in (t) are plotted in Figs. 3 and 4, respectively, as a function of τq for several damping rates. From Fig. 3 we note that the location of the peak scales with τq , in agreement with the scaling obtained in Ref. [10] for the freeze-out time as a function of τq . This is consistent with the fact that tpeak is the time after which the field is able to begin relaxing to its equilibrium state. −p We found that the amplitude of the peak, Apeak , scales as τq , where p = 1.2 ± 0.05 for all but the fastest coolings. We conclude that there are essentially three regimes: (i) the field remains thermalized all through the cooling process, thus approximating an “adiabatic” cooling (for τq 256, with η = 1); (ii) the field remains out of equilibrium throughout the duration of the cooling, thus approximating an instantaneous quench (for τq 16 and η = 1); (iii) the intermediate regime, where the field reaches thermalization during the cooling process, signaled by the appearance of a plateau in (t), which decreases in amplitude as τq is increased, until it cannot be differentiated from noise (for τq 256).
6. Summary We have presented a novel method to measure the departure from thermal equilibrium in nonlinear field theories, based on the computation of the momentum-integrated time-derivative of the structure function, which we called (t). In order to illustrate our method, we have discussed two particular cases: in the first one, we followed the evolution of a field in a quadratic potential prepared initially in a thermal state at T0 and quenched instantaneously
130
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
to T = 0. Due to the simplicity of the model, we were able to produce an analytical solution which was compared with the results from the numerical simulation. In the second case, we studied the Ginzburg–Landau potential with a linear pressure quench as suggested by Laguna and Zurek [10]. After reproducing their results, we showed how nonequilibrium is characterized by the presence of a peak in (t). The amplitude of the peak is a measure of the departure from equilibrium, while its location is related to the equilibration time-scale of the system as compared to the quenching time-scale. The general picture can be described qualitatively as follows: the modes are initially thermalized, and their rate of change increases as they try to “catch up” with the changing potential (or, when appropriate, with the changing heat bath). This may or may not happen during the cooling, determined by the value of τq . As the quench proceeds at the steady rate of τq−1 , (t) increases until the field begins to probe the bottom of the potential (the free-energy minima), when it has essentially “caught up” with the changing potential. If this happens while the potential is still changing, the field then changes at the same rate as the potential ((t) begins to plateau), until the potential stops changing and the field can fully thermalize ((t) → 0). For the fastest coolings (τq 16 for η = 1), the field is not able to “catch up” with the changing potential before the end of the cooling, and so (t) never plateaus; the amplitude of the peak does not obey the scaling in Fig. 4 because the field is still changing after the cooling has ended, which, as we remarked above, approximates an instantaneous quench. Finally, we note that our approach is applicable to any scalar field theory in the Ising universality class, with or without topological defects. It would be interesting to apply it to systems with metastable (stationary) states, as a tool to measure the lifetime and to quantify the approach to equilibrium in those systems. Also, it should be straightforward to obtain an expression of (t) for systems described by a complex scalar order parameter, such as 4 He. The extension of our approach to systems with local gauge invariance, applicable to the description of superconductors, remains open for investigation.
Acknowledgements CG was supported in part by a National Science Foundation grant no. PHY-9453431 and Dartmouth College. MG was supported in part by NSF grants PHY-0070554 and PHY-9453431, and by the Mr. Tomkins Fund for Cosmology and Field Theory at Dartmouth College. CG thanks the Theoretical Condensed Matter Group at Boston University, where parts of this work were developed, for their hospitality. We both thank Sidney Redner, Paul Krapivsky, Wonmuk Hwang, Prashant Sharma, Victor Spirin, Harvey Gould and Bill Klein for many helpful discussions. References [1] N.C. Wong, C.M. Knobler, Phys. Rev. Lett. 43 (1979) 1733; N.C. Wong, C.M. Knobler, Phys. Rev. Lett. 45 (1979) 498. [2] K. Binder, Phys. Rev. B 15 (1977) 4425. [3] N.C. Wong, C.M. Knobler, J. Chem. Phys. 69 (1978) 725. [4] R. Ruiz, Phys. Rev. A 26 (1982) 2227. [5] A. Onuki, Prog. Theoret. Phys. 66 (1981) 1230; A. Onuki, Prog. Theoret. Phys. 67 (1982) 768; A. Onuki, Prog. Theoret. Phys. 67 (1982) 787; A. Onuki, Phys. Rev. Lett. 48 (1982) 753. [6] J.D. Gunton, M. San Miguel, P.S. Sahni, in: C. Domb, J.L. Lebowitz (Eds.), Phase Transitions and Critical Phenomena, vol. 8, Academic Press, London, 1983. [7] J.S. Langer, in: C. Godreche (Ed.), Solids Far from Equilibrium, Cambridge University Press, Cambridge, 1992. [8] C.J. Gagne, M. Gleiser, Phys. Rev. E 61 (2000) 3483; M. Gleiser, H. Müller, Phys. Lett. B 422 (1998) 69.
C.J. Gagne, M. Gleiser / Physica D 181 (2003) 121–131
131
[9] A.J. Bray, Adv. Phys. 43 (1994) 357. [10] P. Laguna, W.H. Zurek, Phys. Rev. D 58 (1998) 085021; P. Laguna, W.H. Zurek, Phys. Rev. Lett. 78 (1997) 2519; W.H. Zurek, Phys. Rep. 276 (1996) 177; W.H. Zurek, Acta Pol. B 24 (1993) 1301; W.H. Zurek, Nature (London) 317 (1985) 505; Los Alamos Report No. LAUR 84-3818, Unpublished, 1984. [11] L.M.A. Bettencourt, N.D. Antunes, W.H. Zurek, Phys. Rev. D 62 (2000) 065005. [12] M. Bowick, A. Momen, Phys. Rev. D 58 (1998) 085014. [13] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed., Cambridge University Press, Cambridge, 1992.