# Three-dimensional Ising model in the fixed-magnetization ensemble: a Monte Carlo study

## I Introduction

Second order phase transitions display many interesting and subtle properties associated with scale invariance and universality at critical points. Some of these, such as power-law singularities of the free energy and other quantities at the critical point, and their critical exponents and amplitudes, have been studied rather thoroughly (see, e. g., ). Among less investigated items are the universal characteristics of finite-size effects. These are important for analysis of experiments with finite samples, as well as for computer simulations, which necessarily have to deal with finite-size systems.

In clear distinction from, for example, critical indices, the finite size effects depend crucially on the nature of the statistical ensemble under consideration. To be concrete, let us consider one of the standard model systems of the phase transition theory — the three-dimensional (3D) Ising model on the simple cubic lattice, with nearest-neighbor interactions.

According to universality, this model describes the critical properties of a wide range of systems of a different physical nature, including second order phase transitions in uniaxial magnetic systems and the liquid-vapor critical point. The statistical ensemble most commonly used (usually denoted as canonical ensemble) is defined by the partition function

$$Z_c(h)=\underset{\{s_i\}}{}\mathrm{exp}\left\{\beta \underset{ij}{}s_is_j+h\underset{i}{}s_i\right\},s_i=\pm 1,$$
(1)   
where the sum includes all the $`2^N`$ possible configurations of a total number of $`N`$ spins. This ensemble is perfectly natural for applications to magnetic phase transitions, where $`s_i=\pm 1`$ corresponds to physical spin at the lattice site $`i`$ pointing up or down, respectively. However, in the language of the lattice gas, $`s_i=\pm 1`$ corresponds to an occupied or unoccupied site, respectively: the total number of particles fluctuates. The canonical partition sum Eq. (1) in the spin language thus corresponds with the grand partition sum in the lattice gas language. To avoid confusion, we emphasize that in the rest of this paper we are using the word ‘canonical’ in the language of the Ising model: it means that the magnetization is allowed to fluctuate freely.

In contrast, many real experiments and simulations of liquid-vapor systems are performed with a fixed number of particles. Fixing the number of lattice-gas particles is equivalent, in the language of magnetic systems, to fixing the total magnetization $`M_{\mathrm{total}}_{i=1}^Ns_i`$ of the system, or, in other words, fixing the average magnetization per spin $`M\frac{1}{N}M_{\mathrm{total}}=\frac{1}{N}_{i=1}^Ns_i`$.

Thus we will be interested in the properties of the 3D Ising model in the fixed-magnetization ensemble,

$$Z_f(M)=\underset{\{s_i\}:_is_i=NM}{}\mathrm{exp}\left\{\beta \underset{ij}{}s_is_j\right\}=\underset{\{s_i\}}{}\delta _{NM,_is_i}\mathrm{exp}\left\{\beta \underset{ij}{}s_is_j\right\}.$$
(2)   
Note that the magnetic field $`h`$ is absent; it can only contribute a constant factor.

One of the main difficulties encountered in computer simulation studies of critical phenomena is the critical-slowing-down phenomenon. For a number of spin models in the canonical ensemble this problem has been largely overcome with the invention of cluster algorithms . Until recently, no similarly useful algorithm has become available for the fixed-magnetization ensemble. This situation has now changed by the development of a geometric cluster algorithm . We have used this algorithm extensively in this work to efficiently simulate systems of fixed magnetization at the critical point.

In the canonical Ising model described by Eq. (1), the magnetic field $`h`$ is an adjustable external parameter. In contrast, the magnetization is a fluctuating observable. For each configuration taken from the ensemble one can sample its magnetization per spin $`M=\frac{1}{N}_{i=1}^Ns_i`$. Having accumulated $`M`$ over a sufficiently large set of configurations, one can construct the probability distribution $`P(M)`$ , and determine various expectation values such as $`M`$, $`M^2`$. Examples of such probability distributions at the critical point are shown in Fig. 1.

