text:
# Phase transitions in a cluster molecular field approximation

## I Introduction

The molecular field approximation (MFA)  is the method of choice to gain a first, usually qualitatively correct, understanding of phase diagrams in cases where exact results are not available. However, quantitative predictions derived in the framework of this approximation compare rather poorly with results obtained by other, more refined, methods. For example, the molecular field approximation systematically overestimates the temperatures at which phase transitions take place. This is of course due to the fact that correlations are completely neglected in this treatment. Furthermore, the critical exponents, describing the behaviour of critical quantities in the vicinity of a critical point, usually do not agree with the molecular field critical exponents. Above the upper critical dimension, however, the molecular field predictions of critical exponents become exact.

The shortcomings of the original MFA led Bethe  (see also Ref.  and ) to consider a more general treatment. In the Weiss theory of ferromagnetism, a single spin is considered in the molecular field which is caused by all the other spins. This is different in the Bethe approach where a small cluster of spins, consisting of the central spin and its nearest neighbours, is treated exactly. The molecular field acts exclusively on the edges of the cluster. Thus, correlations between a certain number of spins are included into the calculation. This is a significant improvement over simple molecular field theory, yielding for example critical temperatures which are much closer to the correct values than those obtained from the Weiss theory.

The general idea of inserting a cluster of spins into a homogeneously magnetised medium led to several applications in the past. Müller-Krumbhaar and Binder, for example, proposed a Monte Carlo method where a molecular field acts only on the surface spins of a cluster. This field has to be determined in a self-consistent way during the course of the simulation . In a series of papers, Suzuki and coworkers established a general theoretical framework for extracting non-classical critical exponents from cluster-mean-field approximations . They considered a systematic series of clusters and showed that the true critical exponents can be obtained from the amplitudes of mean-field-type singularities. Various other generalisation schemes of the one-spin molecular field approximation have been discussed in the literature .

In this work we present a systematic way of improving the molecular field approximation by considering spin clusters of variable sizes embedded into a homogeneously magnetised background. Hereby, the interactions within the clusters are taken into account exactly. To this end we compute the density of states as a function of the internal energy and of the interaction energy with a molecular field of unit strength. This is a hard task, but when it has been completed everything else does no longer demand time-consuming calculations. The actually applied molecular field itself is determined by imposing a condition of homogeneity in the central part of the cluster. This approach is used to evaluate various temperature dependent quantities for Ising models defined on different two- and three-dimensional lattices. A two-dimensional Ising model with nearest and next-nearest neighbour couplings has also been studied. Results obtained for three-state Potts models on square and triangular lattices are discussed, too.

With these data we demonstrate that critical quantities (critical temperatures as well as critical exponents) can be extracted reliably from the cluster molecular field approximation. This is already achieved by inserting very small clusters containing a small number of spins (the largest cluster studied is formed by 181 spins) into the molecular field. Besides analysing effective exponents we also show that the cluster molecular field data exhibit a remarkable finite size scaling behaviour which can be exploited for the determination of the true non-classical values of the critical exponents.

The paper is organised in the following way. In the next Section we present our method of systematically improving the molecular field approximation. We also discuss how to obtain various quantities from this data. In Section III the molecular field finite size scaling theory is discussed. We thereby use both the reduced temperature and the reduced energy as possible scaling variables. Section IV presents our data. We analyse them in different ways and show that critical quantities can be obtained from small clusters in a molecular field. Finally, in the last Section we present our conclusions.

## II Clusters in a molecular field

In the following we propose a systematic scheme for improving the molecular field approximation. This scheme combines the exact computation of the density of states of small clusters with a molecular field acting solely on the borders of the clusters. In Section IV this method is applied to the determination of non-classical critical quantities of various Ising and three-state Potts models. Before discussing our approach in detail we will briefly reiterate the reasoning leading to the molecular field self consistency equation for these two models.

### A The molecular field approximation

In zero field the Ising Hamiltonian reads

$$H_I=J_I\underset{i,j}{}\sigma _i\sigma _j$$
(1)
where $`\sigma _i=\pm 1`$ is an Ising spin on the lattice. Angular brackets denote sums over nearest neighbour pairs and $`J_I`$ is the exchange constant. The Hamiltonian for one selected spin $`i`$, given the spin values $`\sigma _{i+\delta }`$ of its $`z`$ neighbours, is

$$\stackrel{~}{H}_I=J_I\sigma _i\underset{\delta =1}{\overset{z}{}}\sigma _{i+\delta }.$$
(2)
If the spin $`\sigma _i`$ is considered and the rest of the lattice is replaced by a homogeneously magnetised medium with magnetisation $`m`$ then the spin variables $`\sigma _{i+\delta }`$ are replaced by the expectation value $`\sigma _{i+\delta }=m`$ and $`\stackrel{~}{H}_I`$ becomes $`H_{I,MFA}(\sigma _i)=zJ_Im\sigma _i`$ .

At the temperature T the expectation value of $`\sigma _i`$ in the molecular field $`zJ_Im`$ is:

$$\sigma _i=\frac{Z_+Z_{}}{Z}=\mathrm{tanh}(zJ_Im/k_BT)$$
(3)
where $`Z_+=\mathrm{exp}(zJ_Im/k_BT)`$ , $`Z_{}=\mathrm{exp}(zJ_Im/k_BT)`$ and $`Z=Z_++Z_{}`$ is the partition function for the two states of a single spin.

The condition $`\sigma _i=m`$ yields the self consistency equation which can be solved for the temperature

$$k_BT(m)/J_I=\frac{zm}{\text{artanh}(m)}.$$
(4)
Eq. (4) can be inverted or solved graphically if necessary. From the limit $`m0`$ of (4) one finds the MFA critical temperature $`k_BT_{c,MFA}/J_I=z`$ and for small but finite $`m`$ one obtains

$$k_BT(m)/J_I=k_BT_{c,MFA}/J_I\frac{z}{3}m^{1/\beta _{MFA}}+\mathrm{}$$
(5)
with the critical exponent $`\beta _{MFA}=1/2`$.

The derivation of the corresponding self consistency equation for the three-state Potts model  defined by the Hamiltonian (the spins $`\sigma _i`$ taking the values 1, 2, 3)

$$H_P=J_P\underset{i,j}{}\delta (\sigma _i,\sigma _j)$$
(6)
with $`\delta (\sigma _i,\sigma _j)=1`$ if $`\sigma _i=\sigma _j`$ and zero otherwise is somehow more involved, so we only quote the result

$$m=\frac{\mathrm{exp}(zJ_Pm/k_BT)1}{\mathrm{exp}(zJ_Pm/k_BT)+2}.$$
(7)
Inspection of the molecular field free energy reveals that at the temperature

