# 1 Introduction

## 1  Introduction

The scaling equation of state near a critical point provides fundamental information on the critical behaviour of a thermodynamic system. For bulk Ising-like systems accurate predictions have been made recently  on the basis of the $`\phi ^4`$ field theory in three dimensions. Testing these predictions by Monte Carlo simulations  would be of considerable interest. Such simulations, however, are necessarily made only for finite systems and thus phenomenological concepts like finite-size scaling  are needed to perform extrapolations from mesoscopic lattices.

In order to test the theory in a more conclusive way it is desirable to go beyond the phenomenological finite-size scaling concept and to calculate explicitly the finite-size effects on the equation of state, i.e., on the magnetization as a function of the temperature $`T`$ and external field $`h`$ in a finite geometry. Although such field-theoretic calculations  are perturbative and not exact they have been found to be in good agreement with the MC simulations. So far the calculations of thermodynamic quantities were restricted to the case of zero external field $`h`$. At finite $`h`$, only the order-parameter distribution function was calculated . In the present paper we extend these calculations to the equation of state of the three-dimensional Ising model in a finite cubic geometry at finite $`h`$. For simplicity these calculations are performed at infinite cutoff and thus our results neglect lattice effects. The latter have been shown  to yield (exponentially small) non-universal non-scaling contributions. Here we focus our interest on the universal scaling part of the equation of state. Our theory predicts non-monotonic effects in the $`h`$ dependence of the equation of state which we then confirm with surprisingly good agreement by standard Monte Carlo simulations.

The identification of such non-monotonic effects is of great practical importance. If, for example, the critical temperature $`T_c(L)`$ of a lattice with $`L^3`$ sites in three dimensions varies asymptotically as $`T_c(L)T_c(\mathrm{})1/L^y`$ with some correlation length exponent $`y=1/\nu `$, then a plot of the numerical $`T_c(L)`$ versus $`1/L^y`$ gives the extrapolated $`T_c(\mathrm{})`$ as an intercept. If, however, higher order terms make the curve $`T_c(L)`$ non-monotonic, then such a plot for finite $`L`$, where the non-monotonicity is not yet visible, would give a wrong estimate for $`T_c(\mathrm{})`$. Similar effects would make estimates of other quantities unreliable, and there exist examples of such estimates in the literature.

## 2 Field-theoretic calculations

We consider the $`\phi ^4`$ model with the standard Landau-Ginzburg-Wilson Hamiltonian

$$H(h)=_V\left[\frac{1}{2}r_0\phi ^2+\frac{1}{2}(\phi )^2+u_0\phi ^4h\phi \right]$$
(1)   
where $`\phi (𝐱)`$ is a one-component field in a finite cube of volume $`V=L^d`$ and $`h`$ is a homogeneous external field. We assume periodic boundary conditions. Accordingly we have

$$\phi (𝐱)=L^d\underset{𝐤}{}\phi _𝐤e^{i𝐤𝐱}$$
(2)   
where the summation $`_𝐤`$ runs over discrete $`𝐤=\frac{2\pi }{L}𝐦`$ vectors with componets $`k_j=\frac{2\pi }{L}m_j,m_j=0,\pm 1,\pm 2,\mathrm{},j=1,2,\mathrm{},d`$ in the range $`\mathrm{\Lambda }k_j<\mathrm{\Lambda }`$, i.e., with a sharp cutoff $`\mathrm{\Lambda }`$. The temperature enters through $`r_0=r_{0c}+a_0t,t=(TT_c)/T_c`$.

As pointed out recently  the Hamiltonian (1) for periodic boundary conditions with a sharp cutoff $`\mathrm{\Lambda }`$ does not correctly describe the exponential size dependence of physical quantities of finite lattice models in the region $`\xi L`$. Instead of (1), a modified continuum Hamiltonian with a smooth cutoff  would be more appropriate. Even better would be to employ the lattice version of the $`\phi ^4`$ theory to describe the non-scaling lattice effects of finite Ising models in the region $`\xi L`$. In the present paper, however, we shall neglect such effects by taking the limit $`\mathrm{\Lambda }\mathrm{}`$ (see below).

The fluctuating homogeneous part of the order-parameter of the Hamiltonian (1) is $`\mathrm{\Phi }=V^1_Vd^dx\phi (𝐱)=L^d\phi _0`$. As previously  $`\phi `$ is decomposed as