On the other hand, for systems in the fixed-magnetization ensemble described by Eq. (2), the roles of $`h`$ and $`M`$ are interchanged: now $`M`$ is the adjustable parameter, and it is intuitively clear that there should be some way to define an observable, which we denote by $`\stackrel{~}{h}`$ (to avoid confusion with $`h`$ in Eq. (1)), that will correspond to the magnetic field. Thus $`\stackrel{~}{h}`$ will be a fluctuating quantity, that can be sampled on a microscopic level from configurations taken from the fixed-magnetization ensemble. In the limit $`N\mathrm{}`$ in both ensembles (such that the correlation length vanishes in comparison with the system size), the fluctuations in $`M`$ and $`\stackrel{~}{h}`$ become negligible, and the difference between $`h`$ and $`\stackrel{~}{h}`$ vanishes.

In Section II we discuss the definition of $`\stackrel{~}{h}`$ and its properties. In Section III we establish the relation between the function $`\stackrel{~}{h}(M)`$ in the fixed-$`M`$ ensemble and the probability distribution $`P(M)`$ in the canonical ensemble. We conclude with a discussion of the relation between the function $`\stackrel{~}{h}(M)`$ in the fixed-$`M`$ ensemble and the function $`h(M)`$ in the canonical ensemble, and with a summary of our main results.

## II The magnetic observable $`\stackrel{~}{h}(M)`$ for the fixed-$`M`$ ensemble

We will now describe a definition of $`\stackrel{~}{h}(M)`$ that is based on statistical analysis of the local environment of a given spin. By local environment we mean the set of neighbors with which this spin interacts. In our particular case of the Ising model with nearest-neighbor interactions on the simple cubic lattice, the local environment consists of 6 spins on the neighboring sites. The local environment has $`2^6`$ possible configurations which divide in 7 types: 0+6$``$ (zero spins up, 6 down), 1+5$``$, 2+4$``$, 3+3$``$, 4+2$``$, 5+1$``$, 6+0$``$ (six spins up, zero down). The simplest way to Monte Carlo sample $`\stackrel{~}{h}(M)`$ is on the basis of the symmetric case: 3+3$``$ . For every Monte Carlo configuration, go through all sites, and select all spins with the required 3+3$``$ local environment. Then compute the average $`s_0`$ of the selected spins, and define $`\stackrel{~}{h}`$ by

$$s_0=\mathrm{tanh}\stackrel{~}{h}.$$
(3)   
It is also possible to employ, instead of 3+3$``$, any other of the 7 types of local environment. In these non-symmetric cases the definition reads

$$s_0=\mathrm{tanh}\left(\stackrel{~}{h}+\beta \underset{i=1}{\overset{6}{}}s_{0i}\right),$$
(4)   
where the $`s_{0i}`$ are the nearest neighbors of $`s_0`$.

One easily notices that the definition is constructed in such a way that $`\stackrel{~}{h}`$ corresponds, on the mean field level, to the external field $`h`$ in Eq. (1).

Now it is interesting to see what are the results of Monte Carlo simulations for $`\stackrel{~}{h}(M)`$. As has already been demonstrated in , at the critical temperature $`\stackrel{~}{h}(M)`$ practically coincides with the relation $`h(M)`$ in the canonical ensemble as obtained by Monte Carlo simulations, provided $`M`$ is sufficiently large, so that the correlation length is sufficiently small in comparison with the system size and the finite-size effects are suppressed. At the same time, the striking feature of $`\stackrel{~}{h}(M)`$ for not-so-large $`M`$ is its nonmonotonic behavior. First $`\stackrel{~}{h}(M)`$ goes negative at small $`M`$, then it begins to grow, and finally assumes the usual behavior at larger $`M`$ . This is clearly seen in Fig. 2 (diamonds), which shows Monte Carlo results obtained by means of the the geometric cluster algorithm .

In the remaining part of the paper, we will give the explanation of this behavior (which turns out to be a peculiar kind of finite size effect), by establishing a close relation between $`\stackrel{~}{h}(M)`$ in the fixed-$`M`$ ensemble, and the probability distribution $`P(M)`$ in the canonical ensemble.

## III  Connection between $`\stackrel{~}{h}(M)`$ in the fixed-$`M`$ ensemble and the probability distribution of $`M`$ in the canonical ensemble

