- Dynamic critical behavior of the worm algorithm for the Ising model
- Nonequilibrium critical dynamics of the ferromagnetic Ising model with Kawasaki dynamics
- Numerical Study on Aging Dynamics in the 3D Ising Spin-Glass Model. III. Cumulative Memory
- Exact n-spinon dynamic structure factor of the Heisenberg model
- Damage Spreading in a 2D Ising Model with Swendsen-Wang Dynamics
- Low-temperature nucleation in a kinetic Ising model with soft stochastic dynamics
- DROPLET DYNAMICS FOR ASYMMETRIC ISING MODEL
- High speed algorithms and simulations of the critical dynamics of the Ising Model
- The mixing time evolution of Glauber dynamics for the mean-field Ising model
- Dynamic critical exponents of the Ising model with multispin interactions

Dynamic structure factor of the Ising model with purely relaxational dynamics.

arXiv:cond-mat/0301588v1 [cond-mat.stat-mech] 30 Jan 2003

Pasquale Calabrese,a,? Victor Mart? ?n-Mayor,b, ? Andrea Pelissetto,c,? Ettore Vicarid, Scuola Normale Superiore and INFN, Sezione di Pisa, I-56100 Pisa, ITALIA b Departamento de F?sica Teorica I, Universidad Complutense de Madrid, ? ? E-28040 Madrid, ESPANA c Dipartimento di Fisica, Universit` di Roma La Sapienza a and INFN, Sezione di Roma, I-00185 Roma, ITALIA Dipartimento di Fisica, Universit` di Pisa and INFN, Sezione di Pisa, I-56100 Pisa, ITALIA a

a

d

Abstract

We compute the dynamic structure factor for the Ising model with a purely relaxational dynamics (model A). We perform a perturbative calculation in the ? expansion, at two loops in the high-temperature phase and at one loop in the temperature magnetic-?eld plane, and a Monte Carlo simulation in the high-temperature phase. We ?nd that the dynamic structure factor is very well approximated by its mean-?eld Gaussian form up to moderately large values of the frequency and momentum k. In the region we can investigate, k 5, 10, where is the correlation length and the zero-momentum autocorrelation time, deviations are at most of a few percent. PACS Numbers: 64.60.Ht, 05.70.Jk, 75.40.Gb, 75.10.Hk

Typeset using REVTEX 1

I. INTRODUCTION

The dynamic structure factor C(k, ) is a physically interesting quantity that can be directly measured in scattering experiments. Indeed, in neutron-scattering experiments and in Born approximation, C(k, ) is proportional to the cross section for inelastic scattering with momentum transfer k and energy transfer . At a continuous phase transition the structure factor shows a universal scaling behavior that depends on the dynamic universality class of the system. In this paper we consider the dynamic universality class of the threedimensional Ising model with purely relaxational dynamics without conservation laws, which is also known as model A, see Ref. [1]. This dynamic universality class should be appropriate to describe the dynamic critical properties of uniaxial magnetic systems in which the energy is not conserved due to the coupling to phonons and of alloys such as -brass at the orderdisorder transition; see, e.g., Ref. [1]. Note that this universality class does not describe the dynamic behavior of simple ?uids and mixtures at the liquid-vapor or mixing transitions because of additional conservation laws [1]. The model-A dynamics for the Ising universality class may also be relevant for the dynamics of quarks and gluons at ?nite temperature and ?nite baryon-number chemical potential ?. Indeed, using quantum chromodynamics, which is the current theory of strong interactions, one can argue that in the T ? ? plane there exists an Ising-like continuous transition at the endpoint of a ?rst-order phase transition line [2,3]. Model A (or model C if one takes into account the baryon-number conservation [4]) should describe the critical dynamics at this critical point (see also Ref. [5]). In this paper we compute the structure factor C(k, ) for the three-dimensional Ising universality class with purely relaxational dynamics (model A) in equilibrium. We consider the ?-expansion and compute C(k, ) to two loops in the high-temperature phase and to one loop in the whole temperature magnetic-?eld plane. In the high-temperature phase we also perform a Monte Carlo simulation, using the standard Ising model and the Metropolis dynamics. We ?nd that, for moderately large k and , C(k, ) is very well approximated by its mean-?eld (Gaussian) expression. In the high-temperature phase the ?eld-theoretical analysis and the simulation show that corrections to the mean-?eld behavior are less than 1% for k 5 and 10, where is the correlation length and is the zero-momentum autocorrelation time. In the low-temperature phase, on the basis of a one-loop ?eld-theoretical analysis, we expect slightly larger corrections, but still of the order of a few percent. Note that our study concerns the scaling behavior of C(k, ) in equilibrium, but it should be observed that similar conclusions have been obtained for the nonequilibrium dynamics in which one quenches a disordered system at Tc [6]. The paper is organized as follows. In Sec. II we de?ne the quantities that are computed in the following sections. We report a list of de?nitions together with some properties that are used in the calculation. In Sec. III we present our ?eld-theoretical results, obtained using the general formalism of Refs. [7C9]. Sec. IV is devoted to the presentation of the Monte Carlo results. In the Appendix we report some technical details.

2

II. DEFINITIONS AND BASIC OBSERVABLES

In this paper we consider the equilibrium dynamics for an Ising-like theory with scalar order parameter ?(r, t) at temperature T in the presence of a time- and space-independent external (magnetic) ?eld H. We consider the connected two-point correlation function of the order parameter G(r, t1 ? t2 ) ?(r, t1 )?(0, t2 )

conn ,

(1)

where we have assumed to be in equilibrium, so that the correlation function depends only on the di?erence t1 ? t2 . Then, we de?ne its Fourier transform G(k, t) with respect to r, G(k, t) = and the structure factor C(k, ) C(k, ) = dt eit G(k, t). (3) dd r eikr G(r, t), (2)

