# Probability distribution of the order parameter for the 3D Ising model universality class: a high precision Monte Carlo study

\[

## Abstract

We study the probability distribution $`P(M)`$ of the order parameter (average magnetization) $`M`$, for the finite-size systems at the critical point. The systems under consideration are the 3-dimensional Ising model on a simple cubic lattice, and its 3-state generalization known to have remarkably small corrections to scaling. Both models are studied in a cubic box with periodic boundary conditions. The model with reduced corrections to scaling makes it possible to determine $`P(M)`$ with unprecedented precision. We also obtain a simple, but remarkably accurate approximate formula describing the universal shape of $`P(M)`$.

\]

This work is devoted to the study of the following problem. Consider a finite system belonging to the universality class of the three-dimensional (3D) Ising model, exactly at its critical point. Let the system have a non-conserved order parameter, cubic symmetry and periodic boundary conditions. For such a finite-size system the order parameter $`M`$ (for the Ising model, the sum of all spins, divided by the total number of spins in the system) will be a fluctuating quantity, characterized by the probability distribution $`P(M)`$ . In the scaling limit (system size going to infinity) this function is universal (up to rescaling of $`M`$) and can be thus considered a very interesting and informative characteristic of the given universality class. (One should bear in mind that $`P(M)`$ depends on the geometry of the box, and on the boundary conditions; in this study we always consider a cubic box with periodic boundary conditions). For example, $`P(M)`$ contains the information about all momenta $`M^k`$ of $`M`$, including the universal ratios such as the Binder cumulant $`U=1(1/3)M^4/M^2^2`$, which has been a subject of many Monte Carlo studies . Moreover, the precise knowledge of $`P(M)`$ proved to be important for locating and characterizing the critical point in Monte Carlo studies of various systems, including the liquid-gas critical point , the critical point in the unified theory of weak and electromagnetic interactions  and in quantum chromodynamics .

The first Monte Carlo computation of $`P(M)`$ for the 3D Ising model in a cubic box with periodic boundary conditions has been performed in Ref. , where its double-peak shape was established. A more accurate determination of $`P(M)`$ has been done in Ref. , also by Monte Carlo. Results reported for the 3D case in Ref.  appear to be incorrect.

Despite considerable progress in computation of $`P(M)`$ by analytical methods , numerical simulation remains the main source of information about its properties.

Our aim was to compute $`P(M)`$ on a qualitatively new level of accuracy, in comparison to what has been done before , and to put the result into form convenient for further use. We would like to emphasize the following two features of our computation that made this possible.

1. The computation of Ref.  used the 3D Ising model on a simple cubic lattice of the size $`20^3`$ and $`30^3`$. As we will see, the shape of $`P(M)`$ obtained with these relatively small lattice sizes still differs noticeably from its scaling limit, due to nonnegligible corrections to scaling. To get over this difficulty, we employed, besides the 3D Ising model, the more sophisticated model in the same universality class, which was shown to have remarkably small corrections to scaling . This made it possible to determine the scaling limit of $`P(M)`$ with an accuracy far exceeding what would be achievable when one is restricted to the standard 3D Ising model.

2. The existing results for $`P(M)`$ were presented in the form of Monte Carlo-generated histograms . We present a simple 3-parameter formula which is suitable for quantitative applications. Its accuracy is about $`210^3`$ of the maximum value of $`P(M)`$.

We have performed Monte Carlo simulations of two models. The first one is the standard 3D Ising model on the simple cubic lattice, defined by the partition function

$$Z=\underset{\{s_i\}}{}\mathrm{exp}\left\{\beta \underset{ij}{}s_is_j\right\},s_i=\pm 1.$$
(1)   
Here $`ij`$ denotes the pairs of nearest neighbours, and the sum is over the $`2^N`$ possible configurations, where $`N`$ is the total number of spins. We simulate this model at the critical point, which we take to be at $`\beta _c=0.221654`$, using the Swendsen-Wang cluster algorithm , and lattice sizes ranging from $`12^3`$ to $`58^3`$ (with periodic boundary conditions).

The second model (with dramatically reduced corrections to scaling, as was shown in Ref. ) is the spin-1 Blume-Capel model . Here the spins can take 3 discrete values: $`1`$, 0, and $`+1`$. The model is defined by the partition function

$$Z=\underset{\{s_i\}}{}\mathrm{exp}\left\{\beta \underset{ij}{}s_is_jD\underset{m}{}s_m^2\right\},s_i=1,0,+1.$$
(2)   
The sum thus includes $`3^N`$ possible configurations. The parameter $`D`$ is fixed to the special value $`D=\mathrm{ln}2`$ (as explained in Ref. ), and we perform the simulations at the critical point, which is taken to be $`\beta _c=0.393422`$ , using lattice sizes from $`12^3`$ to $`58^3`$. The simulations used a hybrid method, which alternates one Metropolis sweep with 5 or 10 Wolff  steps, depending on the system size, as described in Ref. .

