Recurrence relations for the three-dimensional Ising-like model in the external field

The method for calculation of the partition function of lattice model for the magnet in the external field near critical point (CP) is proposed. The recurrence relations and their explicit solution near the critical point are founded. It is shown that dependence on temperature of thermodynamic functions near CP, when the field value comes down to zero, is in good agreement with the previous results obtained using the collective variable method. The phase transition temperature (when h = 0) is calculated and the dependence on parameters of interaction potential is found.


Introduction
A significant progress in the theoretical description of phase transitions on the microscopical level has been achieved using the collective variables (CV) method [1].This approach allows one to take into account the collective behavior which plays a crucial role near the phase transition point.The effectiveness of the CV method was demonstrated by applying a one-component three-dimensional spin model [2] to the description of critical behavior.The critical temperature T c with explicit expressions for thermodynamic functions near T c have been obtained and their dependence on microscopic parameters of the system was found.However, the influence of external field on the above mentioned system is still unclear.The evaluation of expressions for free energy and other thermodynamic functions (heat capacity and order parameter in particular) near T c in the vanishing external field is of great theoretical and practical interest.The critical point for one-component spin system is defined at the temperature equal to T c and in zero external field h.The critical behavior of this system at h = 0 has been studied quite well.In particular, the description of such a system on the microscopical level is given in [2].We also know some of the characteristics for the one-component spin model at T = T c and with the external field tending to zero.Still unclear are the details of critical behavior of the spin system when T → T c and h → 0, but T = T c and h = 0.The solution of this problem is rather non-trivial, for instance, there is no exact solution even for the two-dimensional Ising model at non-zero external field.
In this study we develop an approximate method for the description of onecomponent 3D spin model near the phase transition point.This method is based on the microscopic theory of phase transitions developed in [2] which uses the collective variables set.The introduction of external field leads to a more general description of critical behavior, but the main ideas and calculation schemes of [1,2] are kept intact.
Similar problems of phase transitions occurring in binary alloys were considered by Gurskii [3,4].In particular, the first principles approach was developed to calculate the partition function and order-disorder phase diagrams for the above mentioned systems [5].
The behavior of uniaxial magnetic systems and some other objects studied by statistical physics can be described quite well by the 3D Ising model.The Hamiltonian of this model takes the following form where Φ(r l j ) is the interaction potential between l-th and j-th lattice sites, r l j = | r l − r j | is an interparticle distance, h = µH is the normalized external field.The variable σ l takes two values ±1.Let us consider a simple cubic lattice with the spacing c.In the following calculations the exponentially decaying potential is used.Here A is a constant, b is an effective interaction radius.The partition function of the system described via the Hamiltonian (1) can be written in the CV ρ k representation in the following form [2] Z = exp The summation in (3) is performed over the wave-vectors k within the first Brillouine zone Here N = N x • N y • N z is the total number of particles, β = 1/k B T is the inverse temperature, Φ(k) is the Fourier transform of the interaction potential, and J(ρ) is the transition Jacobian from the spin variables to the collective variables.If the external field h is present it takes the following form J h (ρ) = Sp e βh l σ l J(ρ − ρ) , where the transition operator J(ρ − ρ) is expressed as [1] J(ρ − ρ) = exp The integration in equation ( 5) is performed over Here the prime means that k > 0.
The operators ρ k are The summation here is performed over the periodic volume with periodic boundary conditions .The explicit expression for transition Jacobian can be found as the result of summation Sp over eigenvalues σ l = ±1 and performing the integration over variables ω k in equation (5).In this way we obtain the known result where the site collective variables are introduced as The volume element in the CV space ρ k and site variables ρ l are related via the following expression where the transition Jacobian j is The change of variables from ω k to ω l leads to where The evaluation of partition function (3), free energy, and other thermodynamic functions near the critical point demands some approximations.It is connected with the fact that we can decompose expression (3) into two parts.The first, energetic part, is diagonal in terms of the variables ρ k , while an entropic part connected with transition Jacobian (7) is diagonal in the space of the site CV ρ l .At present, we do not have adequate mathematical equipment for an exact calculation of expression (3).Usually, the approximate methods of calculation have been used.One of such approximations in studying the critical behavior of the statistical systems consists in using the Gaussian distribution for the order parameter.In such a way, for the scalar ϕ 4 theory in d = 3 [6], the values of the critical exponents have been obtained by means of resummation of the series of Gaussian perturbation theory.At present these values are considered to be the most reliable.They are used as the basis for investigation of other objects, for example, the weakly quenched disorder Ising model [7].
Although the Gaussian distribution of fluctuation prooved to be beneficial in calculating the critical exponents and other universal quantities, it does not permit us to obtain non-universal parameters of the phase transition, for example, its critical temperature T c .The calculation of the non-universal quantities is connected with the use of non-Gaussian distributions in calculating the free energy [8,9], or with the use of some non-perturbative approach describing the critical properties of three-dimensional systems [10] accounting for the non-Gaussian fluctuations of the order parameter.The use of non-Gaussian distribution of fluctuations is especially important near the critical point of the second order phase transition for three-dimensional systems.The peculiarity of the method which uses non-Gaussian distribution of fluctuations is the so-called intermediate integration [1] that allows us to obtain an analytical expression for free energy near the critical point.In the present paper we generalize the method of the works [1,2] to the case of non-zero external field.In this case the system always has a non-zero order parameter.It is interesting to describe the behavior of thermodynamic functions in the case of non-zero external field.