$$k_BT_c/J_P=\frac{z}{4\mathrm{ln}(2)}$$
(8)
a first order transition takes place where the order parameter $`m`$ jumps from 0 to the value 1/2. This discontinuous behaviour in the MFA is in contrast to the fact that the temperature driven phase transition is continuous for the three-state Potts model in two dimensions.

These highly familiar results have been repeated also for the purpose of analysing the main defect of the approach. It is the lack of any correlation between the selected spin $`\sigma _i`$ and all the other spins in the lattice. This can be partially repaired if a cluster of $`N`$ spins is inserted into the homogeneously magnetised medium. In our approach, as discussed in the following, interactions within such clusters are taken into account exactly whereas the effect of the magnetised background is treated in the ordinary Weiss approximation. Therefore the correlations of the central spin with a number of neighbouring spins are included into the calculation thus leading to greatly improved results.

### B Cluster molecular field approximation

In this Section we present our method to calculate physical quantities within an extended molecular field approximation based on the idea of Bethe.

Consider a classical spin system on a lattice defined by the Hamiltonian

$$H=\underset{i,j}{}I(\sigma _i,\sigma _j)$$
(9)
where $`I(\sigma _i,\sigma _j)`$ is the interaction energy of the two neighbouring spins $`\sigma _i`$ and $`\sigma _j`$. For the Ising model this is just the negative product of the spin variables, i.e. $`I_I(\sigma _i,\sigma _j)=J_I\sigma _i\sigma _j`$, for the $`q`$-state Potts model ($`\sigma _i=1,\mathrm{},q`$) the interaction energy is given by $`I_P(\sigma _i,\sigma _j)=J_P\delta (\sigma _i,\sigma _j)`$. Suppose one has constructed a cluster of size $`N`$ by some method (see below) which is embedded in the original lattice. For this cluster the Hamiltonian (9) is replaced by the molecular field Hamiltonian

$$H_N=\underset{i,j}{}I(\sigma _i,\sigma _j)+\mu \underset{\{i,k\}}{}I(\sigma _i,\sigma _k)\delta (\sigma _k,1)$$
(10)
where $`i,j`$ indicates a summation over nearest neighbour pairs within the cluster. The second term in (10) represents the interaction of the cluster with the homogeneously magnetised background. The sum runs over all pairs of nearest neighbour spins where the spin $`\sigma _i`$ is from inside the cluster and its partner $`\sigma _k`$ is from outside. This is denoted by the curly brackets. For convenience the $`\sigma _k`$ are set equal to 1 and their common expectation value $`\sigma _k=\mu `$ is placed outside the summation. The Hamiltonian (10) of the extended molecular field approximation treats the statistical mechanics of the spin system within the cluster exactly and the interaction of the cluster with the background is modelled by a molecular field that has to be determined in a self consistent way. In order to calculate physical quantities in the cluster molecular field approximation it is advantageous to work out the density of states

$$\mathrm{\Omega }_{\sigma _0,\sigma _1}(E,B)=\underset{\mathrm{\Gamma }_N(\sigma _0,\sigma _1)}{}\delta (H_C,E)\delta (H_B,B),$$
(11)
with the internal energy of the cluster

$$H_C=\underset{i,j}{}I(\sigma _i,\sigma _j)$$
(12)
and the interaction energy

$$H_B=\underset{\{i,k\}}{}I(\sigma _i,\sigma _k)\delta (\sigma _k,1)$$
(13)
of the cluster with the spins outside. The summation in (11) runs over all $`\mathrm{\Gamma }_N(\sigma _0,\sigma _1)`$ states of the cluster with fixed spin value $`\sigma _0`$ of the central spin and $`\sigma _1`$ of one of its $`z`$ equivalent neighbours. For classical discrete spin models the density of states is just the number of configurations with specified values $`E`$, $`B`$, $`\sigma _0`$ and $`\sigma _1`$.

To work out physical quantities within this approximation scheme the density of states is evaluated exactly so that these quantities can be computed with arbitrary accuracy. The technical details of computing the density of states (11) for a given cluster size $`N`$ are summarised in the appendix. The size of the clusters which can be handled exactly is limited by the available computational capacities. To treat bigger clusters a recently proposed, highly efficient Monte Carlo algorithm to compute the density of states could be used . However, a Monte Carlo calculation has to be performed with great care to achieve the desired accuracy.

Clusters of increasing particle numbers are selected in two different ways. A possible systematic series of clusters is constructed by choosing a spin to be the central spin and successively including shells of nearest neighbours. Hence, the first cluster of $`z+1`$ spins is made up of the central spin and its $`z`$ nearest neighbours. The next cluster consists of the first cluster and the shell of nearest neighbours of the border spins of the first cluster. By this way a new cluster is made up of the previous one together with the nearest neighbours of the boundary sites of the original cluster. This method has the advantage that assemblies of the same shape are build up. Fig. 1 displays the cluster of 25 spins emerging from this method for a two-dimensional square lattice. For the topology of the square lattice a cluster consists of $`L`$ rows of length $`L`$ perpendicular to the bisector of a quadrant of the plane and of $`L1`$ short rows of length $`L1`$ between the longer ones. Hence the number of spins in a cluster is given by $`N=L^2+(L1)^2`$. The construction of successive clusters leads to a series of clusters with 5, 13, 25, 41, 61, 85, 113, 145, and 181 spins for the two-dimensional square lattice which can be handled with the available computational capacities. For the triangular lattice clusters of 7, 19, 37, 61, 91, and 127 spins emerge. In this geometry the basic row and the central row have the lengths $`L`$ and $`2L1`$, respectively. However to avoid additional computational effort due to the length of the central row for the treatment of the triangular lattice we alternatively use the clusters of the square lattice by inserting additional bonds between the sites of the rows of lengths $`L`$ and $`L1`$ so that all sites have the coordination number 6 of the triangular topology. In three dimensions the equivalent procedure would be restricted to only one, or at most two system sizes by the available storage capacities. This imposes a less systematic approach. The clusters are constructed by adding successive shells of neighbours of the central spin. With the nearest neighbours of the central spin one obtains a cluster of 7 spins, with the next nearest neighbours the size is 19, then 27, 33, 57, and so on. The cluster of 57 spins that consists of the central spin and includes all spins up to the 5th neighbours was the largest one which could be handled in three dimensions. This way of selecting clusters for the three-dimensional case has the disadvantage that the various assemblies do not have the same geometry.

Once the density of states is known physical quantities are easily evaluated. For example the expectation value of the central spin is given by

$$\sigma _0=\underset{\sigma _0,\sigma _1}{}\sigma _0Z_{\sigma _0,\sigma _1}/Z$$
(14)
with the restricted partition functions

$$Z_{\sigma _0,\sigma _1}=\underset{E,B}{}\mathrm{\Omega }_{\sigma _0,\sigma _1}(E,B)\mathrm{exp}\left((E+\mu B)/k_BT\right)$$
(15)
and the total partition function