$$\phi (𝐱)=\mathrm{\Phi }+\sigma (𝐱)$$
(3)   
where $`\sigma (𝐱)`$ includes all inhomogeneous modes

$$\sigma (𝐱)=L^d\underset{𝐤\mathrm{𝟎}}{}\phi _𝐤e^{i𝐤𝐱}.$$
(4)   
The order-parameter distribution function $`P(\mathrm{\Phi })P(\mathrm{\Phi },t,h,L)`$ is defined by functional integration over $`\sigma `$,

$$P(\mathrm{\Phi },t,h,L)=Z(h)^1D\sigma e^{H(h)},$$
(5)   
where

$$Z(h)=_{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Phi }D\sigma e^{H(h)}$$
(6)   
is the partition function of system. This distribution function depends also on the cutoff $`\mathrm{\Lambda }`$ which implies non-scaling finite-size effects . From the order-parameter distribution function $`P(\mathrm{\Phi })`$ we can calculate the magnetization

$$M=<|\mathrm{\Phi }|>=_{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Phi }|\mathrm{\Phi }|P(\mathrm{\Phi }).$$
(7)   
The functional integration over $`\sigma `$ in Eqs. (5) and (6) can only be done perturbatively . Recently a novel perturbation approach was presented . Using this approach the order-parameter distribution can be written in the form

$$P(\mathrm{\Phi })=e^{H^{eff}(\mathrm{\Phi })}/_{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Phi }e^{H^{eff}(\mathrm{\Phi })}$$
(8)   
where the (bare) effective Hamiltonian of the Ising-like system reads

$`H^{eff}(\mathrm{\Phi })`$ $`=`$ $`H_0(\mathrm{\Phi },h){\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐦\mathrm{𝟎}}{}}\mathrm{ln}\{{\displaystyle \frac{Z_1[y_{0𝐦}(r_{0L})]}{Z_1[y_{0𝐦}(r_0)]}}\},`$ (9)   
$`H_0(\mathrm{\Phi },h)`$ $`=`$ $`L^d({\displaystyle \frac{1}{2}}r_0\mathrm{\Phi }^2+u_0\mathrm{\Phi }^4h\mathrm{\Phi }),`$ (10)   
with $`r_{0L}=r_0+12u_0\mathrm{\Phi }^2`$ and $`y_{0𝐦}(r)=(2L^d/3u_0)^{1/2}(r+\frac{4\pi ^2}{L^2}𝐦^2)`$. The function $`Z_1[y]`$ in Eq.(9) is defined as

$$Z_1[y]=_0^{\mathrm{}}𝑑ss\mathrm{exp}(\frac{1}{2}ys^2s^4).$$
(11)   
After minimal renormalization at fixed dimension $`d<4`$  in the limit $`\mathrm{\Lambda }\mathrm{}`$, we obtain the finite-size scaling form of the order-parameter distribution function 

$$P(\mathrm{\Phi },t,h,L)=L^{\beta /\nu }p(hL^{\beta \delta /\nu },tL^{1/\nu },\mathrm{\Phi }L^{\beta /\nu }),$$
(12)   
where

$$p(x,q,z)=\frac{\mathrm{exp}[F(x,q,z)]}{_{\mathrm{}}^{\mathrm{}}𝑑z\mathrm{exp}[F(x,q,z)]},$$
(13)   
with $`x=hL^{\beta \delta /\nu }`$, $`q=tL^{1/\nu }`$, $`z=\mathrm{\Phi }L^{\beta /\nu }`$ and