Representation of the partition function
Let us write the functional representation for the partition function of the model (1), which will be useful in our subsequent calculation of thermodynamic functions near the critical point.We start with expression (3), for which (7) holds.Fourier transform of the interaction potential (Φ(k)), appearing in (3), has, according to (2), the following representation [1] where We are interested in the long-wave limit of Φ(k) because the critical behavior is determined by the long-range correlations.Therefore, we use the approximation which corresponds to the parabolic approximation of (1.1) for small wave vectors with subsequent averaging of the Fourier transform of (1.1) over the wave vectors near the boundary of the Brillouine zone (4).In the investigation of universal parameters of the model (1) such as its critical exponents, the value Φ 0 is unessential and may be put to zero.But the value Φ 0 is essential in calculating non-universal quantities, for example, the critical temperature [11].The definition of "small" values of the wave vector is ambiguous and depends on the form of the interaction potential.For the exponentially decreased potential (2) the region B 0 , where the parabolic approximation (1.3) holds, has the form where N 0x N 0y N 0z = N 0 , N 0 = N • s −d 0 , and s 0 1.The parameter s 0 determines the period of some effective block lattice c 0 = c • s 0 .
The parameter s 0 is determined differently for different interaction potentials (exchange interaction, nearest-neighbor interaction, etc.) provided that in the region k ∈ B \ B 0 , the dependence Φ(k) on the wave vector should be the weakest.In any case, we can consider that we investigate the critical behavior of the system with interaction potential (1.3) representing the long-range type (in particular, exponentially decreased) of the interparticle interaction.
In virtue of (1.3), the partition function ( 3) is presented in the form Similarly to the previous results ( 9), (11), we have where Here l belongs to the volume of periodicity (V = N 0 c 3 0 ) with periodic boundary conditions.Expression (1.5) for the partition function allows us to perform integration over variables ρ k for k ∈ B. To this end, we have a transit to the cite CV ρ l defined by (8).Let us introduce variables ω k for k ∈ B: Integrating (1.5) over the variables ρ l gives where and ω l is defined by (12).Now we use the cumulant series for ch(. . .).In virtue of [1] we have ch where h = βh.The quantities D n (ω l ) are given by where cumulants M n (h ) have the form (1.13) Making use of ( 12) and (1.13), we obtain an explicit form of the expression (1.10) for the Z(ω) The summation over wave vectors in (1.14) is performed for k ∈ B. But the variables ω k , in virtue of (1.8), are different from zero only for k ∈ B 0 .Therefore, the sums in (1.14) must be calculated only for k ∈ B 0 .Then, we replace Kronecker symbol In accordance with the above mentioned simplifications, one may recast the expression (1.9) for the partition function in the form where We have the following expression for the I l (η l ) Here The n 0 determines the number of terms of the exponent in (1.14) and defines the type of "model" -the order of approximation used for a concrete calculation.The case n 0 = 2 corresponds to the Gaussian approximation.In this case, since M 2 (h ) is positive for all values of the h, the integrals over ω k -variables are finite for all values of the field.When n 0 decreases (n 0 = 4, 6, 8, . . .), the type of the model ρ n 0 complicates.For an exact calculation we have put n 0 → ∞.However, in the real calculation we use finite n 0 .It is important that for small h all M 2n (h) have such signs which ensure finiteness of the ω k -integrals in (1.16).When the field decreases, the cumulants M 2n (h) with n 2 change their signs (see Appendix 1).This indicates the nonstability of the model ρ 2n .For example, the model ρ 4 is stable for h ∈ (−h c , h c ), where h c = 0.658.The value of magnetic field which corresponds to the h c , is given by H c ≈ h T c • 10 4 oersteds.Here T c is dimensionless number which is equal to the value of absolute temperature of the phase transition 10 2 ÷ 10 3 .Then H c ≈ h • (10 6 ÷ 10 7 ) oersteds.Comparing the value of H c with the field of saturation magnetization for iron, H Fe = 1.99 • 10 4 oersteds, we obtain h Fe ≈ 0.01.Therefore, the value h c ≈ 0.658 corresponds to very strong magnetic field and the model ρ 4 can be applied to the description of the critical properties for many real objects.
For convenience of presentation we perform in (1.16) the change of variables where and the following denotations have been used For compact writing (1.17) we define where the coefficients a n have to be calculated by means of the formulas [12] e In calculating the quantities L n , some approximation will be used since the expressions for ϕ(ν l ) and U (ν l ) are actually infinite series.These series may be considered formally as expansions in some small parameter s −d 0 (s 0 2).It is known from the previous publications [1,2,11] that every term in the function U (ν l ) leads to a concrete approximation.So, the choice g = f = • • • = 0 gives us the Gaussian approximation with classical critical exponents.The condition g = 0 and f = • • • = 0 corresponds to the model "ϕ 4 ", g = 0 and f = 0 leads to the model "ϕ 6 " etc.
We may expect that the similar situation takes place for the function ϕ(ν l ).In virtue of (1.19), its simplest approximation is as follows: (1.23) Hereinafter this approximation will be called first odd cumulant approximation since the coefficient a is proportional to the cumulant M 1 (h ).The expression will be called second odd cumulant approximation etc.Now we will perform an approximate calculation of the coefficients a n from (1.21).It should be noted that all a n are real numbers.Hereinafter for the quantity n 0 from (1.20) we put n 0 = 4.It means that our subsequent results will correspond to the model "ϕ 4 ".The generalization to the case of "ϕ 6 "-model etc. one may perform by means of the results [11,13].Concerning the function ϕ(ν l ) for simplicity we will use the first odd cumulant approximation, for which (1.23) holds.In this case we obtain (1.25) The method for calculating the coefficients a n is given in [13].The prime in the denotation of the coefficients a n indicates the use of the first odd cumulant approximation 1 .We note that coefficients a 1 and a 3 are proportional to the external field h . 1 The coefficients a n in the second odd cumulant approximation are given in Appendix 2.
The partition function of the model is written in the form where while for the coefficient ã2 we have ã2 = a 2 + βΦ 0 . (1.28) Expression (1.26) is our starting point in step-by-step calculation of the free energy for the Ising model with potential (1.3) near the critical point.In contrast to the works [1,2,[11][12][13] here the external field appears explicitly and leads to the appearance of odd powers of the variables η k in the exponent.The coefficients (a 1 , a 3 , . . . ) near odd powers of η k tend to zero in the limit h → 0. In the subsequent calculation, especially in obtaining recurrence relations, the approximation used for initial coefficients a n is unessential.However, the fact of the appearance of even and odd powers of η k in the exponent (1.26) is important.Initial values of the coefficients a n and their relation with the field h become significant only in calculating the observable quantities.

Method of calculating the partition function
We perform a step-by-step calculation of partition function (1.26), beginning with integration over the variables η k with wave vectors k near their maximum value B 0 and ending with integration over η k with | k| → 0. We use the method which has been proposed in [1].
Let us introduce Brillouine zone where Here N 1 is number of sites of the block lattice, and N 1 = N 0 s −d .In the same way we define the block lattices with periods c 2 = c 0 s 2 , c 3 = c 0 s 3 , . . ., c n = c 0 s n , which contain N 2 = N 0 s −2d , N 3 = N 0 s −3d ,. . ., N n = N 0 s −3n sites, respectively.The parameter s controls the growth of the block structures with periodicity volumes where Accordingly to [2] we select in (1.26) the variables η k with k ∈ B 0 \ B 1 and average the quantity Φ(k) over such values of wave vectors.Then Here Φ(B 0 , B 1 ) denotes the mean value of the potential Φ(k) in the region k ∈ B 0 \B 1 .For the mean-arithmetical averaging we have After such transformations the partition function takes the form where The set of N 1 variables ν k (ν 0 , ν c k , ν s k ) determines intermediate integration which allows one to perform the integration in (2.4) over N 0 variables η k .To this end, we transit to the site variables with transition Jacobian equal to j −1 0 , which cancels the corresponding factor in (2.4), and introduce the quantity where The result of integration of (2.4) can be written as where J(ν l ) is given by Hereinafter we shall use approximate expressions (1.25) for the coefficients a n .
For the subsequent calculations it is convenient to perform in (2.9) the change of variables [11] η l = 24 a 4 1/4 x.
Then we obtain for the J(ν l ) where

11)
with (2.12) Now reexpress the quantity T (ν l ) into the form The coefficients S n may be calculated in much the same way as the coefficients a n from (1.20).We have a system of equations The left side of the equation (2.14) has the form where The derivatives of the right side of (2.14) are , where e −S 0 = K 0 (h 2 , h 3 ). (2.18) The coefficients S n from (2.13) may be obtained from the relations