$$Z=\underset{\sigma _0,\sigma _1}{}Z_{\sigma _0,\sigma _1}.$$
(16)
Similarly the expectation value of the spin $`\sigma _1`$ is computed from (11). The unknown homogeneous magnetisation $`\mu `$ of the sites outside the cluster must be determined self consistently. When the expectation value $`\sigma _0`$ of the central spin is calculated as a function of $`\mu `$ and $`T`$ the temperature might be adjusted such that $`\sigma _0`$ is equal to $`\mu `$:

$$\sigma _0(T)=\mu (T).$$
(17)
Thereby the rigid external magnetisation is moved further away from the central spin and an improvement of the plain molecular field approximation is observed . This is the usual condition to determine the molecular field discussed in the literature. Nevertheless only the central spin is directly related to the rigid background magnetisation through the self-consistency equation so that the expectation values of the central spin and one of its neighbour spins are different. In our approach the condition

$$\sigma _0(T)=\sigma _1(T)$$
(18)
is used instead. This condition, which implicitly determines $`\mu `$ as a function of $`T`$, ensures that the spin in the centre of the cluster and its neighbouring spins have the same expectation value so that at least in the central part of the cluster the equivalence of the spins is implemented. The spins in the centre of the cluster are not directly coupled to the rigid background magnetisation. In the extended molecular field approximation based on condition (18), fluctuations of the expectation value of a spin near the border of the cluster are suppressed due to the direct coupling to the homogeneous background magnetisation whereas the central spin can develop comparatively large fluctuations. The behaviour of the central spin in the approximation used is consequently most closely related to the corresponding behaviour of a spin in the infinite system. The same holds for the interaction energies of bonds in the centre compared to those of bonds near the boundary of the cluster. These aspects suggest that the calculation of physical quantities should be restricted to the central part of the cluster where the condition of homogeneity is rather well fulfilled.

After solving the implicit equation (18) for the molecular field all physical quantities can be evaluated using the density of states (11). The order parameter of the system is determined from $`\sigma _0`$. Note that in the limit $`T0`$ the interaction of the border spins of the cluster with the homogeneously magnetised background forces all spin variables to be in the spin state $`\sigma =1`$ of the spins outside the cluster. For the Ising model the order parameter is given by

$$m_I(T)=\sigma _0(T)$$
(19)
and for the $`q`$-state Potts model one has

$$m_P(T)=\frac{q\sigma _0(T)1}{q1}.$$
(20)
In the cluster molecular field approach the order parameter is exactly zero at high temperatures. When decreasing the temperature the order parameter of the Ising model exhibits a sudden onset at the temperature $`T_{c,N}`$ indicating a transition to a magnetised phase. The abrupt appearance of a finite order parameter defines the pseudo-critical temperature of the finite system. Fig. 2 illustrates this temperature dependent behaviour of the order parameter in the case of the square lattice Ising model. It is instructive to compare the pseudo-critical temperatures resulting from equations (17) and (18). For a cluster of $`N=145`$ spins defined on the square lattice our approach (18) leads to the value $`k_BT_{c,N}/J_I=2.37096`$ ($`T_{c,N}/T_c=1.045`$, where $`k_BT_c/J_I=2/\mathrm{ln}(1+\sqrt{2})`$ is the critical temperature of the two-dimensional Ising model defined on the square lattice), which should be compared to the value 2.57524 ($`T_{c,N}/T_c=1.135`$) obtained when using the self-consistency equation (17) . A similar remarkable improvement is also achieved for three-dimensional clusters: for $`N=57`$ spins our pseudo-critical temperature is $`k_BT_{c,N}/J_I=4.6627`$ ($`T_{c,N}/T_c=1.034`$), much closer to the value $`4.5115`$ of the infinite system than the value $`5.0018`$ ($`T_{c,N}/T_c=1.109`$) resulting from Eq. (17) . The pseudo-critical temperatures of the Ising model for various lattice topologies are summarised in Table I and Table II for the two-dimensional and the three-dimensional systems, respectively.

In Fig. 3 we show the order parameter for the three-state Potts model defined on a square lattice as a function of temperature for a series of cluster sizes. For all clusters a jump in the order parameter is observed, thus disclosing the discontinuous character of the phase transition in the MFA. This jump gets smaller when the number of cluster spins is increased reflecting the softening of the first order character of the phase transition with increasing cluster size. In the limit $`N\mathrm{}`$ the jump will disappear completely in accordance with the continuous phase transition taking place in the bulk three-state Potts model. The temperature at which the jump occurs is deduced from the behaviour of the free energy. The values of these transition temperatures are displayed in Table III for the square and the triangular lattice.

We estimate the internal energy from the bonds of the central spin with its $`z`$ neighbour spins which yields

$$u(T)=\frac{z}{2}\underset{\sigma _0,\sigma _1}{}I(\sigma _0,\sigma _1)Z_{\sigma _0,\sigma _1}/Z.$$
(21)
From the internal energy one now gets the specific heat

$$c(T)=\frac{du(T)}{dT}.$$
(22)
The specific heat is also related to the entropy of the central spin through the relation

$$c(T)=T\frac{ds(T)}{dT}.$$
(23)
This allows the computation of the free energy

$$f(T)=u(T)Ts(T)$$
(24)
from the knowledge of the internal energy. Solving equation (23) for the entropy $`s(T)`$ one obtains the free energy

$$f(T)=u(T)T_{T_0}^T𝑑\tau \frac{c(\tau )}{\tau }$$
(25)
with $`T_0`$ being some arbitrarily chosen temperature representing the constant of integration. The free energy of the Potts model for the cluster with 13 spins in the vicinity of the transition point is shown in the inset of Fig. 3. The solution of (20) that minimises the free energy is physically stable. This leads to the kink of the free energy at the transition temperature again indicating the discontinuity of the phase transition for finite $`N`$.

Next we consider the response of the central spin to an external field that acts exclusively on the central spin itself. In order to evaluate this zero field susceptibility one has to differentiate the magnetisation with respect to a magnetic field that acts only on the central spin variable $`\sigma _0`$:

$$\chi (T)=\underset{h0}{lim}\frac{m(T,h)}{h}.$$
(26)
To obtain an analytic expression for the susceptibility of the Ising model we add a magnetic field term to the Hamiltonian:

$$H_IH_Ih\sigma _0.$$
(27)
When evaluating the susceptibility one has to bear in mind that one gets a field-dependent molecular field $`\mu =\mu (T,h)`$ due to this additional external field $`h`$. This has the consequence that (26) depends on the derivative of $`\mu `$ with respect to the magnetic field $`h`$. However, the order parameter $`m`$ is given by the two equivalent equations $`m=\sigma _0`$ and $`m=\sigma _1`$ so that one of them can be used to eliminate $`\mu /h`$ from (26). Differentiating $`\sigma _0`$ with respect to $`h`$ one obtains the expression