Here, we do not write explicitly the dependence on T and H that is always understood in the notation. Near the critical point correlations develop both in space and time. They can be characterized in terms of the second-moment correlation length and of the zero-momentum integrated autocorrelation time de?ned by 2 1 2d 1 2 dd r |r|2 G(r, 0) = ? dt G(0, t) = 1 ? G(k, 0) ?k 2 ,

k 2 =0

(4) (5)

1 C(0, 0), 2

where G(0, 0) is the static magnetic susceptibility. As is well known, for T Tc (Tc is the critical temperature) and H 0, and diverge. In the absence of a magnetic ?eld, ? |T ? Tc |? , ? |T ? Tc |?z ? z , (6)

where is the usual static exponent and z a dynamic exponent that depends on the considered dynamics. The static exponents for the three-dimensional Ising universality class are very well determined [10C15], see Ref. [16] for a review. Present-day lattice studies give estimates that can be summarized as follows [16]: = 1.2372(5), = 0.6301(4), = 0.0364(5), = 0.110(1). The exponent z depends on the dynamics. For model-A dynamics, estimates of z in three dimensions have been obtained by employing several methods. There exist ?eld-theoretical perturbative calculations in di?erent schemes [17C19] and Monte Carlo analyses that determine z by studying the equilibrium dynamics at Tc in ?nite volume [20,21], damage spreading [22,23], the critical relaxation from an ordered state [24,25], hysteresis scaling [26], and the short-time critical dynamics [27]. For experimental determinations see, e.g., Refs. [28,29]. The exponent z turns out to be slightly larger than two. For instance, 3

z 2.017 from the ?xed-dimension ?eld-theoretical expansion [19], z 2.02 from an analysis interpolating the 4 ? ? and 1 + ? expansions [17], and z 2.04 from Monte Carlo simulations. Near the critical point correlation functions show a scaling behavior. For the static structure factor, neglecting scaling corrections, we have [30,31] G(k, 0) gstat (y; x), (7)

where x a0 (T ? Tc )M ?1/ , y k 2 2 , M ? is the time-independent (we only consider the equilibrium dynamics) static magnetization, and a0 is a normalization factor that is ?xed by requiring that x = ?1 corresponds to the coexistence line. The magnetization M is related to T and H by the equation of state, which, in the critical limit, can be written in the scaling form H = b0 M f (x), (8)

where b0 is a nonuniversal constant which is ?xed by the condition f (0) = 1. The function gstat (y; x) has been extensively studied, both in the high-temperature [32,33] and in the low-temperature phase [34]; see Ref. [16] for an extensive review. In the high+ + temperature phase, the static function gstat (y) is known to O(?3 ) [35], and satis?es gstat (y) = 1 + y + O(?2 y 2 ). Its small-momentum expansion in three dimensions has been accurately determined using high-temperature expansion techniques, see, e.g., Refs. [15,16], ?nding

+ gstat (y) = 1 + y ? 0.000390(6) y 2 + 0.0000088(1) y 3 + O(y 4).

(9)

There are also precise estimates of the equation-of-state scaling function f (x) [36,11,15,37]. Eq. (7) can be extended to ?nite values of t. In the critical limit we can write G(k, t) g(y, s; x) (10)

with s t/ . We can also de?ne a scaling function for the structure factor: C(k, ) C(y, w; x), 2 where w and C(y, w; x) = 1 2 ds e?iws[g(y, s; x)]?1. (12) (11)

We also de?ne an integrated autocorrelation time at momentum k, (k) 1 2 dt G(k, t) G(k, 0) = C(k, 0) 2G(k, 0) , (13)

and an exponential autocorrelation time 4

exp (k) ? lim

|t|

ln G(k, t)

|t|

,

(14)

which controls the large-t behavior of G(k, t): G(k, t) ? exp[?|t|/exp (k)] for |t| . In the scaling limit, neglecting scaling corrections, (k) T (y; x) = C(y, 0; x)gstat(y), 1 exp (k) Texp (y; x) = , (k) |w0 (y; x)|C(y, 0; x)gstat(y) (15) (16)

where iw0 (y; x) are the zeros of [C(y, w; x)]?1 at ?xed y and x on the imaginary w-axis that are nearest to the origin w = 0. For a Gaussian free theory we have C(k, )|Gaussian = 2?m2 , ?2 (m2 + k 2 )2 + 2 (17)

where ? is an Onsager transport coe?cient and m 1/. It follows [C(y, w; x)]?1 = (1 + y)2 + w 2 , [T (y; x)]?1 = 1 + y, Texp (y; x) = 1.

(18)

For y 0 and w 0 the above-de?ned scaling functions have a regular behavior and one can write [C(y, w; x)]?1 = (1 + y)2 + w 2 +

m,n=0

cn,m (x)y n w 2m ,

[T (y; x)]?1 = 1 + y +

n=0

tn (x)y n , texp,n (x)y n , (19)

Texp (y; x) = 1 +

n=0

with c0,0 (x) = 0 because of the de?nition of . The expansion coe?cients cn,m (x), tn (x), and texp,n (x) parametrize the deviations from the Gaussian behavior (18) in the low-frequency and low-momentum regime. At the critical point, T = Tc , H = 0, the structure factor obeys the scaling law C(k, ) = 1 (2?+z)/z fC (k ?z ), (20)

with fC () ?nite, which implies that, for y , w keeping u wy ?z/2 = u0 k ?z ?xed, we have C(y, w; x) = f0 f (u/u0), (2?+z)/z C w 5 (21)

where f0 is a normalization constant. For large w at y and x ?xed we have C(y, w; x) c (y; x)w ?(2?+z)/z , (22)

where c (y; x) is ?nite and x-independent for y . The large-frequency behavior of the structure factor allows us to compute the nonanalytic small-s behavior of g(y, s; x) at y and x ?xed. We obtain for s 0+ [38]