.19)
Approximate explicit expressions for the coefficients S n may be obtained assuming small values for the quantities h 2 , h 3 , which are arguments of the functions K n (h 2 , h 3 ) from (2.16).That h 2 is small in the critical region has been proven in [2,13], while h 3 is proportional to the field and tends to zero near the critical point.Taking this into account we obtain an approximate expression for (2.18) where γ = Γ 3 4 2 /π √ 2 ≈ 0.337989.Here and henceforth we do not consider the terms which are proportional to h n 2 with n 2, and to h n 3 with l 3, since they are inessential near the critical point.
In the approximation when h 2 and h 3 are small we have: where To simplify the calculation we will use the linear approximation of (2.21) neglecting the terms which are proportional to h 2 3 and h 2 h 3 .After integration (2.8) over the variables η k the partition function takes the form where use is made of (2.7) and the definitions where The final step of the calculation consists in integrating over the N 1 variables ν k .To perform this we transit to the site variables in (2.22) and integrate.The result is (2.26) where For the coefficients P n we have the following expressions where The change of variables x.
in (2.28) allows one to obtain the expressions for L m (ρ) where where We put (2.29) into the form where 2 = (a 1 ) 2 + e a (1) (1) 1 a (1) 2 − (a 2 Similarly to our calculation of the approximate expressions for the coefficients a n (see (1.21)-(1.25)),we can write approximate expressions for a n .The following correspondence should be noted: Moreover, the difference consists in the change of sign near the cross term x • ρ m in the expression (2.30) and near the corresponding term η l ν l in (1.17).Approximate expressions for R n are: Inserting the obtained expressions into (2.34),we find approximate formulae for the coefficients a (1) n : where e −a (1) 0 . (2.38) The coefficients d 1 (k) and ã(1) 1 have the form (2.39) Comparing (2.37) with (1.26) we can see that the functional form of the partition function did not change.The number of integration variables decreases (from N 0 to N 1 = N 0 s −d ) and the coefficients d(k) and a n change their values.

Recurrence relations
Now we shall write an explicit form of the recurrence relations (2.36) expressing the coefficients a (1) n in the term of their initial values a n from (1.21) or from approximate results (1.25).Making use of (2.32) and (2.28) we obtain from (2.36) where = f 00 (a 4 ) Here we denote: Taking into account (2.12), we obtain from (3.2) 2 = f 00 a Let us perform in (2.37) the change of variables As a result, the k-dependent part of the Fourier transform of the interaction potential in (2.37) will change from k 2 to (sk) 2 .Taking into account that k ∈ B 1 , we can see that s k belong to the set B 0 , as it was for the initial expression (1.26).Therefore now we may compare other coefficients before and after the integration.
The change of variables (3.5) leads (2.37) to the form where the coefficients are given by 1 + s d/2 a 1 ], 2 − βΦ(0)], ; Taking into account (3.4), we obtain an explicit form of the recurrence relations where The quantity q is given by Now we find the fixed point of the RR (3.7) from the conditions The last equation of (3.7) gives where h * 2 = √ 6(r * + q)(u * ) −1/2 .Since f 01 < 1, there always exists such s = s * , for which h * 2 vanishes.As s is near s * , the quantity h * 2 is small 2 .Calculating the partition function under s = s * = 3.3783, we achieve an essential simplification since we need only a few terms in the series expansion on the powers of h 2 .By virtue of third equation (3.7) we find for (3.8)The quantity h * 3 = h 30 v * /(u * )3/2 is proportional to v * .Hence, at the fixed point with s = s * we have Using the second equation (3.7) we obtain where the quantity q is defined in (2.6).First equation (3.7) leads to the condition w * = 0.
Therefore, we have the following coordinates of the fixed point of RR (3.7): Here q = βΦ(0) • q and q = π 2 (b/c) 2 s −2 0 (1 + s −2 ).Now we return to the expression (2.37).In much the same way as (1.26) it may be integrated over the variables ρ k with indices from the domain where B 2 is given by Here The calculation is analogous to the equations (2.4)-(2.39).Now the coefficients a n have to be replaced by a n , for instance: After performing (n + 1)-th step of integration (1.26) we obtain The sums in (3.14) are taken over k ∈ B n+1 , where n+1) .The partial partition functions Q n are of the form where Then we obtain the following recurrence relations where the quantities f 00 and f 01 are given in (3.3), (3.20) In the case of absence of the external field (h = 0) RR (3.19) contain only second and fourth equations.This case has been investigated by us in detail in [2,13].