$$\chi (T,h)=\frac{1}{k_BT}\left[\sigma _0^2\sigma _0^2+\frac{\mu }{h}(\sigma _0B\sigma _0B)\right].$$
(28)
On the other hand differentiation of $`\sigma _1`$ yields

$$\chi (T,h)=\frac{1}{k_BT}\left[\sigma _0\sigma _1\sigma _0^2+\frac{\mu }{h}(\sigma _0B\sigma _1B)\right].$$
(29)
Combining these two relations in order to eliminate the derivative of the molecular field results in the desired expression of the zero-field susceptibility for temperatures below $`T_{c,N}`$:

$`\chi (T)`$ $`=`$ $`{\displaystyle \frac{1}{k_BT}}[\sigma _0^2\sigma _0^2`$ (31)
$`{\displaystyle \frac{(\sigma _0B\sigma _0B)(\sigma _0^2\sigma _0\sigma _1)}{\sigma _0B\sigma _1B}}].`$
The first two terms in the square brackets represent the ordinary fluctuation formula for the susceptibility. The additional third term arises from the coupling to the molecular field. Its denominator vanishes in the limit $`TT_{c,N}`$ yielding a diverging susceptibility for $`TT_{c,N}`$ as expected for the molecular field approach. Fig. 4 displays the zero field susceptibility of the two-dimensional Ising model defined on the square lattice for different cluster sizes.

For the $`q`$-state Potts model the magnetic field acts on the central spin only if it is in the spin state $`\sigma `$ of the sites outside the cluster:

$$H_PH_Ph\delta (\sigma _0,\sigma ).$$
(32)
The susceptibility of the Potts model for temperatures below the pseudo-critical point of the finite cluster is then given by

$`\chi (T)`$ $`=`$ $`{\displaystyle \frac{q}{q1}}{\displaystyle \frac{1}{k_BT}}[\sigma _0\delta (\sigma _0,\sigma )\sigma _0\delta (\sigma _0,\sigma )`$ (34)
$`{\displaystyle \frac{(\sigma _0B\sigma _0B)(\sigma _0\delta (\sigma _0,\sigma )\sigma _1\delta (\sigma _0,\sigma ))}{\sigma _0B\sigma _1B}}].`$
This expression for the Potts model has the same structure as the corresponding expression (31) of the Ising model. However the susceptibility of the Potts model does not diverge as the molecular field displays a jump at the transition temperature.

## III Molecular field finite size scaling theory

In this Section we present a phenomenological finite size scaling theory for our cluster molecular field approximation. We consider both the reduced temperature and the reduced energy as scaling variables and discuss the relative size of the corrections to scaling. Our data are then analysed in the following Section using the framework of the molecular field finite size scaling to be exposed in the following subsection.

### A Finite size scaling for the temperature picture

The spontaneous magnetisation of the cluster molecular field approximation displays the characteristics of spontaneous symmetry breaking. This allows an unambiguous definition of the pseudo-critical temperature $`T_{c,L}`$ of the finite system of linear extension $`L`$ taken to be $`N^{1/d}`$ where $`d`$ is the dimensionality of the lattice. The spontaneous magnetisation near this pseudo-critical temperature then behaves like

$$m_{sp,L}(t^{})=A_L(t^{})^\beta ^{}$$
(35)
where the reduced temperature $`t^{}`$ is defined as

$$t^{}=\frac{TT_{c,L}}{T_{c,L}}$$
(36)
and where $`\beta ^{}=0.5`$ is the molecular field critical exponent. For increasing clusters the influence of the molecular field which acts only on the boundary becomes weaker so that the approximation eventually becomes exact in the limit $`L\mathrm{}`$. In the infinite system the spontaneous magnetisation vanishes like

$$m_{sp}(t)=A(t)^\beta $$
(37)
for small reduced temperatures $`t`$ which is related to the critical temperature $`T_c=T_{c,\mathrm{}}`$ of the infinite system by

$$t=\frac{TT_c}{T_c}.$$
(38)
The fractional shift $`\tau _L`$ is given by

$$\tau _L=\frac{T_cT_{c,L}}{T_c}DL^\lambda $$
(39)
with the shift exponent $`\lambda =1/\nu `$ where $`\nu `$ is the critical exponent of the correlation length. Equation (39) and (36) are combined to yield

$$t^{}=t+\tau _L+O(\tau _L^2)$$
(40)
for large $`L`$. Equation (39) can be used to obtain an approximation of the transition temperature in the thermodynamic limit. This is exploited in Section IV.

For sufficiently small $`t^{}`$ we assume that the spontaneous magnetisation of the finite system of extension $`L`$ is given by

$$m_{sp,L}(t^{})=L^p\varphi (t^{}L^{1/\nu }).$$
(41)
The function $`\varphi (x)`$ is universal showing no further $`L`$-dependence. The exponent $`p`$ in (41) is obtained by taking the limit $`L\mathrm{}`$ for fixed temperature $`TT_c`$. In this limit equation (41) reduces to (37) yielding

$$\varphi (x)x^\beta \text{ as }x\mathrm{}.$$
(42)
The residual $`L`$ dependence is cancelled out if $`p=\beta /\nu `$. Hence the finite size scaling relation

$$L^{\beta /\nu }m_{sp,L}(t^{})=\varphi (t^{}L^{1/\nu })$$
(43)
follows.

At this point we have to mention that an unambiguous definition of the linear extension of our finite clusters is not possible. This crucial fact was already pointed out in . For our investigations we choose $`L=N^{1/d}`$ but other definitions are also possible. For the three-dimensional systems our choice is especially problematic as the clusters do not have the same geometry. We will come back to this point in Section IV when discussing our numerical data.

From equation (35) follows that the universal function $`\varphi `$ behaves like a square root for small scaling variables, i.e.

$$\varphi (t^{}L^{1/\nu })=A_\varphi \sqrt{t^{}L^{1/\nu }},$$
(44)
for $`t^{}0`$ and fixed $`L`$. Taking the limit $`t^{}0`$ and comparing with relation (41) the finite size scaling law

$$A_L=A_\varphi L^{\frac{\beta ^{}\beta }{\nu }}$$
(45)
for the amplitudes of the magnetisation of the finite systems is derived.

In an analogous way we obtain a finite size scaling relation for the susceptibility. In the molecular field approach the susceptibility diverges with the molecular field critical exponent $`\gamma ^{}=1`$. We then have the scaling relation

$$L^{\gamma /\nu }\chi _L(t^{})=\psi (t^{}L^{1/\nu }),$$
(46)
where for small scaling variables the universal function $`\psi `$ is given by