The probability distribution $`P(M)`$ is obtained as follows. For each configuration generated by the Monte Carlo algorithm, we determine the order parameter $`M=\frac{1}{N}_{i=1}^Ns_i`$, and increment the population of the corresponding bin of the histogram by 1.

We have found that the following ansatz gives a surprisingly good approximation to $`P(M)`$:

$$P(M)\mathrm{exp}\left\{\left(\frac{M^2}{M_0^2}1\right)^2\left(a\frac{M^2}{M_0^2}+c\right)\right\}.$$
(3)   
At the same time, the simpler ansatz

$$P(M)\mathrm{exp}\left\{c\left(\frac{M^2}{M_0^2}1\right)^2\right\}$$
(4)   
is clearly insufficient. This is illustrated in Fig. 1, using our highest-statistics data set for the $`20^3`$ lattice. One observes that the accuracy of approximation (3) is approximately 20 times higher than that of Eq. (4), and the residual discrepancy of Eq. (3) is comparable to the statistical noise, even with the high statistics used.

The ansatz (3) was motivated by the observation that $`M^6`$ plays an important role in the effective potential of the models in the 3D Ising universality class, while higher powers of the order parameter can usually be neglected . That is, the effective potential can in many cases well be approximated by a polynomial consisting of $`M^2`$, $`M^4`$ and $`M^6`$ terms. This is exactly what appears in the exponent in Eq. (3).

The approximate nature of the ansatz (3) manifests itself by its failure to correctly reproduce the large-$`M`$ behaviour of the tails of $`P(M)`$, which is governed by the critical index $`\delta `$,

$$P(M)M^{(\delta 1)/2}\mathrm{exp}\{constM^{\delta +1}\}$$
(5)   
(see Ref. ; for the discussion of the preexponential factor in a more general setting, see Ref. ). However, due to the fact that for the 3D Ising universality class the exponent $`\delta +15.8`$ is close to 6, this does not prevent the ansatz (3) from accurately describing the main part of $`P(M)`$ (excluding extremely-far-tail region).

The polynomial in the exponent of Eq. (3) has three parameters. Instead of simply parametrizing it by the coefficients in front of $`M^2`$, $`M^4`$ and $`M^6`$, we have chosen the parametrization so as to separate the scale-invariant parameters ($`a`$ and $`c`$) and the scale-dependent parameter $`M_0`$ (which parametrizes the position of the peak of the order parameter). The values of $`a`$ and $`c`$ in the scaling limit are universal and determine the “universal shape” of $`P(M)`$.

The results of our Monte Carlo simulations are collected in Tables I and II and shown in Figures 1, 2. For the spin-1 model, no deviations from scaling are observed on lattices $`16^3`$ and larger, while the simple cubic Ising model demonstrates pronounced corrections to scaling, which are, even on our largest lattices, much higher than both statistical errors of our spin-1 simulations and the accuracy of approximation (3). Corrections to scaling make it difficult to extract accurate scaling limit values of $`a`$ and $`c`$ from the simple cubic Ising model data, even if statistical errors are reduced by a higher-statistics simulation, due to necessity to extrapolate to $`L\mathrm{}`$.

There is no such problem with the spin-1 model, and we obtain the universal parameters of Eq. (3),

$$a=0.158(2),c=0.776(2).$$
(6)   
Here the errors take into account both statistical uncertainties and the systematic deviations inherent in the approximation (3). The latter are estimated from the lower right plot in Fig. 1.

From Eqs. (3) and (6) one can easily obtain any required property of $`P(M)`$. For example, one immediately learns that the ratio of the peak value of $`P(M)`$ to its value at $`M=0`$ is

$$e^c=2.173(4).$$
(7)   
Summarizing, we have computed, with a higher accuracy than previously available, the scaling limit form of the probability distribution $`P(M)`$ of the order parameter $`M`$ of systems with 3D Ising universality, in a cubic box with periodic boundary conditions. A convenient description of $`P(M)`$ is given by Eqs. (3) and (6), which deviates from the actual $`P(M)`$ by no more than $`210^3`$ times its maximum value (Fig. 1, right).

###### Acknowledgements.

We thank INTAS (grant CT93-0023) and DRSTP (Dutch Research School for Theoretical Physics) for enabling one of us (M.T.) to visit Delft University.