(2?)/z [g(y, s; x)]?1 , nonanalytic = g0 (y; x)s

(23)

where 1 g0 (y; x) = ? c (y; x) sin(§¸/2)(1 ? ), (24)

with 1 + (2 ? )/z. Notice that, since (2 ? )/z 0.96, the nonanalytic small-t behavior of G(k, t) turns out to be practically indistinguishable from the analytic background.

III. FIELD-THEORETICAL RESULTS A. Field-theoretical approach

In order to determine the critical behavior of a purely relaxational dynamics without conservation laws, the so-called model-A dynamics, one may start from the the stochastic Langevin equation [1] H(?) ??(r, t) = ?? + (r, t), ?t ?(r, t) (25)

where ?(r, t) is the order parameter, H(?) is the Landau-Ginzburg-Wilson Hamiltonian H(?) = dd x 1 1 1 (??)2 + r?2 + u?4 ? H? , 2 2 4! (26)

? a transport coe?cient (cf. Eq. (17)), and (t) a Gaussian white noise with correlations (r, t) = 0, (r1 , t1 )(r2 , t2 ) = ?(r1 ? r2 )(t1 ? t2 ). (27)

The correlation functions generated by the Langevin equation (25) and averaged over the noise , can be obtained starting from the ?eld-theoretical action [7C9] S(?, ?) = ? dtdd x ? ? ?? H(?) + ?? ? ? ?? 2 ? lnJ(?) . ? ?t ? (28)

The last term in the action is an appropriate Jacobian term that compensates the contributions of self-loops of response propagators [8,9]. 6

In order to perform the ?eld-theoretical calculation it is useful to introduce the response function Y (r, t)it gives the linear response to an external magnetic ?eldde?ned by Y (r, t1 ? t2 ) = ?(r, t1 )?(0, t2 ) , ? (29)

(again we have assumed to be in equilibrium so that time-translation invariance holds), its Fourier transform Y (k, t) with respect to r, and its double Fourier transform R(k, ) with respect to r and t, de?ned as C(k, ) in Eq. (3). The response function and the two-point correlation function are strictly related. First, the zero-frequency response functions are related to the static correlation functions, G(k, 0) = ?R(k, 0). (30)

Moreover, because of the ?uctuation-dissipation theorem that holds for the equilibrium dynamics, we have C(k, ) = 2? Im R(k, ). Also the response function R(k, ) shows a scaling behavior and one can write r(y, w; x), ?R(k, ) neglecting scaling corrections. The function r(y, w; x) is such that r(y, 0; x) = 1 + y + O(y 2), r(0, w; x) = 1 ? iw + O(w 2), [r(y, ?w; x)]? = r(y, w; x). Then, it is easy to show by using Eqs. (30) and (31) that r(y, 0; x) = gstat (y; x), Im r(y, w; x) . C(y, w; x) = ? w|r(y, w; x)|2 For a Gaussian theory r(y, w; x) = 1 + y ? iw. (35) (32) (31)

(33)

(34)

The behavior of r(y, w; x) for small w and large w is similar to that of C(y, w; x). For small frequencies and momenta, the scaling function has a regular expansion in powers of w and y: r(y, w; x) = gstat (y; x) ? iw 1 + rn,m (x)(iw)m y n ,

n,m

(36)

where the coe?cients rn,m (x) are real and parametrize the w-dependent deviations from the Gaussian behavior (35). For w at ?xed y we have

+ r(y, w; x) r (y; x)(?iw)(2?)/z .

(37)

7

TABLE I. Numerical values of the coe?cients rn,m for 0 m, n 3. ?+ n\m 0 1 2 3 0 0 ?1.0487610?3 4.2337510?5 ?2.4853910?6 1 1.0312210?3 ?8.7216310?5 7.9321110?6 ?7.5141610?7 2 6.1941610?5 ?1.1584410?5 1.6810410?6 ?2.216810?7 3 6.5131610?6 ?1.9246610?6 3.8623610?7 ?6.5447110?8

B. Correlation functions in the disordered phase

In this Section we consider the equilibrium dynamics in the high-temperature phase H = 0, T > Tc . In order to determine the two-point correlation function, we have computed the scaling function r + (y, w) (here and in the following we will not indicate x and add instead a superscript + to remind the reader that we refer to the high-temperature phase) and we have then used Eq. (34). A two-loop calculation in the framework of the ? expansion gives

+ r + (y, w) = gstat (y) ? iw 1 + ?2 A(y, w) + O(?3 ) ,

(38)

where A(y, w) is reported in App. A. Note that A(0, 0) = 0 and A(y, ?w)? = A(y, w), as + expected from Eq. (33). The static function gstat (y) is known to O(?3 ) [35], and at order ?2 it reads

+ gstat (y) = 1 + y + 10?3 ?2 [?3.76012y 2 + 0.095966y 3 ? 0.00407101y 4 + O(y 5)] + O(?3 ) . (39) + Expanding A(y, w) in powers of y and w one obtains the coe?cients rn,m de?ned in Eq. (36). + We have rn,m = ?2 rn,m , where rn,m are reported in Table I for n, m 3. The coe?cients ?+ ?+ rn,m are rather small, the largest ones being of order 10?3 , and decrease quite rapidly. The ?+ analysis of the coe?cients of the expansion of A(k 2 , ) in powers of k 2 (at ?xed ) shows the presence of a singularity for w = ?3i. Therefore, we expect asymptotically

1 + rn,m rn,m?1 . ?+ ? 3

(40)

We have veri?ed numerically this relation, although quantitative agreement is observed only for quite large values of m: for n = 0, this relation is satis?ed at the 10% level only for m 41. Analogously, the coe?cients of the expansion of A(k 2 , ) in powers of become singular for k 2 = ?9, so that asymptotically 1 + ? rn,m ? rn?1,m . ?+ 9 (41)