Considering the fixed-$`M`$ ensemble, Eq. (2), one notices that it can be obtained by taking the canonical ensemble (1), and cutting from it the subset satisfying the constraint, $`_is_i=NM`$. Within this subset we still have the usual Boltzmann probabilities $`\mathrm{exp}\{\beta _{ij}s_is_j\}`$ for individual configurations.

This makes it possible to establish a relation between $`\stackrel{~}{h}(M)`$ in the fixed-$`M`$ ensemble, and the properties of the system in the canonical ensemble. The definition of $`\stackrel{~}{h}(M)`$ described in Sect. II is equivalent to the following. Let us take the fixed-$`M`$ ensemble and concentrate our attention on one particular lattice site, and on the spin located there. Let us perform the following measurement. For every configuration consider the local environment of our selected site. If it is not 3+3$``$, do not measure anything for this configuration. If it is 3+3$``$, measure the selected spin $`s_0`$ and store it. Finally, find $`s_0`$, and use Eq. (3) to determine $`\stackrel{~}{h}`$.

One notices that, as long as we are performing a thought experiment, we need not care about the Monte Carlo statistics. We can just stick to one site and get the same $`s_0`$ without averaging over all sites, because they are equivalent.

Up to now we have distinguished between 7 types of local environments, such as 3+3$``$. Let us go a bit further and treat separately all $`2^6`$ possible local environments. In other words, the measurement of $`s_0`$ is now performed in an even smaller subset of the fixed-$`M`$ ensemble: also the six spins forming the local environment of $`s_0`$ are fixed. In the case that the predetermined local environment is of the 3+3$``$ type, it may seem that there is no interaction between $`s_0`$ and the remaining system of $`N7`$ spins. Nevertheless, the fixed-$`M`$ ensemble probabilities that $`s_0`$ is $`+1`$ or $`1`$, which we denote, respectively, by $`P_+`$ and $`P_{}`$, are not equal in general. These probabilities may still depend on the magnetization of the remaining system, which is coupled to $`s_0`$ by the overall magnetization constraint:

$$s_0+\underset{i=1}{\overset{6}{}}s_{0i}+\underset{iRS}{}s_i=\underset{i=1}{\overset{N}{}}s_iNM.$$
(5)   
The total magnetization of the system is thus expressed as the sum of three terms: the local spin $`s_0`$, the sum of its six neighbors and the magnetization of the remaining $`N7`$ spins denoted as $`_{iRS}s_i`$, where $`RS`$ stand for “remaining system”.

The conditional probabilities $`P_\pm `$ can be more explicitly written as

$$P_\pm =P(s_0=\pm 1|NM,s_{01}\mathrm{}s_{06})$$
(6)   
The two conditional arguments specify the total magnetization $`NM`$ and the states of the 6 neighbor spins. We now make the connection with the canonical probabilities $`P_c`$ which include the magnetization as an unconditional argument. We use the zero-field canonical probabilities i.e., $`h=0`$ in Eq. (1).

$$P_\pm =P_c^1(NM|s_{01}\mathrm{}s_{06})P_c(s_0=\pm 1,NM|s_{01}\mathrm{}s_{06})$$
(7)   
We may slightly rewrite this by substitution of the probability $`P_c`$ by $`\widehat{P}_c`$ which is equal but uses the magnetization of the $`N7`$ remaining spins as its second argument:

$$P_\pm =P_c^1(NM|s_{01}\mathrm{}s_{06})\widehat{P}_c(s_0=\pm 1,\underset{iRS}{}s_i=NM\underset{i=1}{\overset{6}{}}s_{0i}s_0|s_{01}\mathrm{}s_{06})$$
(8)   
Let us first consider the simplest case $`_{i=1}^6s_{0i}=0`$. Thus the canonical probability $`\widehat{P}_c`$ does not depend on its first argument, which can thus be skipped:

$$P_\pm =\frac{1}{2}P_c^1(NM|s_{01}\mathrm{}s_{06})\widehat{P}_c(\underset{iRS}{}s_i=NM1|s_{01}\mathrm{}s_{06})$$
(9)   
Therefore,