Solution to the recurrence relation near critical point
Comparing RR (3.19) with (3.7), we note that they both have the same fixed point (3.12).Next we put s = s * , where s * is given by (3.9).It is convenient to write (3.19) in a matrix form The matrix R has the following entities, R ij , The matrix R possesses four different real eigenvalues: We note that eigenvalues E 1 and E 3 are connected with the presence of the external field.Both these eigenvalues are greater than unity In the case h = 0 recurrence relations (3.19) simplify so that E 1 and E 3 vanish and a saddle fixed point appears.We obtain the following value for the critical exponent of the correlation length ν = 0.609.
It is slightly lower than the result of numerical estimate of this exponent for the Ising model, ν c = 0.630 (see, e.g., [6,15]).The difference in the values of ν and ν c is considered to be rather our limitation of the ρ 4 -model than the approximation during its calculation.Indeed, the Ising model may be really described by the model ρ 2m with m → ∞ [2].To achieve a real correspondence between ν and ν c we have to use at least the model ρ 6 .
Eigenvectors W ik of the matrix R are determined by the system of equations Using (4.2) gives where It is known from the matrix theory [16], that the nonsymmetric matrix R with different eigenvalues can be expressed in the form where the raws of the matrix T are eigenvectors of the matrix R, T −1 is matrix inverse to T , so that where I is a unit matrix, and Λ is diagonal matrix with eigenvalues of the R on the main diagonal.
The n-th power of the matrix R is given by where the matrix Λ n is also diagonal with the n-th powers of the eigenvalues of R on the main diagonal.Making use of (4.7)-(4.9)allows us to rewrite (4.1) in the form where the raw x n has the form To determine x 0 we have put n = 0 in (4.11) .Taking into account (4.9) and (4.10) yields 12) The latter expression allows us to obtain an expression for coefficients w n , r n , v n and u n (as well as a (n+1) n and d n ) from (3.20) in terms of the initial values a 1 , d(0), a 3 and a 4 given by (1.26).To this end we have obtained an inverse matrix T −1 .It is built with eigenvectors of the matrix transposed to the matrix R which are determined by the equations Therefore we obtain where the quantities T ij are the same as in (4.5).The normalization conditions j Now, using (4.12) we find an explicit form of the vector x n .Noting that and performing other calculations from (4.12) allows us to find where the denotations have been used.The quantity V 2 is given by (4.16).It is easy to verify that in the case n = 0 equations (4.17) hold identically.When h = 0, the system of equations (4.17 Therefore we obtain an equation for the temperature of phase transition: Figure 1 shows the dependence of inverse temperature β c Φ(0) on the parameter b/c.For convenience a scaling factor has been used since Φ(0) = 2dJ, where J corresponds to the constant of nearest-neighbor interaction.For the Ising model with nearest-neighbor interaction it has been obtained [17,18] that  It is easy to see that the value (4.21) is approached at b/c ≈ 0.3.Figure 3 gives the dependence β c Φ(0)/6 versus s 0 at b/c = 0.3 and Φ = 0.092.
Our calculation of the values T c was not intended to obtain an "exact" temperature of phase transition.Rather, it was necessary for determining the coordinates of critical point (T = T c , h = 0), near which we wish to investigate the thermodynamic functions of the system with the appearance of the external field.Now we introduce the denotations Then the solutions to RR (4.17) take the form

Conclusions
The paper gives the expression (3.14) for the partition function of the onecomponent magnet in the external field near the critical point.Explicit expressions for the partial partition functions Q n (3.15)-(3.17)are obtained.This allows us to calculate the free energy and other thermodynamic functions of the system near critical point by means of the methods of [2,8,9] .They will depend on the temperature and field and the form of these dependences as well as the corresponding critical exponents will be determined by the recurrence relations (3.19) and their solutions (4.22) near the fixed point (3.12).

Appendix 1
Let us consider the dependence of the cumulants M n (βh) from (1.12) on the value h = βh.We rewrite the expression (1.15) for the partition function in the form: The parameter s 0 1 determines the form of potential (1.3), which is used in calculations.To any value of s 0 there is a corresponding model system with its own values of parameters.The curves of dependence of M n (h) on the field h are shown in the figures 4 to 10 at the value s 0 = 2.For brevity the primes near M n and h are omitted.
The convergence of the integrals in (A1.1) is determined by the signs of the even cumulants M 2l (h ).The cumulant M 2 (h ) is positive for all values of the h .The M 4 (h ) is positive for h ∈ (−0.658; 0.658) and negative for |h| > 0.658.The quantity M 6 is positive everywere except the values |h| ∈ (0.421; 1.575).

Figure 1 .
Figure 1.The dependence of the β c Φ(0)/6 on the ratio of the range of interaction b to the lattice constant c under s = 2 and Φ = 0.092.
) reduces to two equations (w n = v n = 0).Since E 4 < 1, the E n 4 rapidly decreases as n grows.The values of r n and u n will tend to their fixed values under the condition r 0 − r * − T 24 (u 0 − u * ) = 0, which determines the temperature T c of the phase transition.An explicit form of this equation is given by ã2 .21) This value can be recovered by means of direct calculations with some set of parameters Φ, b/c, and s 0 .Figure2shows the dependence of the β c Φ(0)/6 on the parameter Φ.

Figure 3 .
Figure 3.The dependence of inverse temperature of phase transition on the parameter s 0 .

Figure 4 .
Figure 4. Dependence of the coefficient M 0 on the field h.

Figure 5 .Figure 6 .
Figure 5. Dependence of the first odd cumulant M 1 on the field h.

Figure 7 .Figure 8 .
Figure 6.Dependence of the second cumulant M 2 on the field h.Figure 7. The cumulant M 3 versus the field h.

Figure 9 .Figure 10 .
Figure 9. Dependence of the M 5 on the field h.
.36)Here the quantities G and G h are defined in (2.32) while P 2 is given in (2.28).