The behaviors (40) and (41) can be interpreted in terms of the analytic structure of R+ (k, ). If one considers the structure factor, it is well known [39,35] that the nearest singularity [40] appearing in [G(k, 0)]?1 is the three-particle cut at k = 3imexp , where mexp is the mass 8

TABLE II. Numerical values of the coe?cients c+ for 0 n, m 3. ?n,m n\m 0 1 2 3 0 0 0.00104876 ?0.0170494 0.00106254 1 0.00212438 0.000951544 ?0.000075777 3.4320110?6 2 ?0.000075868 9.7305910?7 1.1333410?6 ?2.2254810 ?7 3 1.22037 10?6 ?1.8765710?7 ?1.1313110?8 1.0345710 ?8

gap of the theory. Since in the critical limit mexp 1 [16] with very small corrections (more precisely mexp ? 1 = ?2.00(3) 10?4 , see Ref. [15]), the nearest singularity to the origin + appearing in gstat (y) corresponds to y ?9. In view of relation (41), it is natural to conjecture that the same behavior holds for R+ (k, ), so that Eq. (41) should approximately hold + for the three-dimensional coe?cients rn,m and not only for their two-loop approximation. Relation (40) is consistent with the idea that the three-particle cut also controls the smallw behavior. In this case it is natural to conjecture that the coe?cients of the expansion of [R+ (k, )]?1 in powers of k 2 have a singularity for = ?3i/exp (0). Thus, turning to the scaling function r + (y; w), we expect a singularity at w = ?3i /exp (0) ?3i, since, as we shall see, in the critical limit /exp (0) 1. Therefore, we expect relation (40) to be a + general property of the three-dimensional coe?cients rn,m. This discussion indicates that C(y, 0) and its w-derivatives at w = 0 have a convergent expansion in y for |y| 9 and analogously that C(0, w) and its y-derivatives at y = 0 have a convergent expansion for |w| 3. Mathematically, this does not tell us much about the convergence of the double expansion which requires to know the singularity structure for both y, w = 0. At two loops, one can easily verify from the exact expression that C(y, w) has a convergent double expansion in the whole region |w| < 3, |y| < 9, and it is sensible to conjecture that the same is true for the exact expansion. From the results of Table I, one sees quite clearly that the response function R+ (k, ) is well described by the Gaussian approximation for |w| 3 and |y| 9. Deviations should be smaller than 1% in this region. This result is very similar to that obtained for the static structure factor: in that case hightemperature expansions and Monte Carlo simulations [33] show that the deviations from the Gaussian behavior are less that 0.3% for y 9. + We now consider the large-frequency behavior. At order ?2 the function r (y) de?ned in Eq. (37) turns out to be constant and given by

+ r (y) = 1 + 0.00538992?2 + O(?3 ).

(42)

Again the correction term is quite small. Using the ?uctuation-dissipation theorem (34), we obtain for the scaling function + C (y, w): [C + (y, w)]?1 = gstat (y)2 + w 2 + ?2 E(y, w) + O(?3 ), where 9 (43)

E(y, w) = 2w(1 + y)Im A(y, w) + w 2 ? (1 + y)2 Re A(y, w).

(44)

We can then obtain the small-w and small-y behavior. For the coe?cients c+ , see Eq. (19), n,m we obtain c+ = ?2 c+ , where the constants c+ are reported in Table II for n, m 3. ?n,m ?n,m n,m Again, we should note that the coe?cients c+ are very small and show the same pattern ?n,m observed for rn,m . We expect that C + (y, w) has singularities at y = ?9 and w = 3i, so ?+ that |cn,m /cn+1,m | |cn,m /cn,m+1 | 9. Thus, in complete analogy with what observed for the static structure factor and R+ (k, ), the dynamic C + (k, ) is essentially Gaussian in the region y < 9 and |w| < 3. ? ? + We also compute the large-frequency behavior. For the coe?cients c+ (y) and g0 (y), see Eqs. (22) and (23), we obtain c+ (y) = 1 ? 0.00538992?2 + O(?3 ), + 2g0 (y) = 1 + 0.00136716?2 + O(?3 ), (45) (46)

where, since the corrections are very small, one may simply set ? = 1 to obtain a threedimensional numerical estimate. Therefore, for large w we predict C + (y, w) 0.995/w 1.95, + which is not very di?erent from the purely Gaussian behavior CGauss (y, w) 1/w 2 . Thus, the Gaussian approximation should be a reasonably good approximation even outside the small-w region, w < 3, discussed above. Trusting the above estimate of c+ (y) we ?nd that ? + C + (y, w)/CGauss (y, w) = 1.12, 1.25, 1.41 respectively for w = 10, 100, 1000. Thus, quite large values of w are needed in order to observe a signi?cant di?erence. + Finally we compute the scaling function Texp (y) de?ned in Eq. (16). For this purpose we need to compute exp (k) and therefore the large-t behavior of G(k, t). Because of the ?uctuation-dissipation theorem, it is equivalent to consider Y (k, t). For y < 3 we obtain Y (k, t) e?s(1+y) {1 ? ?2 sA[y, ?i(1 + y)]}, where s t/ , while for y > 3 we have Y (k, t) e?s(1+y) + ?2 27 e?s(y/3+3) [1 + O(s?1/2 )] + O(?3 ) 4 (y ? 3)2 (y + 9)3 8 s (48) (47)

For y < 3 the correction term exponentiates as expected, and as a consequence we obtain

+ + C + (y, 0)[gstat(y)]2 Texp (y) = 1 + ?2 A[y, ?i(1 + y)] + O(?3 ).

(49)