$$\psi (t^{}L^{1/\nu })=A_\psi (t^{}L^{1/\nu })^1.$$
(47)
One also derives, in complete analogy with Eq. (45), a finite size scaling law involving the finite size amplitudes of the susceptibility. Note that this relation was already used in references  to get estimates for the critical exponent $`\gamma `$. All these relations show that the critical exponents of the infinite system are encoded in the amplitudes of the corresponding quantity of the finite system. It is precisely this fact which will allow the determination of critical exponents from the cluster molecular field approximation. It is also worth noting that in the present approach the critical temperature of the infinite system is not needed for the determination of critical exponents.

### B Finite size scaling for the energy picture

The internal energy of a system is given by the expectation value (21). The expression $`u_L=u_L(T)`$ can be inverted in order to obtain the temperature as a function of the internal energy of the system. All physical quantities can then be expressed as functions of the internal energy. In the following subsections this transformation from the temperature picture to the energy picture is discussed for systems displaying either an algebraically or a logarithmically diverging specific heat in the thermodynamic limit.

#### 1 Algebraically diverging specific heat

Consider the spontaneous magnetisation near the pseudo-critical point of the finite system:

$$m_{sp,L}(t^{})=A_\varphi L^{\frac{\beta ^{}\beta }{\nu }}\sqrt{t^{}}(1+B_Lt^{}+\mathrm{})$$
(48)
where a correction term to the leading behaviour with the $`L`$-dependent amplitude $`B_L`$ is explicitly taken into account. From finite size scaling theory one obtains a similar expression for the (intensive) zero-field specific heat near $`t^{}=0`$ for negative $`t^{}`$:

$$c_{h=0,L}(t^{})=aL^{\frac{\alpha \alpha ^{}}{\nu }}(1+2b_Lt^{}+\mathrm{})$$
(49)
with an $`L`$-dependent amplitude of the linear correction term. Here it is assumed that the singularity of the specific heat of the infinite system is characterised by $`\alpha 0`$. The case of a logarithmic divergence is deferred to the next subsection. Note also that the specific heat in molecular field approximation does not diverge. It exhibits a jump singularity, hence one has the exponent $`\alpha ^{}=0`$. At the pseudo-critical point of the finite system, i.e. $`t^{}=0`$, equation (49) reduces to

$$C_L=c_{h=0,L}(0)=aL^{\frac{\alpha \alpha ^{}}{\nu }}.$$
(50)
This relation can be used to estimate the true critical exponent $`\alpha `$ from the maxima of the specific heat of the various finite systems.

As

$$c_{h=0,L}=\frac{du(T)}{dT}=\frac{du}{dt^{}}\frac{dt^{}}{dT}=\frac{1}{T_{c,L}}\frac{du}{dt^{}}$$
(51)
integration of (49) gives

$$\epsilon ^{}(t^{})=u_L(t^{})u_{c,L}$$
(52)
with $`u_{c,L}`$ being the constant of integration. This constant is just the pseudo-critical energy of the finite system. The expression for $`\epsilon ^{}(t^{})`$ can be inverted for small $`t^{}`$ to give

$$t^{}=\frac{1}{T_{c,L}a}L^{\alpha /\nu }\epsilon ^{}\frac{b_L}{T_{c,L}^2a^2}L^{2\alpha /\nu }\epsilon ^2+\mathrm{}$$
(53)
in the asymptotic limit $`\epsilon ^{}0`$. Plugging this result in equation (48) one finally obtains

$`m_{sp,L}(\epsilon ^{})`$ $`=`$ $`A_\varphi L^{\frac{\beta ^{}\beta _\epsilon }{\nu _\epsilon }}\sqrt{{\displaystyle \frac{\epsilon ^{}}{aT_{c,L}}}}`$ (55)
$`\left(1+L^{\alpha /\nu }{\displaystyle \frac{1}{aT_{c,L}}}(B_Lb_L)\epsilon ^{}+\mathrm{}\right)`$
with the renormalised exponents $`\beta _\epsilon =\frac{\beta }{1\alpha }`$ and $`\nu _\epsilon =\frac{\nu }{1\alpha }`$. These renormalised critical exponents describe the singular behaviour of the various thermostatic quantities of the infinite system in the energy picture . The reformulation also yields the finite size scaling relation

$$L^{\beta _\epsilon /\nu _\epsilon }m_{sp,L}(\epsilon ^{})=\stackrel{~}{\varphi }\left(\frac{\epsilon ^{}L^{1/\nu _\epsilon }}{T_{c,L}}\right)$$
(56)
for the energy picture with a scaling variable that contains the transition temperature of the finite system. In a completely analogous way the reformulation can be done for other physical quantities like the susceptibility. For large systems the pseudo-critical temperature $`T_{c,L}`$ is close to the critical temperature of the infinite system so that is can be replaced by $`T_c`$. However, as we are dealing with very small clusters the additional $`L`$-dependence in the scaling variable due to $`T_{c,L}`$ can not be neglected. Similar to equation (39) one obtains the shift relation for the pseudo-critical energy:

$$e_L=\epsilon _{c,\mathrm{}}\epsilon _{c,L}D_\epsilon L^{1/\nu _\epsilon }.$$
(57)
The proposed scaling variable for the energy picture and the corresponding scaling relations are tested in Section IV.

In equation (55) the amplitude of the linear correction term contains the amplitudes $`B_L`$ and $`b_L`$ of the temperature picture. Moreover a prefactor $`L^{\alpha /\nu }`$ appears which will lead to a suppression of the leading correction term compared to the corresponding expression (48) in the temperature picture if the critical exponent $`\alpha `$ is positive. This suppression is determined by the prefactor together with the $`L`$-dependence of the amplitudes $`B_L`$ and $`b_L`$. If the amplitudes grow strongly enough with the system size the suppression becomes weaker for increasing $`L`$. Note however that a reduction of the linear corrections in the energy picture will eventually always take place even if this reduction becomes weaker for larger cluster sizes. For a system with a negative exponent the correction term in the energy picture is blown up relative to the temperature picture. The investigation of this interesting consequence in the $`3d`$ Heisenberg model as an example of a system with $`\alpha <0`$ is left to future studies.

#### 2 Logarithmically diverging specific heat

With slight modifications the above discussed transfer to the energy picture is carried out for the case of a logarithmically diverging specific heat $`c=A\mathrm{ln}t`$ of the infinite system. A prominent example of such a system is the two-dimensional Ising model. In this case the amplitude $`C_L`$ (i.e. the maximum) of the finite system specific heat varies as

$$C_L=A\mathrm{ln}L^{1/\nu }+\xi =A^{}\mathrm{ln}L+\xi .$$
(58)
Fig. 5 shows equation (58) for the $`2d`$ Ising model in the extended molecular field approximation. The depicted linear behaviour is characterised by the slope $`A^{}=0.562`$ and the intercept $`\xi =1.76`$ with the vertical axis. For sufficiently large system sizes the $`L`$-independent constant $`\xi `$ can be neglected relative to the increasing contribution from the logarithmic term. As relatively small systems are considered in the molecular field approximation discussed here this constant is explicitly taken into consideration. The arguments used to replace the variable $`t^{}`$ by the internal energy still apply and give rise to the relation