$$\frac{P_+}{P_{}}=\frac{\widehat{P}_c(\underset{iRS}{}s_i=NM1|s_{01}\mathrm{}s_{06})}{\widehat{P}_c(_{iRS}s_i=NM+1|s_{01}\mathrm{}s_{06})}.$$
(10)   
The condition $`s_{01}\mathrm{}s_{06}`$ in effect introduces a defect in the remaining system: an octahedron-shaped bubble with six spins at its vertices fixed, while the spin $`s_0`$ in the middle is decoupled and plays no role any more.

Obviously, the ratio (10) could be obtained by performing a usual canonical ensemble simulation of such a system with a defect, and measuring the probability distribution for its overall magnetization $`_{iRS}s_i`$. The value of the ratio (10) would then be given by the ratio of the heights of the two neighboring bins in the corresponding histogram.

In all cases of practical interest for the study of the scaling limit (sufficiently large systems, sufficiently small magnetization) the ratio (10) is close to 1. Otherwise a difference of one unit in the total magnetization would lead to a large change of probability: this would obviously be far from the scaling limit. Thus we always work with $`\stackrel{~}{h}1`$.

It is convenient to introduce a shorter notation $`P_{RS}(x)\widehat{P}_c(_{iRS}s_i=Nx|s_{01}\mathrm{}s_{06})`$ where $`x`$ may be read as the magnetization of the remaining system if the factor $`N`$ instead of $`N7`$ makes a negligible difference, i.e. for large systems. We have to keep in mind that the notation $`P_{RS}`$ refers to a system with a defect whose type is not explicitly shown. Again restricting ourselves to local environments of the type 3+3$``$, we obtain

$$s_0=\frac{P_+P_{}}{P_++P_{}}=\frac{P_{RS}(M\frac{1}{N})P_{RS}(M+\frac{1}{N})}{P_{RS}(M\frac{1}{N})+P_{RS}(M+\frac{1}{N})}.$$
(11)   
Thus we arrive at

$$s_0\frac{1}{N}\frac{1}{P_{RS}(M)}\frac{dP_{RS}(M)}{dM}.$$
(12)   
Also, due to $`\stackrel{~}{h}1`$, Eq. (3) reduces to

$$s_0=\stackrel{~}{h},$$
(13)   
and we get

$$\stackrel{~}{h}=\frac{1}{N}\frac{d}{dM}\mathrm{log}P_{RS}(M).$$
(14)   
Defining the effective potential $`V_{\mathrm{eff}}^{(RS)}(M)`$ (i.e., the Ginzburg – Landau fixed-$`M`$ free energy) of the present system (with defect) by

$$P_{RS}(M)\mathrm{exp}\{NV_{\mathrm{eff}}^{(RS)}(M)\},$$
(15)   
we get immediately

$$\stackrel{~}{h}=\frac{dV_{\mathrm{eff}}^{(RS)}(M)}{dM}.$$
(16)   
For large systems, the relative contribution of the defect is small, and thus $`\stackrel{~}{h}(M)`$ is well approximated by $`V_{\mathrm{eff}}(M)`$ for the finite system without a defect:

$$P(M)\mathrm{exp}\{NV_{\mathrm{eff}}(M)\},$$
(17)   
$$\stackrel{~}{h}=\frac{dV_{\mathrm{eff}}(M)}{dM}+\mathrm{}$$
(18)   
where the ellipsis stands for corrections vanishing at large $`N`$.

As is well known, for finite 3D Ising models in a cubic box with periodic boundary conditions, the distribution $`P(M)`$ has a double-peak structure  at the critical point. Thus $`V_{\mathrm{eff}}(M)`$ has a double-well shape, which immediately explains why $`\stackrel{~}{h}`$ goes negative for small values of $`M`$. In Fig. 2 we show the quantitative comparison of $`\stackrel{~}{h}(M)`$ (depicted by diamonds) and $`dV_{\mathrm{eff}}/dM`$ (solid line). One observes that the correspondence between the points and the line clearly improves with increasing lattice size. To extract $`V_{\mathrm{eff}}`$ from the Monte-Carlo-generated $`P(M)`$ (Fig. 1) we have exploited the fact, reported in , that for the system under consideration, $`P(M)`$ can be well approximated by the ansatz