On the other hand, for y > 3 the correction term decreases with a di?erent exponential factor which dominates for large values of t, suggesting that, at leading order in ?, exp (k)/ = (y + 9)/3. In other words, the interaction turns on a new singularity (a three-particle cut) that becomes the leading one for y large enough. However, this is not the end of the story. Indeed, by considering graphs in which one recursively replaces each line with a two-loop watermelon graph one obtains contributions to Y (k, t) decreasing as exp[?s(3?n y + 3n )] (3n -particle cut) which would be more important for y large enough. These singularities will not probably be the only ones, since we also expect a 5-particle cut, a 7-particle cut, etc. + On the basis of these results, we expect Texp (y) to have several singularities on the positive 10

?exp,n for n 2 and m 3. TABLE III. Numerical values of the coe?cients rn,m , c? , t? , and t? ?? ?n,m ?n m r0,m ?? r1,m ?? r2,m ?? c? ?0,m c? ?1,m c? ?2,m ?m t? ? ? texp,m 0 0

5 ? 192 11 1920

1

1 48 3 ? 320 29 8960 3 64 7 1920 743 ? 107520 5 192 3 ? 65 + 2 ln 2 64

2

1 192 7 ? 1920 37 21504 17 ? 1920 221 71680 1 ? 4032 23 1920 1551 7 640 ? 2 ln 2

3

1 640 1 ? 672 115 129024 69 71680 251 ? 322560 949 2838528 697 ? 215040 ? 422211 + 17 71680 2

0

5 192 19 640

0

3 8

?

1 2

ln 2

ln 2

real y axis and to become eventually in?nite as y . This is not unexpected since, for y , R+ (k, ) behaves as ?(2?)/z and therefore has a branch cut starting at = 0. For y < 3, we can use Eq. (49) to compute the coe?cients t+ exp,n de?ned in Eq. (19). We + + 2 2 2 obtain, at order ? , texp,0 = 0.00110075? , texp,1 = 0.00337789? , t+ = 0.000217173?2, etc. exp,2 The coe?cients decrease as t+ /t+ exp,n exp,n+1 3, which re?ects the presence of a singularity at y = 3. Again, for y < 3 the deviations from a purely Gaussian behavior are very small.

C. Correlation function in the (t, H) plane

In the presence of an external magnetic ?eld H, a one-loop calculation gives r(y, w; x) = gstat (y; x) ? iw 1 + ? where B(y, w) is de?ned in App. A and 4+y+ y y 2? 4+y ?1 ? + ln gstat (y; x) = 1 + y + + O(?2 ). 3+x 12 2 y 4+y? y (51) 2 B(y, w) + O(?2 ) , 3+x (50)

Note that the O(?) correction vanishes for x in agreement with the results of the previous Section. Moreover, the x-dependence is very simple and in Eqs. (50) and (51) is always given by the prefactor 2/(3 + x) that becomes 1 on the coexistence curve x = ?1. As a consequence, such a prefactor will always appear in this Section, multiplying the lowtemperature results that will be speci?ed by adding a superscript ? to all de?nitions. Of course, such a simple x-dependence does not hold at higher loops, as it can be seen for instance from the two-loop results of Ref. [34] for the static structure factor. One can easily derive the small-momentum and small-frequency behavior by expanding the function B(y, w). The coe?cients rn,m (x), see Eq. (36), are given by rn,m ?? + O(?2 ), rn,m (x) = 2? 3+x 11 (52)

where rn,m are given in Table III for m 3 and n 2. ?? Again, we note that the corrections to the Gaussian behavior are small, although a factor-of-ten larger than the corresponding high-temperature ones. For instance, r0,1 0.02 ?? to be compared with r0,1 0.002. Moreover, the coe?cients decrease slower with n and m. ?+ This fact can be understood in terms of the singularities of the function B(y, w). A simple analysis shows the presence of singularities for w = ?2i and y = ?4, so that asymptotically 1 ? rn,m rn,m?1 , ?? ? 2 1 ? rn,m ? rn?1,m . ?? ? 4 (53)

This behavior can be understood on general grounds. Considering the static structure factor, it is known that the nearest singularity in the low-temperature phase is the two-particle cut ? k = 2imexp , so that gstat (y) has a singularity for y = ?4(mexp )2 ?4, where we have used the fact that in the critical limit mexp 1 (more precisely, mexp 0.96(1) [41,11]). As we did for the high-temperature phase we can thus conjecture that also the singularities of the dynamic functions are controlled by the two-particle cut. Therefore, we expect singularities for y = ?4(mexp )2 ?4 and w = ?2i /exp (0) ?2i, where we have used the fact that /exp (0) 1, with corrections of order a few percent as discussed below, in the critical limit. Therefore, Eq. (53) should also approximately apply to the three-dimensional coe?cients ? rn,m . The above-reported discussion shows that in the region |y| 4, |w| 2 the response function can be reasonably approximated by a Gaussian form. Note however that, while in the high-temperature phase corrections are expected to be less than 1%, here deviations should be larger. We have also studied the large-frequency behavior. The coe?cient r (y; x) turns out to be y-independent at one loop: r (y; x) = 1 ? ? + O(?2 ). 3+x (54)

Note that the correction is quite large and thus signi?cant deviations for the Gaussian behavior should be observed as soon as w is large. Using the ?uctuation-dissipation theorem, we can compute at one loop the scaling function C(y, w; x). For the small-w, small-y coe?cients, we obtain cn,m (x) = 2? ? c ? + O(?2 ), 3 + x m,n 2? ?? tn (x) = t + O(?2 ). 3+x n

(55)

?n The coe?cients c? and t? are reported in Table III for n 2 and m 3. ?m,n We have also investigated the large-frequency behavior. It is very simple to show, using the above-reported formulas, that at this order c (y; x) = 1/r (y; x) and g0 (y; x) = c (y; x)/2. Finally, we consider Texp (y; x). For this purpose we need to compute the large-t behavior of Y (k, t). We observe a behavior analogous to that observed in the high-temperature phase. For y < 2, 12