$`m_{sp,L}(\epsilon ^{})`$ $`=`$ $`A_\varphi L^{\frac{\beta ^{}\beta _\epsilon }{\nu _\epsilon }}\sqrt{{\displaystyle \frac{\epsilon ^{}}{(A^{}\mathrm{ln}L+\xi )T_{c,L}}}}`$ (60)
$`\left(1+{\displaystyle \frac{1}{(A^{}\mathrm{ln}L+\xi )T_{c,L}}}(B_Lb_L)\epsilon ^{}+\mathrm{}\right).`$
The linear correction term is again diminished by an additional $`L`$-dependent prefactor as already observed for systems with $`\alpha >0`$. This form of the magnetisation suggests the finite size scaling law

$$L^{\beta _\epsilon /\nu _\epsilon }m_{sp,L}(\epsilon ^{})=\stackrel{~}{\varphi }\left(\frac{\epsilon ^{}L^{1/\nu _\epsilon }}{(A^{}\mathrm{ln}L+\xi )T_{c,L}}\right)=\stackrel{~}{\varphi }\left(\frac{\epsilon ^{}L^{1/\nu _\epsilon }}{C_LT_{c,L}}\right)$$
(61)
with a slightly modified scaling variable. This modified scaling variable in the energy picture is investigated for the two-dimensional Ising model in the next Section. Note however that the factor causing the suppression of the linear correction term of the magnetisation is now integrated into the scaling variable. This has the consequence that the reduction of the corrections to scaling in (61) are solely controlled by the relative size of the amplitudes $`B_L`$ and $`b_L`$. Therefore it has to be expected that the reduction of the corrections are weaker than those in the case of an algebraically diverging specific heat.

## IV Critical quantities

In this Section we extract critical quantities like critical exponents or transition temperatures of the infinite system from the molecular field data of the corresponding finite systems. Apart from the molecular field finite size scaling relations presented in the previous Section the notion of the effective critical exponent provides another way to obtain critical exponents from finite system quantities .

The effective critical exponent $`\beta _{eff,N}`$ of the order parameter is defined by expressing the temperature dependent spontaneous magnetisation $`m_N`$ in the form

$$m_N=At^{\beta _{eff,N}(T)},$$
(62)
where $`t=(T_cT)/T_c`$ and $`T_c`$ is again the critical temperature of the infinite system. The temperature dependent effective exponent is obtained by differentiating relation (62):

$$\beta _{eff,N}(T)=\frac{d\mathrm{ln}(m_N)}{d\mathrm{ln}t}.$$
(63)
Fig. 6 shows the effective exponent of the two-dimensional Ising model for the different clusters on the square lattice together with the effective critical exponent $`\beta _{eff,\mathrm{}}`$ of the infinite system . As the pseudo-critical temperatures of the finite systems are above the critical temperature $`T_c`$ of the infinite system the spontaneous magnetisation $`m_N(T)`$ is finite at $`T_c`$ and therefore the effective exponent has to vanish at $`T_c`$. It is also observed that the effective exponents of the largest clusters are closer to the corresponding curve of the infinite system than those of the small clusters. From the effective exponents estimates for the true critical exponent $`\beta `$ may be extracted by means of linear extrapolation. To this end an extrapolated exponent $`\beta _{ext,N}`$ is defined by

$$\beta _{ext,N}(T)=\beta _{eff,N}(T)+\frac{d\beta _{eff,N}(T)}{dT}(T_cT).$$
(64)
This extrapolated exponent is displayed in Fig. 7 for the square lattice again together with $`\beta _{ext,\mathrm{}}`$ of the infinite system. Note that the effective exponent depends on the critical temperature of the infinite lattice. For models with an unknown value of $`T_c`$ a reliable estimate has to be obtained from the finite system pseudo-critical temperatures. This can be achieved by using the finite size scaling relation (39).

In the next Subsections we investigate our molecular field data for the various systems using both finite size scaling relations and effective exponents.

### A Two-dimensional Ising models

Using the cluster molecular field approximations we studied two-dimensional nearest neighbour Ising models defined on the square and on the triangular lattices as well as an Ising model with both ferromagnetic nearest and next-nearest neighbour interactions defined on the square lattice. The extrapolated exponent from Eq. (64) allows us to obtain reliable values for the order parameter critical exponent $`\beta `$, as shown in Fig. 7 for the square lattice model with only nearest neighbour interactions. The best estimate for $`\beta `$ is obtained from the value of the temperature where the extrapolated exponents for the two largest systems ($`N=145`$ and $`N=181`$) still coincide, yielding $`\beta 0.126`$ which should be compared to the exact value $`\beta _{ex}=1/8`$. This overestimation of only $`0.5\%`$ is rather remarkable considering the smallness of the systems under consideration. Note however that the slight overshooting above the true value 1/8 is not an artifact of our approximation, it is due to the peculiar behaviour of the exact solution of the infinite system, see the inset of Fig. 7.

A similar behaviour is also observed for the other two-dimensional Ising models. In these cases the overshooting is somewhat larger, yielding a slightly increased value for $`\beta `$. Thus we obtain the value $`\beta 0.128`$ for the triangular lattice, still in very good agreement with the exact result.

Other critical exponents, as for example the critical exponent $`\gamma `$ of the susceptibility, might be derived in the same way. In the present study, however, we restrict ourselves to the small system sizes where the density of states can be computed exactly. For these systems the extrapolated exponents obtained from the susceptibility do not yet present a plateau, making an estimation of $`\gamma `$ rather tedious. Here, the molecular field finite size scaling theory of Section III is an interesting alternative for extracting critical quantities, as discussed in the following.

Let us start by probing the finite size scaling relations (43) and (61) of the spontaneous magnetisation in the temperature as well as in the energy picture. Note that relation (61) has to be used instead of (56) as the specific heat of the two-dimensional model presents a logarithmic divergence with $`\alpha =0`$. Plotting the scaled magnetisation as a function of the scaling variable and setting $`L=N^{1/d}`$, $`d`$ being the number of space dimensions, the best data collapse as derived from a least-square fit procedure or from a more sophisticated approach  yields the values of the critical exponents. The resulting data collapses are shown in Fig. 8 for the nearest neighbour Ising model defined on the square lattice where the smallest systems with 5 and 13 spins have been discarded. Hereby the critical exponents in the temperature picture take the values $`\beta =0.127`$ and $`\nu =1.01`$ (see Fig. 8a), while in the energy picture (see Fig. 8b) the best data collapse is achieved for $`\beta _\epsilon =0.124`$ and $`\nu _\epsilon =1.02`$. This convincingly shows that the finite size scaling relations derived in Section III are indeed fulfilled by our scaled data. Again, one has to stress the remarkable agreement between our estimates of the order parameter critical exponents and the exact value 1/8, based on data obtained for very small system sizes.