$$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\},$$
(19)   
which applies to the finite-size regime, i.e. the finite size is small compared to the bulk correlation length. We have fitted the Monte Carlo generated $`P(M)`$ data accordingly, determined the parameters $`a`$ and $`c`$, and thus obtained $`dV_{\mathrm{eff}}(M)/dM`$ in a simple polynomial form.

It is also worth mentioning that the shape of $`P(M)`$ for a given geometry (in our case, a cubic box with periodic boundaries) is universal at the critical point. That is, the parameters $`a`$ and $`c`$ have well-defined scaling limits when the system size grows to infinity. These values, $`a=0.158(2)`$, $`c=0.776(2)`$, have been determined in  by making use of a special model in the 3D Ising universality class, which has almost no corrections to scaling . The corresponding scaling form of $`dV_{\mathrm{eff}}(M)/dM`$ is plotted by the dashed line in Fig. 2. One observes that deviations from scaling (between the solid and the dashed lines) go down with increasing size, as they should.

The results in Fig. 2 confirm the relation between the observable $`\stackrel{~}{h}(M)`$ as defined above in the fixed-$`M`$ ensemble, and the probability distribution $`P(M)`$ in the canonical ensemble. The remaining discrepancy (between the diamonds and the solid line in Fig. 2) is due to the “defect” discussed above. The question arises whether it is possible to modify our definition of $`\stackrel{~}{h}`$ in order to suppress this discrepancy. We have found that this is indeed the case. Up to this point, we restricted ourselves to symmetric local environments (3+3$``$) to define $`\stackrel{~}{h}`$ via Eq. (3). As has already been mentioned, using Eq. (4) one may use other types of local environments as well. In those cases the magnetization $`k_{i=1}^6s_i`$ enters the definition:

$$s_0=\mathrm{tanh}(\stackrel{~}{h}+\beta k).$$
(20)   
Following the same arguments as before, we decompose the system in the local spin $`s_0`$, its fixed neighbors, and the remaining system ($`RS`$). This leads to the following generalization of Eq. (11):

$`s_0`$ $`=`$ $`{\displaystyle \frac{e^{\beta k}P_{RS}(M\frac{k}{N}\frac{1}{N})e^{\beta k}P_{RS}(M\frac{k}{N}+\frac{1}{N})}{e^{\beta k}P_{RS}(M\frac{k}{N}\frac{1}{N})+e^{\beta k}P_{RS}(M\frac{k}{N}+\frac{1}{N})}}`$ (21)   
$``$ $`\mathrm{tanh}\beta k{\displaystyle \frac{1}{\mathrm{cosh}^2\beta k}}{\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{P_{RS}(M\frac{k}{N})}}{\displaystyle \frac{dP_{RS}(M\frac{k}{N})}{dM}}`$ (22)   
$``$ $`\mathrm{tanh}\left(\beta k{\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{P_{RS}(M)}}{\displaystyle \frac{dP_{RS}(M)}{dM}}\right).`$ (23)   
Thus we arrive once again at Eqs. (14-16). But we now have a different type of defect in the remaining system, and a shift of $`k/N`$ in the magnetization of the remaining system; we neglect the latter effect. Now it seems plausible that one can suppress the influence of the defect by averaging over all configurations of the defect, weighted with their natural occurrence probabilities. Such an averaging should more faithfully reproduce the characteristics of a system without a defect. The modified determination of $`\stackrel{~}{h}`$ is as follows. Sample configurations from the fixed-$`M`$ ensemble. For each spin determine its orientation ($`+`$ or $``$) and the type of its local environment (type 0, …, 6 for 0+6$``$, …, 6+0$``$). Accumulate these data by incrementing one out of 14 bins $`N_{q,+}`$, $`N_{q,}`$, where $`q=0\mathrm{}6`$ denotes the type of local environment, and $`+`$ or $``$ denotes the local spin. The resulting population numbers satisfy $`_{q=0}^6(N_{q,+}+N_{q,})=N`$. Then, for each $`q`$, find $`s_0_q=(N_{q,+}N_{q,})/(N_{q,+}+N_{q,})`$ and compute $`\stackrel{~}{h}_q`$ according to Eq. (4). Finally,