TABLE IV. Results of the Monte Carlo simulations. L Nit (a) 64 0.215 30106 4.4598(9) 19.38(11) (b) 64 0.219 8106 8.081(6) 64.9(9) (c) 128 0.2204 9106 13.050(7) 176(4) (d) 128 0.2210 4106 19.739(14) 420(23)

C(y, 0; x)[gstat (y; x)]2Texp (y; x) = 1 +

2? B[y, ?i(1 + y)] + O(?2 ), 3+x

(56)

while for y > 2 the two-particle cut contribution dominates so that exp (k)/ = 2 + y/2. The discussion reported in Sec. III B can be repeated also here. One can easily identify diagrams that decrease as exp[?s(2?n y + 2n )], indicating that Texp (y; x) has an in?nite number of singularities on the y axis and that it diverges for y . For small y, we can use Eq. (56) to compute the small-y expansion coe?cients texp,n (x). We have texp,n (x) = 2? ?? t + O(?2 ). 3 + x exp,n (57)

Numerical values are reported in Table III. Note that exp (0)/ = 1 + t? 1 + 0.0284?, exp,0 and thus we expect this ratio to be 1 with corrections of order of a few percent.

IV. MONTE CARLO RESULTS

We determine the dynamic structure factor C(k, ) and the scaling function G(k, t) in the high-temperature phase H = 0, T > Tc for small values of kas we shall see, we are able to reach k 10/¦±by means of a large-scale Monte Carlo simulation. We consider the Ising model on a cubic lattice, i.e. the Hamiltonian H = ? i j ,

ij

(58)

where 1/T , i = 1, and the summation is over all nearest-neighbor pairs ij . We measure the correlation function G(k, t) = 1 (eiqx + eiqy + eiqz ) 0,0,0 (t = 0)x,y,z (t) , 3 x,y,z (59)

for four di?erent values of L and : (a) L = 64, = 0.215; (b) L = 64, = 0.219; (c) L = 128, = 0.2204; (d) L = 128, = 0.221. Of course, in Eq. (59) q = 2n/L where n is an integer. For each and L we ?rst reached equilibrium by running 20000 Swendsen-Wang iterations, then we collected Nit iterations using the Metropolis algorithm [42]. The results of the simulations are reported in Table IV. There we report the number of iterations Nit , the second-moment correlation length (for the L = 128 lattices we report 13

the more precise results of Ref. [33]) and the autocorrelation time . Note that all lattices have L/ 6, a condition that usually ensures that ?nite-size e?ects are reasonably small (for static quantities corrections are less than 1%). The correlation length has been determined by using a discretized form of Eq. (4): 2 = /F ? 1 , 4 sin2 (/L) (60)

where F = G(k, 0) with k = (2/L, 0, 0). The integrated autocorrelation time , and also the autocorrelation times (k) considered below, have been determined using the self-consistent method of Ref. [43]: 1 (k) = + 2

M (k)

G(k, t) G(k, 0)

,

(61)

t=1

where t is the Monte Carlo time in sweeps and the cuto? M(k) is chosen self-consistently so that 6 (k) < M(k) 6 (k) + 1. Since G(k, t) decays exponentially, this choice makes the systematic error due to the truncation small, keeping the statistical variance small at the same time; see Ref. [43] for a discussion. First, we check that z |T ?Tc |?z . Using the precise estimate c = 0.22165459(10) of Ref. [12], we obtain from a least-square ?t z = 2.10(2) including all data and z = 2.11(5) discarding the estimate of for lattice (a). This result is in reasonable agreement with the estimates reported in Sec. II, if we take into account that we quote here only the statistical error. The systematic error due to corrections to scaling and to neglected ?nite-size e?ects is probably larger. Then, we determine the correlation function G(k, t). In Fig. 1 we report the function f + (y, s) G(k, t) g + (y, 0) , + (y, s) g G(k, 0) (62)

for three di?erent values of y k 2 2 , y = 0, 4, 16, as computed from lattices (a), (b), and (c). We have not included the results for lattice (d), because they have much larger errors. In order to obtain G(k, t) for a given k = 2n/L, we have performed a linear interpolation, using two nearby values of k. First, we observe reasonable scaling: corrections due to the ?nite values of and L are under control, although they increase as y increases. For y = 0 the results for the three di?erent lattices agree within a few percent, while for larger values of y we observe larger discrepancies. In particular, for y = 4 and y = 16, the estimates of f + (y, s) obtained from lattice (b) are always larger than those obtained from (a) and (c), the discrepancy being of order 20% when f + (y, s) 10?1 and 80% when f + (y, s) 10?2 . These di?erences are probably ?nite-size e?ects, since (a) and (c) have L/ 10, while L/ 8 for (b). It is also remarkable that the plot of ln f + (y, s) is a straight line, indicating that f + (y, s) is quite precisely a pure exponential. No deviations can be observed in Fig. 1. Therefore, G(k, t) G(k, 0) exp [?t/exp (k)] , 14 (63)

0

log10f

y=0, =0.2204, L=128 y=0, =0.2190, L=64 y=0, =0.2150, L=64 y=4, =0.2204, L=128 y=4, =0.2190, L=64 y=4, =0.2150, L=64 y=16, =0.2204, L=128 y=16, =0.2190, L=64 y=16, =0.2150, L=64

+

-1

-2

0

1

2

3

4

FIG. 1. The scaling function f + (y, s). We report results for lattices (a), (b), and (c) and for three di?erent values of y.

t/

15

0

-1

=0.2150, L=64 =0.2190, L=64 =0.2204, L=128 =0.2210, L=128 Gaussian

log10T

-2 -3 0

+

50

100

150

200

y

FIG. 2. Scaling plot of T (y) vs y k2 2 .