Fig. 9 displays the scaled susceptibilities as function of the scaled reduced temperature. The values of the critical exponents $`\gamma =1.77`$ (to be compared with the exact value $`\gamma _{ex}=1.75`$) and $`\nu =1.01`$ again result from the best collapse of the scaled data. The same can be done in the energy picture where we obtain $`\gamma _\epsilon =1.80`$ and $`\nu _\epsilon =1.02`$, in reasonable agreement with the exact values. This shows that the finite size scaling theory enables us to obtain reliable values for critical exponents also in cases where an analysis of effective exponents does not succeed.

At this stage we have to emphasise again that no non-universal quantities of the infinite system enter the molecular field finite size scaling theory. Especially, the exact location of the critical point is not needed. On the other hand, using the value of $`\nu `$ obtained from finite size scaling together with the pseudo-critical temperatures of the three largest systems listed in Table I, a remarkable good estimation of the critical temperature follows from Eq. (39): $`k_BT_c/J_I=2.266`$ in excellent agreement with the exact value $`k_BT_c/J_I=2/\mathrm{ln}(1+\sqrt{2})2.269`$.

Before closing this subsection, we note that almost identical estimates are obtained for critical quantities when studying the triangular model: $`\beta =0.127`$, $`\beta _\epsilon =0.125`$, $`\gamma =1.78`$, $`\gamma _\epsilon =1.80`$, $`\nu =1.03`$, $`\nu _\epsilon =1.03`$, and $`k_BT_c/J_I=3.631`$ (exact: $`3.634`$). The discrepancies are somewhat larger for the square lattice clusters where the next-nearest neighbour interactions have also been taken into account. We obtain $`\beta =0.120`$ and $`\gamma =1.86`$, for example, analysing systems with up to 81 spins.

### B Three-dimensional Ising model

When trying to extract critical quantities for the three-dimensional Ising model from the cluster molecular field approximation, we are faced with two difficulties: (1) only systems with a very small number of spins ($`57`$) can be handled exactly and (2) the various clusters have different geometries, as explained in Section II. Especially the second point makes a finite size scaling analysis much more problematic, as discussed below. But let us start again by looking at the effective exponent derived from the spontaneous magnetisation. Fig. 10 shows the extrapolated exponent obtained for the different system sizes. The estimated value of 0.316 for the critical exponent $`\beta `$ is slightly smaller than the expected value $`0.326`$. This is due to the finite-size effects which for our small systems set in already at rather low temperatures. Fig. 10 also displays an oddity due to the different cluster geometries: on the dashed line the extrapolated exponents of two different system sizes, $`N=27`$ and $`N=33`$, are almost indistinguishable. Looking at the raw data, one observes that, by coincidence, due to the difference in the geometries, the cluster molecular field approximation yields virtually identical spontaneous magnetisations for these two clusters.

The difficulties originating from the difference in geometry are enhanced when looking at the finite size scaling behaviour in three dimensions. In fact, the rather naive identification $`L=N^{1/d}`$, which works very well in two dimensions with identical geometries, is inadequate when different cluster geometries are involved. Sticking to this definition of $`L`$, we are only able to obtain simultaneously a good data collapse for all the quantities of interest (magnetisation and susceptibility as function of both reduced temperature and reduced energy) when the sizes $`N=27`$ and $`N=33`$ are omitted. This procedure then yields the following values for the critical exponents: $`\beta =0.32`$, $`\beta _\epsilon =0.41`$, $`\gamma =1.23`$, $`\gamma _\epsilon =1.45`$, $`\nu =0.63`$ and $`\nu _\epsilon =0.69`$. Whereas in the temperature picture the values of the critical exponents agree nicely with the values given in the literature ($`\beta =0.326`$, $`\gamma =1.24`$, $`\nu =0.63`$), some discrepancies show up in the energy picture (for example, one expects 0.367 for $`\beta _\epsilon `$), even if the overall agreement is still very satisfactory with regard to the smallness of our systems.

A possible way to circumvent some of the difficulties arising from the difference in cluster geometry would be to eliminate $`L`$ altogether in Eq. (43) resp. (56) by exploiting  the relation (39) resp. (57). However, this would introduce the necessity to determine precisely the location of the critical point beforehand.

An estimate of the critical exponent $`\alpha `$ of the infinite three-dimensional Ising system may be obtained by applying relation (50) for the maxima of the specific heat of the finite cluster systems. Using the system sizes 19 and 57 one obtains the value $`\alpha =0.132`$, again in reasonable agreement with the literature value $`\alpha =0.109`$.

In Section III we already discussed briefly the possible reduction of the correction to scaling terms when going from the temperature to the energy picture. Denoting the amplitude of the linear correction term in (55) by

$$K_L=L^{\alpha /\nu }\frac{1}{a}(B_Lb_L)$$
(65)
the ratio $`B_L/K_L`$ may serve as a measure of the suppression of the linear correction term in (55) compared to the corresponding term in (48). Note that the pseudo-critical temperature is not integrated into this ratio as $`T_{c,L}`$ appears in the scaling variable of the energy picture. For the magnetisation of the three-dimensional clusters this ratio is of the order 8. Similarly the ratio $`B_L/(B_Lb_L)`$ describes the reduction of the correction term for a system with a logarithmically diverging specific heat. Again the factors appearing in the scaling variable of (61) are excluded from the ratio. For the two-dimensional Ising system in the cluster molecular field approximation this ratio is of the order 6. In both cases the corrections to scaling are suppressed in the energy picture as expected from the considerations in Section III. It is also observed in our data that the suppression decreases slightly with increasing system size for the range of extensions $`L`$ investigated in our cluster molecular field approximation.

Finally, we again get a very good estimate of the critical temperature from Eq. (39). Inserting the value $`\nu =0.63`$ coming from the best data collapse of the scaled data as well as the pseudo-critical temperatures of the systems with 7, 19, and 57 spins, we find $`k_BT_c/J_I=4.5208`$ which should be compared to the literature value $`k_BT_c/J_I=4.5115`$.

### C Three-state Potts models

Let us close this Section by discussing three-state Potts models on square and triangular lattices. In these cases the cluster molecular field approximation displays a discontinuous transition for all finite clusters, as illustrated in Fig. 3 by the jump of the order parameters. At first sight this fact seems to rule out the possible determination of critical quantities from our small systems. However, we show in the following that our methods of analysing the data in principle still allow us to obtain reliable values for the critical quantities.