$$\stackrel{~}{h}_{\mathrm{improved}}=\frac{1}{N}\underset{q=0}{\overset{6}{}}\stackrel{~}{h}_q(N_{q,+}+N_{q,}).$$
(24)   
Applying this definition to our simulation data, we observe that, within the statistical accuracy, the discrepancy between $`\stackrel{~}{h}(M)`$ and $`dV_{\mathrm{eff}}(M)/dM`$ is indeed eliminated (Fig. 2).

## IV Discussion and conclusions

The relation (18) looks exactly the same as the standard relation between the field and magnetization in the canonical ensemble:

$$h=\frac{d\stackrel{~}{V}_{\mathrm{eff}}(M)}{dM}.$$
(25)   
The observed differences between the properties of $`\stackrel{~}{h}(M)`$ in the fixed-$`M`$ ensemble and $`h(M)`$ in the canonical ensemble, the most prominent of which is the nonmonotonic behavior of $`\stackrel{~}{h}(M)`$ instead of the monotonic behavior of $`h(M)`$, can be traced to the different definitions of the effective potential. The one that occurs in Eq. (18) is the fixed-$`M`$ free energy,

$$V_{\mathrm{eff}}(M)=(1/N)\mathrm{log}Z_f(M),$$
(26)   
while the one that enters Eq. (25) is defined via a Legendre transformation:

$$\stackrel{~}{V}_{\mathrm{eff}}(M)=(1/N)\mathrm{log}Z_c(h)+hM,$$
(27)   
where

$$M=M_h$$
(28)   
is the canonical average of the magnetization in an external field $`h`$. The partition functions $`Z_f(M)`$ and $`Z_c(h)`$ were defined in Section I. In a situation where fluctuations become negligible, the term $`hM`$ in $`\stackrel{~}{V}_{\mathrm{eff}}`$ cancels the field dependence of the Boltzmann weights. Then both definitions of the effective potential become equivalent, and both effective potentials approach the bulk form so that the difference between $`h`$ and $`\stackrel{~}{h}`$ vanishes.

In a finite system, due to fluctuations, $`\stackrel{~}{V}_{\mathrm{eff}}`$ differs from $`V_{\mathrm{eff}}`$. For instance, at the Ising critical point, the double-well form of $`V_{\mathrm{eff}}`$ is absent: $`\stackrel{~}{V}_{\mathrm{eff}}`$ has a single-well form.

Returning to Eq. (16), there is another finite-size effect: the difference between $`V_{\mathrm{eff}}^{(RS)}(M)`$ and $`V_{\mathrm{eff}}(M)`$ due to the presence of the defect. The relative contribution of the defect becomes small for large systems, and it is further suppressed by the improved definition of $`\stackrel{~}{h}`$, Eq. (24).

Thus for large systems and sufficiently high $`M`$, when the correlation length is small in comparison with the system size, the finite size effects are suppressed, the difference between $`V_{\mathrm{eff}}(M)`$ and the bulk effective potential disappears, and our definition of $`\stackrel{~}{h}(M)`$ reproduces the expected bulk behavior.

In conclusion, we have studied the critical three-dimensional Ising model in the fixed-magnetization ensemble, in a cubic geometry with periodic boundary conditions. This was done by means of the recently developed geometric cluster Monte Carlo algorithm. We have defined a magnetic field-like observable $`\stackrel{~}{h}`$ for this ensemble, studied its dependence on the magnetization $`M`$ and explained its counter-intuitive nonmonotonic behavior: $`\stackrel{~}{h}`$ first becomes negative and then positive with increasing $`M`$ (Fig. 2). We have provided a quantitative description of $`\stackrel{~}{h}(M)`$ by establishing a close relation with $`P(M)`$ — the probability distribution of the magnetization in the canonical ensemble. The nonmonotonic behavior of $`\stackrel{~}{h}(M)`$ can be understood as a manifestation of the same finite-size effect that is responsible for the double-peak shape of $`P(M)`$ at the critical point. Furthermore we have shown that, when fluctuations are negligible, our definition reduces to the standard canonical relation $`M(h)`$. Finally, we note that in the different context of the simulation of a system of particles whose number is fixed, a similar line of reasoning enables the determination of the chemical potential of the particles .

###### 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.