within the precision of our results. This behavior appears to be well satis?ed in the region that we can safely investigate, i.e. 1/10 t/ (k) 4 and k 5. Therefore, the dynamic structure factor is well approximated by a Lorentzian in the region of not too large frequencies, i.e. for (k) < 10. ? Then, we consider the scaling function T (y) that encodes the k-dependence of (k). In Fig. 2 we report our numerical results. Again, we observe good scaling up to quite large values of y. In the ?gure, we also report the Gaussian prediction T (y) = 1/(1 + y). It can be seen that the Gaussian approximation describes very well the numerical data. This result should have been expected on the basis of the results of Sec. III where we showed that the deviations from a Gaussian behavior are very small in the small-y regime y < 9, and should ? remain small even for larger y. For instance, using the data with largest y reported in Fig. 2, we estimate T (y) = 0.0053(3) for y = 181, to be compared with the Gaussian prediction 0.0055. Thus, in the range y 200 the discrepancy should be at most 4-10%. + Finally, we consider the function Texp (y). In order to compute exp (k) we de?ne an e?ective quantity e? (t; k) ? ln G(k, t + 1) G(k, t)

?1

.

(64)

The exponential autocorrelation time exp (k) is obtained from e? (t; k) by letting t go to in?nity. In practice, we can only compute e? (t; k) up to t of the order of (1-2) (k) since 16

1.08

=0.2204, L=128

1.06

eff((k);k)/(k)

1.04

1.02 1.00 0.98 0.96

0

10

20

30

40

50

y

FIG. 3. Ratio e? (t; k)/ (k) vs y = 0.2204. k2 2 , for t = (k). Results for lattice (c), L = 128,

17

errors increase rapidly. In Fig. 3 we report the ratio e? (t; k)/ (k) for t = (k) for lattice (c) which is the only one that allows us to reach large values of y. We observe e? (t; k) (k) within the precision of our results. It is tempting to conclude that exp (k) (k) for y < 50, but this is in contrast with the theoretical results of Sec. III B. Indeed, we showed there that exp (k) (k) with very small corrections for y < 3, but we noticed that this relation breaks down for larger values of y. For instance, for y > 3, our two-loop calculation gives exp (k)/ (k) = (3 + 3y)/(9 + y) which is signi?cantly larger than 1 for y > 3. As we already discussed this prediction should not be taken seriously, unless y is close to 3, since other singularities should be present, and indeed we expect exp (k)/ (k) to diverge as k . Therefore, our numerical data show that the asymptotic large-t behavior sets in only for large values of t, i.e. for t ? (k), where the correlation function G(k, t) is very small. Therefore, even if Eq. (63) breaks down for y 3 and t large, it still represents a very good approximation (even for y 50) for the values of t for which G(k, t) is sizeable.

ACKNOWLEDGEMENTS

V.M.-M. is a Ramon y Cajal research fellow and is partly supported by MCyT (Spain), project Nos. FPA2001-1813 and FPA2000-0956.

APPENDIX A: INTEGRALS ENTERING THE FIELD-THEORETICAL CALCULATIONS

In this appendix we report some integrals that enter the perturbative ?eld-theoretical calculations. In the two-loop calculation of the response scaling function in the high-temperature phase, cf. Eq. (38), one needs to compute the function A(k 2 , ) = 2 ?2 N [I(k, ) ? I(0, 0)] , 27 d

3

(A1)

with (dimensional regularization near four dimensions is understood)

I(k, ) =

0

dte

it

d d p1 d d p2 (2)d (2)d

p2 i=1 i

1 e?t +1

2 i (pi +1)

,

(A2)

where p3 = k ? p1 ? p2 , and Nd = 2/[(4)d/2 (d/2)]. The integral I(k, ) can be written in the form

?2 Nd I(k, )

1 1 = 2 3 (4)d Nd

1

t dt

0 0

2

dse s

?s 3?d is 1?t 3

1 0

e

e?s ? k ududv d/2 , ?

Q 2

(A3)

where 1 ? = t2 u[1 ? u + uv(1 ? v)] + (1 ? t2 ), 3 (1 ? t)2 (1 ? t)3 1?t 2 t u[1 ? u + uv(1 ? v)] + t+ . Q = t3 u2 (1 ? u)v(1 ? v) + 3 9 27 18 (A4) (A5)

We will also need the singularity structure of the A(k 2 , ). For this purpose, we will de? termine the large-t behavior of A(k 2 , t) which is the Fourier transform with respect to of 2 A(k , ). This behavior can easily be derived from Eq. (A2). Setting p1 = k/3 + q1 / t, p2 = k/3 + q2 / t, p3 = k/3 + q3 / t, we obtain that for t , 2 ?2 dd q1 dd q2 ?(q1 +q2 +q3 ) 1 2 2 2 2 ? A(k 2 , t) Nd t?d e?k t/3?3t 2 e [1 + O(t?1/2 )] 3 d (2)d 27 (k /9 + 1) (2) 1 (1/3)d/2 2 2 (d/2)t?d e?k t/3?3t [1 + O(t?1/2 )] = 2 3 (k /9 + 1) 54 1 3e?k t/3?3t 2 [1 + O(t?1/2 )] + O(?) (k + 9)3 2t4

2

(A6)

This result implies the presence of a branching cut in A(k 2 , ) starting at = ?i(3 + k 2 /3). The one-loop expression of the response function in the ordered phase, cf. Eq. (50), is written in terms of the function

?1 B(k 2 , ) = Nd [J(k, ) ? J(0, 0)] ,

(A7)

where

J(k, ) =

0

dte

it

d d p1 (2)d

2

p2 i=1 i

1 e?t +1

2 2 i (pi +1)

,

(A8)

with p2 = k ? p1 . The function B(k 2 , ) can be written in the form 1 B(k , ) = 2

2 1 0