In the present situation the analysis of effective critical exponents has the great advantage that only data computed for temperatures lower than the critical temperature of the infinite system are involved. Therefore the discontinuous behaviour of our small systems observed at temperatures above the critical temperature does in this case not interfere with the determination of the critical exponents. Thus the extrapolated critical exponents computed for square lattice systems with up to 85 spins yield the value $`\beta =0.104`$ for the order parameter critical exponent, in good agreement with the exact value $`\beta _{ex}=1/9`$. On the other hand the molecular field finite size scaling analysis is deeply affected by the discontinuity of the phase transitions as here data obtained in close vicinity to the transition points play a prominent role. In fact, the main problem lies in the values of the pseudo-critical temperatures which enter the definition of the scaling variables. Using the transition temperatures given by the kinks in the cluster molecular field free energies, see inset of Fig. 3, does not lead to a collapse of the scaled data. A better result is obtained by fitting the spontaneous magnetisations of the finite systems to polynomials of second order and then inserting the temperatures where these polynomial functions hit the temperature axis into the scaling variables. This approach leads to an acceptable data collapse for the order parameter, as shown in Fig. 11, the critical exponents thereby taking the values $`\beta =0.11`$ and $`\nu =0.76`$. The value of $`\beta `$ again nicely agrees with the exact result but a slight discrepancy prevails for $`\nu `$ as the exact value is 5/6. With regard to the crudeness of our approach, this is still a very satisfactory result. However, we did not succeed in extracting the critical exponent $`\gamma `$ from the susceptibility data in a similar way. Finally, inserting the value $`\nu =0.76`$ into Eq. (39) together with the temperatures obtained from the polynomial fits, we end up with the estimate $`k_BT_c/J_P=0.998`$ for the critical temperature, very close to the exact value $`k_BT_c/J_P=1/\mathrm{ln}(\sqrt{3}+1)0.995`$.

## V Conclusion

The advantages and shortcomings of molecular field theory are well known. On the one hand this approximation provides a rather simple way of obtaining analytical results even in complicated systems. Usually, the results of the molecular field approximation are qualitatively correct. On the other hand, however, large quantitative discrepancies are found when comparing these predictions with exact results or with estimations obtained by numerical methods and/or more sophisticated analytical approaches. This is especially the case when considering systems with a continuous phase transition as the molecular field approximation always yields classical values for the critical exponents irrespective of, for example, the dimensionality of the system. Furthermore, critical temperatures are usually overestimated by a considerable amount.

In this work we have shown that the correct values of critical quantities can nevertheless be obtained from systematic improvements of the molecular field approximation. We thereby insert spin clusters of variable sizes into a homogeneously magnetised background and determine the molecular field acting exclusively on the border spins by imposing a condition of homogeneity in the central part of the cluster. The interactions inside the clusters are treated exactly by computing the density of states as function of the internal energy and of the interaction energy with the molecular field.

The data computed in this cluster molecular field approximation can then be analysed in various ways in order to extract the values of critical quantities. We have discussed in this work two rather complementary approaches, one approach involving effective critical exponents where only data located at temperatures below the critical temperature of the infinite system enter, the other being a molecular field finite size scaling theory (which in the temperature picture is related to that developed by Suzuki ) where the true critical exponents result from the best collapse of the rescaled data. In the latter approach the data located in close vicinity to the pseudo-critical temperatures are of main interest.

With these methods data resulting from the cluster molecular field treatment of various two- and three-dimensional Ising models as well as of two-dimensional three-state Potts models have been investigated. In all these cases we have shown that critical temperatures and critical exponents can be determined with a surprisingly high precision with regard to the smallness of our spin clusters. In two dimensions we have the nice situation that clusters of various sizes with the same shape can be handled by the available computational capacities. This is different in three dimensions where our computational resources do not allow us to handle spin clusters of various sizes having the same geometry. In fact, these restrictions are imposed by our choice to compute the density of states in a numerically exact manner. In principle larger system sizes could be investigated by Monte Carlo simulations, using one of the highly efficient methods for computing the density of states which have been developed recently . However, very good numerical data would be needed in order to obtain the strength of the molecular field as a function of the temperature from Eq. (18). Especially the accurate determination of the pseudo-critical temperature requires a density of states of extremely high precision not easily achievable in a computer simulation. We therefore expect that data obtained from Monte Carlo simulations may not permit a molecular field finite size scaling analysis. However, they should be very useful for the investigations of effective critical exponents, especially when dealing with quantities like the susceptibility or the specific heat for which the small systems considered in this work did not allow a determination of the values of the critical exponents.

Further progress along the lines developed here can be achieved by noting that the homogeneity condition in the central part of the cluster renders our cluster molecular field approximation superior to more traditional approaches. However, the spin expectation value is still not identical in the whole cluster. Here a scheme developed by Galam  could be very useful. In that cluster molecular field scheme a cluster consisting of the central spin and of its $`z`$ neighbours is constructed such that the homogeneity is extended to the whole cluster. It is, however, not immediately obvious how to generalise this scheme to larger clusters containing not only the central spin and its nearest neighbours. Nevertheless we expect that the construction of homogeneously magnetised clusters (be it by the Galam scheme or by a different one) would further improve the systematic cluster molecular field approximation and yield even better values for the critical quantities.

## Appendix

In this appendix we describe how the density of states is evaluated in a numerically exact way. For simplicity consider a square lattice which consists of $`L`$ rows of length $`L`$. In our case the density of states depends on two extensive quantities, the interaction energy $`E`$ of the spins inside the cluster and the interaction energy $`B`$ of border spins with spins $`\sigma _k=1`$ outside the cluster. In other physical contexts the density of states can depend on other extensive quantities, e.g. the interaction energy and the magnetisation. The density of states $`\mathrm{\Omega }(E,B)`$ of the cluster is worked out by successively performing the trace over the $`q^L`$ states of a row of spins. A similar method has already been discussed in the literature for a density of states depending only on one extensive parameter . Let $`k=1,\mathrm{},q^L`$ denote the configurations of the row of spins and $`l=1,\mathrm{},L`$ the rows of the cluster. Suppose $`l`$ rows of the cluster are already built up and the density of states of the incomplete cluster is given by $`\omega _l(k,E,B)`$. If the next row is added additional internal and external interactions of this row with the already built up cluster and the spins outside have to be included. Let $`_{k^{},k}`$ denote the interaction energy of the microstate $`k^{}`$ of the last row of the yet incomplete cluster with the microstate $`k`$ of the row to be added in the next step. Similarly $`_k`$ denotes the additional external energy. Of course $`_k`$ does not contain the contributions of those bonds of the new row which will become internal bonds when the next row is added. Then the density of states of the extended cluster consisting now of $`l+1`$ rows is given by

$`\omega _{l+1}(k,E,B)`$ $`=`$ $`{\displaystyle \underset{k^{},E^{},B^{}}{}}\omega _l(k^{},E^{},B^{})`$ (67)
$`\delta (E,E^{}+_{k^{},k})\delta (B,B^{}+_k).`$
The density of states of the complete cluster is then obtained as

$$\mathrm{\Omega }(E,B)=\underset{k}{}\omega _L(k,E,B).$$
(68)
The described method is easily generalised to densities of states depending on more than two parameters.