$`F(x,q,z)`$ $`=`$ $`c_2(x,\widehat{q})\widehat{z}^2+c_4(x,\widehat{q})\widehat{z}^4xz`$ (14)   
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐦\mathrm{𝟎}}{}}\mathrm{ln}\{{\displaystyle \frac{Z_1[y_𝐦(\stackrel{~}{r}_L(x,\widehat{q},\widehat{z}))]}{Z_1[y_𝐦(\stackrel{~}{r}_L(x,\widehat{q},0))]}}\}.`$   
Here $`\widehat{q}=Q^{}t(L/\xi _0)^{1/\nu }`$ and $`\widehat{z}=(2Q^{})^\beta (\mathrm{\Phi }/A_M)(L/\xi _0)^{\beta /\nu }`$ are dimensionless variables normalized by the asymptotic amplitudes $`\xi _0`$ and $`A_M`$ of the bulk correlation length $`\xi =\xi _0t^\nu `$ at $`h=0`$ above $`T_c`$ and of the bulk order-parameter $`M_{bulk}=A_M|t|^\beta `$ at $`h=0`$ below $`T_c`$. The bulk parameter $`Q^{}`$ is known . The coefficients $`c_2(x,\widehat{q})`$ and $`c_4(x,\widehat{q})`$ read for $`d=3`$

$`c_2(x,\widehat{q})`$ $`=`$ $`(64\pi u^{})^1\widehat{q}\stackrel{~}{\mathrm{}}(x,\widehat{q})^{3(2\beta +1)/\nu }(1+12u^{}),`$ (15)   
$`c_4(x,\widehat{q})`$ $`=`$ $`(256\pi u^{})^1\stackrel{~}{\mathrm{}}(x,\widehat{q})^{34\beta /\nu }(1+36u^{}),`$ (16)   
where $`u^{}`$ is the fixed point value of the renormalized coupling . In three dimensions we have

$`y_𝐦(\stackrel{~}{r}_L(x,\widehat{q},\widehat{z}))`$ $`=`$ $`[6\pi u^{}\stackrel{~}{\mathrm{}}(x,\widehat{q})]^{1/2}[\stackrel{~}{r}_L(x,\widehat{q},\widehat{z})\stackrel{~}{\mathrm{}}(x,\widehat{q})^2+4\pi ^2𝐦^2],`$ (17)   
$`\stackrel{~}{r}_L(x,\widehat{q},\widehat{z})`$ $`=`$ $`\widehat{q}\stackrel{~}{\mathrm{}}(x,\widehat{q})^{1/\nu }+(3/2)\stackrel{~}{\mathrm{}}(x,\widehat{q})^{2\beta \nu }\widehat{z}^2.`$ (18)   
The auxiliary scaling function $`\stackrel{~}{\mathrm{}}(x,\widehat{q})`$ of the flow parameter is determined by

$`\stackrel{~}{\mathrm{}}(x,\widehat{q})^{3/2}`$ $`=`$ $`(4\pi u^{})^{1/2}[\stackrel{~}{y}(x,\widehat{q})+12\vartheta _2(\stackrel{~}{y}(x,\widehat{q}),\widehat{x})],`$ (19)   
$`\stackrel{~}{y}(x,\widehat{q})`$ $`=`$ $`(4\pi u^{})^{1/2}\stackrel{~}{\mathrm{}}(x,\widehat{q})^{3/21/\nu }\widehat{q},`$ (20)   
$`\widehat{x}`$ $`=`$ $`A_M(2Q^{})^\beta \xi _0^{\beta /\nu }(4\pi u^{})^{1/4}\sqrt{8}\stackrel{~}{\mathrm{}}(x,\widehat{q})^{\beta /\nu 3/4}x,`$ (21)   
where

$$\vartheta _2(\stackrel{~}{y},\widehat{x})=\frac{_0^{\mathrm{}}𝑑ss^2\mathrm{cosh}(\widehat{x}s)\mathrm{exp}(\frac{1}{2}\stackrel{~}{y}s^2s^4)}{_0^{\mathrm{}}𝑑s\mathrm{cosh}(\widehat{x}s)\mathrm{exp}(\frac{1}{2}\stackrel{~}{y}s^2s^4)}.$$
(22)   
From Eqs. (7) and (12) we obtain the scaling form

$`M(h,t,L)`$ $`=`$ $`L^{\beta /\nu }f_M(hL^{\beta \delta /\nu },tL^{1/\nu }),`$ (23)   
$`f_M(x,q)`$ $`=`$ $`{\displaystyle \frac{_{\mathrm{}}^{\mathrm{}}𝑑z|z|\mathrm{exp}[F(x,q,z)]}{_{\mathrm{}}^{\mathrm{}}𝑑z\mathrm{exp}[F(x,q,z)]}}.`$ (24)   
Correspondingly the asymptotic equation of state for a finite Ising-like system in the limit of zero lattice spacing  can be written as

$$h/M^\delta =f(hL^{\beta \delta /\nu },t/h^{1/\beta \delta })$$
(25)   
where

$$f(x,y)=x/f_M(x,yx^{1/\beta \delta })^\delta .$$
(26)   
In the comparison with the MC data of the simple-cubic (sc) Ising model in Section 4, the quantities $`h`$, $`M`$ and $`L`$ are used in a dimensionless form in units of the lattice constant, see Sect. 4 of Ref.. We have taken the bulk parameters $`u^{}=0.0412,Q^{}=0.945`$ from Ref.. From Ref. we have taken the bulk amplitudes $`\xi _0=0.495,A_M=1.71`$ in units of the lattice constant (of the sc Ising model) and the bulk critical exponents $`\beta =0.3305,\nu =0.6335`$. Thus our determination of the scaling function $`f(x,y)`$ does not require a new adjustment of nonuniversal parameters.

Taking the limit $`hL^{\beta \delta /\nu }\mathrm{}`$ at fixed $`t/h^{1/\beta \delta }`$, we obtain the scaling form of the bulk equation of state

$$h/M^\delta =f_b(t/h^{1/\beta \delta })=f(\mathrm{},t/h^{1/\beta \delta }).$$
(27)   
At $`T=T_c`$ we find from Eqs. (23)-(27)

$$h/M^\delta D_c=f(\mathrm{},0)=0.202$$
(28)   
in three dimensions (in units of the lattice constant).

## 3 Monte Carlo simulation

Standard heat bath techniques were used for the Glauber kinetic Ising model, with multi-spin coding (16 spins in each 64-bit computer word). Since the effect could be seen best in lattices of intermediate size $`L80`$ for $`L\times L\times L`$ spins, memory requirements were tiny and only trivial parallelization by replication, not by domain decomposition, was used. However, thousands of hours of processor time were needed since we see the non-monotonic effects clearly in our figure if $`h/M^\delta `$ is plotted. The exponent $`\delta `$ is nearly five and thus five percent accuracy in $`h/M^\delta `$ requires one percent accuracy in the directly simulated magnetization $`M`$.

Earlier simulations by the same author and algorithm  did not show the non-monotonic effects at $`T=T_c`$ because they were not searched for; at that time the field-theoretical predictions presented above were not yet known. But even if they had been known it is doubtful that with the Intel Paragon used in  instead of the Cray-T3E now the non-monotonic trends would have been seen in about the same computer time.

Because of the limited system size, errors in the magnetization of order 1 + const/$`L^{\beta /\nu }`$ 1 + const/$`\sqrt{L}`$ are expected; in light of these errors the agreement to be presented now is surprisingly good.

## 4 Results and discussion

Figure 1 compares our Monte Carlo results for $`32^3`$ and $`80^3`$ spins with our theoretical predictions, Eqs. (23)-(26), and shows good quantitative agreement; below $`T_c`$ a pronounced peak is seen in the scaling function (Fig.1a), at $`T_c`$ it is somewhat weaker (Fig.1b), above $`T_c`$ it is barely visible (Fig.1c), and far above $`T_c`$ (Fig.1d) it has vanished. The simulations, which at $`T=T_c`$ have an accuracy of the order of one and five percent for $`L=32`$ and 80, respectively, agree nicely with the theoretical predictions. (Far above $`T_c`$ therefore no simulations were made.)

Thus, if one tries to determine the bulk critical amplitude of $`h/M^\delta `$ at the critical isotherm, then with varying $`L`$ at a fixed field one first gets a too small value (left border of the figures), then a too high value (peak), and then a roughly correct value (plateau in the right part of the figure).

These non-monotonicities do $`not`$ vanish if we take the lattices large enough. They are part of the asymptotic scaling function and thus whatever lattice size $`L`$ we take there will be a field $`h1/L^{\beta \delta /\nu }`$ where $`h/M^\delta `$ has a maximum and near which it thus varies non-monotonically with $`L`$.

We also tested the prediction of  that for fixed $`T<T_c`$ the leading finite-size deviation from the bulk value of $`M`$ should vanish exponentially in $`L`$, and not with a power law $`1/L^d`$ (as predicted by perturbation theory based on the separation of the zero-mode \[4-7\]). Fig. 2a shows again non-monotonic behaviour, in both three and five dimensions. But only in three dimensions these data are accurate enough to distinguish between a tail varying exponentially and one $`1/L^d`$; Fig.2b clearly supports the theoretically predicted  exponential variation. (These three-dimensional data were taken with Ito’s fast algorithm .)

In order to make contact with bulk properties we also used simulations at the critical isotherm with $`1292^3`$ spins . From the simulations we obtain the bulk value of $`h/M^\delta `$ as $`D_c=0.21\pm 0.02`$. Our field-theoretic result in Eq. (28) is in very good agreement with this value.

It is interesting to compare our simulation result also with other bulk theories. From series expansion Zinn and Fisher  obtained $`C^c=0.299`$ for the amplitude of the bulk susceptibility at the critical isotherm $`\chi =C^c|h|^{\gamma /\beta \delta }`$ with $`\gamma =1.2395,\nu =0.6320`$. This leads to $`D_c=0.182`$ according to the relation $`D_c=(C^c\delta )^\delta `$.

$`D_c`$ is also contained in the universal combination of amplitudes 

$$R_\chi =\mathrm{\Gamma }D_cA_M^{\delta 1}$$
(29)   
where $`\mathrm{\Gamma }`$ is the amplitude of the bulk susceptibility $`\chi =\mathrm{\Gamma }t^\gamma `$ above $`T_c`$ at $`h=0`$ and $`A_M`$ is the amplitude of the spontaneous magnetization $`M_{bulk}=A_M|t|^\beta `$ below $`T_c`$. Using $`\phi ^4`$ field theory at $`d=3`$ dimensions Guida and Zinn-Justin  have obtained $`R_\chi =1.649`$. They also used the $`ϵ=4d`$ expansion and obtained $`R_\chi =1.674`$ at $`ϵ=1`$. Using these values for $`R_\chi `$ and the high-temperature series expansion results  $`\mathrm{\Gamma }=1.0928`$, $`A_M=1.71`$ and $`\delta =(d\nu +\gamma )/(d\nu \gamma )`$ with $`\gamma =1.2395,\nu =0.6335`$, we obtain from Eq. (29) $`D_c=0.205`$ ($`d=3`$ field theory) and $`D_c=0.202`$ ($`ϵ`$-expansion), respectively, in good agreement with our theoretical result, Eq. (28), and with our MC simulation.

In summary, our simulations confirmed a posteriori our theoretical predictions of Sect. 2 for the asymptotic finite-size effects in the three-dimensional Ising model. The agreement in Fig.1 is remarkable in view of the fact that the non-universal parameters of the theory were adjusted only to bulk parameters of the Ising model at $`h=0`$ and not to any finite-size MC data.

Acknowledgments

Support by Sonderforschungsbereich 341 der Deutschen Forschungsgemeinschaft and by NASA under contract numbers 960838, 1201186, and 100G7E094 is acknowledged. X.S.C. thanks the National Science Foundation of China for support under Grant No. 19704005, D.S. thanks the German Supercomputer Center in Jülich for time on their Cray-T3E.

Figure Captions 

Fig. 1.  Scaling plot (in units of the lattice constant, see Ref.) of $`h/M^\delta `$ versus $`x=hL^{\beta \delta /\nu }`$ below $`T_c`$ with $`t/h^{1/\beta \delta }=1.0`$ in part a, $`t/h^{1/\beta \delta }=0`$ (that means $`T=T_c`$) in part b, $`t/h^{1/\beta \delta }=1.0`$ in part c, and $`t/h^{1/\beta \delta }=1.6`$ in part d. Monte Carlo data for $`L=32`$ and $`80`$. Solid line is the theoretical prediction Eqs. (23)-(26); no Monte Carlo data are shown in part d, where $`h/M^\delta `$ is a monotonic function of external field $`x=hL^{\beta \delta /\nu }`$.

Fig.2.  a) Monte Carlo data for the spontaneous magnetization in units of the lattice constant in three (diamonds and plusses) and five (square) dimensions at $`T/T_c=0.99`$, versus linear lattice size $`L`$. The horizontal line $`M=0.3671`$ for three dimensions is determined from $`L=2496`$.

b) Selected three-dimensional data from part a) are shown as $`M(L=\mathrm{})M(L)`$ versus $`L`$ in a semilogarithmic plot. The straight line represents an exponential decay , the two curved lines are power law decays $`1/L^2`$ and $`1/L^3`$ which fail to fit the data. The five-dimensional data of part a) were not accurate enough to distinguish between an exponential and a $`1/L^5`$ decay.