i2(1 ? t) ? [1 ? t2 + 4t2 u(1 ? u)]k 2 t dt du . 4 + [1 ? t2 + 4t2 u(1 ? u)]k 2 ? 2i(1 ? t)

(A9)

Such an integral can be computed exactly obtaining 1 1 4 + 2i + k 2 2i 4 + k 2 4 + k2 ? k2 B(k , ) = ? + 2 ln + ln 2 k 4 ? 2i + k 2 k2 4 + k2 + k2 4 2 + (4 + k 2 )2 2 F ? F + + ik 2 1 + F ln + ln ? 2 ln k 16 k 2 F + F ? ? ik 2

2

(A10)

with F 2 + 2ik 2 ? k 2 (4 + k 2 ). It is easy to see using this exact expression or repeating the argument presented for A(k 2 , ) that B(k 2 , ) is singular for = ?i(2 + k 2 /2).

19

REFERENCES

? ? ?

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

Email address: Pasquale.Calabrese@df.unipi.it Email address: Victor@lattice.?s.ucm.es Email address: Andrea.Pelissetto@roma1.infn.it Email address: Ettore.Vicari@df.unipi.it P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977). M. A. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. Verbaarschot, Phys. Rev. D 58, 096007 (1998). J. Berges and K. Rajagopal, Nucl. Phys. B 538, 215 (1999). B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000). Sz. Bors?nyi, A. Patk?s, D. Sexty, and Zs. Sz?p, Phys. Rev. D 64, 125011 (2001). a o e P. Calabrese and A. Gambassi, Phys. Rev. E 65, 066120 (2002); E 66, 066101 (2002). P. C. Martin, E. D. Siggia, and H. A. Rose, Phys. Rev. A 8, 423 (1973). R. Bausch, H. K. Janssen and H. Wagner, Z. Physik B 24, 113 (1976). C. De Dominicis and L. Peliti, Phys. Rev. B 18, 353 (1978). R. Guida and J. Zinn-Justin, J. Phys. A 31 (1998) 8103. M. Campostrini, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. E 60, 3526 (1999). H. W. J. Blte, L. N. Shchur, and A. L. Talapov, Int. J. Mod. Phys. C 10, 1137 (1999). o M. Hasenbusch, Int. J. Mod. Phys. C 12, 911 (2001). P. Butera and M. Comi, Phys. Rev. B 65, 144431 (2002). M. Campostrini, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. E 65, 066127 (2002). A. Pelissetto and E. Vicari, Phys. Rep. 368, 549 (2002). R. Bausch, V. Dohm, H. K. Janssen, and R. P. K. Zia, Phys. Rev. Lett. 47, 1837 (1981). N. V. Antonov and A. N. Vasilev, Teor. Mat. Fiz. 60, 59 (1984) [Theor. Math. Phys. 60, 671 (1984)]. V. V. Prudnikov, A. V. Ivanov, and A. A. Fedorenko, Pisma Zh. Eksp. Teor. Fiz. 66, 793 (1997) [JETP Lett. 66, 835 (1997)]. S. Wansleben and D. P. Landau, Phys. Rev. B 43, 6006 (1991). H. Heuer, J. Phys. A 25, L567 (1992). U. Gropengiesser, Physica A 213, 308 (1995). P. Grassberger, Physica A 214, 547 (1995); (E) A 217, 227 (1995). D. Stau?er and R. Knecht, Int. J. Mod. Phys. C 7, 893 (1996). D. Stau?er, Physica A 244, 344 (1997). G. P. Zheng and J. X. Zhang, Phys. Rev. E 58, R1187 (1998). A. Jaster, J. Mainville, L. Sch lke, and B. Zheng, J. Phys. A 32, 1395 (1999). u D. P. Belanger, B. Farago, V. Jaccarino, A. R. King, C. Lartigue, F. Mezei, J. Phys. (France) Colloque C8, 49, 1229 (1988). M. T. Hutchings, M. P. Schulhof, and H. J. Guggenheim, Phys. Rev. B 5, 154 (1972). M. E. Fisher and R. J. Burford, Phys. Rev. 156, 583 (1967). H. B. Tarko and M. E. Fisher, Phys. Rev. Lett. 31, 926 (1973); Phys. Rev. B 11, 1217 (1975). M. Campostrini, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. E 57, 184 (1998). V. Mart? ?n-Mayor, A. Pelissetto, and E. Vicari, Phys. Rev. E 66, 026112 (2002). M. Combescot, M. Droz, and J. M. Kosterlitz, Phys. Rev. B 11, 4661 (1975). 20

[35] A. J. Bray, Phys. Rev. B 14, 1248 (1976). [36] R. Guida and J. Zinn-Justin, Nucl. Phys. B 489, 626 (1997). [37] J. Engels, L. Fromme, and M. Seniuch, Numerical equation of state from an improved three-dimensional Ising model, e-print cond-mat/0209492. [38] Here we use the fact that the integral d ||xexp(it) is given by ? ?1 sin(x/2)(1 + 2 x)|t|?1?x for x > ?1. [39] R. A. Ferrell and D. J. Scalapino, Phys. Rev. Lett. 34, 200 (1975). [40] Note that, beside the three-particle cut, G(k, 0) shows an additional singularity, the one-particle pole at k = imexp . However, such a point is a simple zerotherefore not a singularityof the inverse [G(k, 0)]?1 . [41] M. E. Fisher and S.-Y. Zinn, J. Phys. A 31, L629 (1998). [42] We used the Metropolis algorithm with sequential updating of the sites. At variance with the Metropolis algorithm with random updating, this dynamics is not reversible (it does not satisfy detailed balance). Nonetheless, it is commonly accepted that in the critical limit these two dynamics should both belong to the dynamic universality class of model A that describes a reversible di?usive dynamics without conservation laws. [43] N. Madras and A. D. Sokal, J. Stat. Phys. 50, 109 (1988).

21