- Email: [email protected]
Cluster approach to dilute magnetism
Phycisa A 152 (1988) 431-450 North-Holland, Amsterdam
CLUSTER
APPROACH
Paulo Departamento
R.C.
TO DILUTE
HOLVORCEM
and Roberto
MAGNETISM OS6RIO
de Fisica, Universidade de BrasrXa, 70910 Brasilia, DF, Brazil
Received 5 February manuscript received
Revised
1988 22 April
1988
A cluster algebra is developed for the definition of independent correlation functions in the cluster-variation method (CVM) for the spin-l Ising model. A scheme is then introduced for the study of site-dilute spin-i Ising models by means of the CVM. The procedure regards the site-dilute spin-i model as the spin-l model with additional constraints due to dilution. The Desjardins-Steinsvoll algorithm is used for the transformation of the CVM equations into a set of differential equations for the independent correlation functions with the inverse temperature as parameter. The evolution of the correlation functions with temperature and the behavior of response functions such as the specific heat and the susceptibility are then obtained for any degree of dilution. As an introduction to this scheme, its detailed application is presented here for the simple case of the pair approximation.
1. Introduction In spite of the success of series-expansion and renormalization-group methods for the treatment of critical behavior in dilute magnets’), closed form approximations for these systems are still an important subject of investigation due to their potential to lead to correct phase diagrams and to the general (non-critical) behavior of thermodynamic functions with little computational effort. We are concerned in this work with the problem of quenched site dilution in the Ising model, defined in a given crystal lattice by the following Hamiltonian: H = -J
2
UjUo;‘lirlj
-
(ii)
where
J is the coupling
external
rli =
field,
1, 0,
cj = tl
C uiirli
(1.1)
>
I
constant between nearest-neighbor (nn) spins, and vi is a random variable with values
with probability with probability
0378-4371/88/$03.50 0 (North-Holland Physics
B
p , 1- p .
Elsevier Science Publishers Publishing Division)
B is an
(1.2) B.V.
P.R.C. HOLVORCEM
432
The external formally tibility. Among
taken
field will be taken into
the closed
consideration
as zero
throughout
in the discussion
form approximations
methods were the first to lead pioneering works used a two-site
AND R. OS6RIO
this paper,
but it will be
of the zero-field
for site-dilute
Ising systems,
suscepcluster
to a non-zero percolation threshold. The cluster approach in the spirit of the Bethe-
Peierls2.‘) or the constant coupling4) approximations. Several extensions of these methods have been later carried ou?). More recently, effective field theories. which use approximations based on Callen’s identity, have been extensively applied to this class of problems”~‘). The results of such calculations, however, have the limitation of depending only on the coordination number of the lattice, not taking into account its dimension or other geometric properties adequately. Kikuchi’s cluster-variation method (CVM), on the other hand, is a hierarchy of closed form approximations for the configurational entropy of classical lattice systems that take into account the geometry of the problem in a more complete way “‘-I*). When applied to pure (non-dilute) systems, the CVM results are seen to converge to the exact solutions when the size of the basic cluster is increased. In its simplest forms the CVM reduces to the mean field (one-site cluster) and Bethe-P eier . 1s ( two-site cluster) approximations. Only very recently approximations beyond the two-site cluster have been carried out in the CVM for the study of dilute systems. Osorio et al.“) have discussed the bee antiferromagnet with nn and next-nearest-neighbor interactions within the tetrahedron approximation. Moran-Lopez and Sanchez”) have given a phase diagram for the fee dilute ferromagnet also within the tetrahedron approximation and discussed the results in the context of a model of ‘He-“He mixtures with a frozen chemical short-range order parameter. The starting point of the CVM is the definition of the basic cluster, whose size limits the description of correlations between sites and determines the accuracy of the approximation. We define {a} as the set of all independent concentrations of configurations of the basic cluster. For instance, in the case of the pair approximation, the different basic cluster configurations are given by {(M, N)}, where M, N = +1 for the spin-$ Ising model. In this case, { cr} is the set of independent pair concentrations x2(M, N). Notice that the site concentrations x,(M) are not independent variables, since x,(M) = C,,,x,(M, N). The normalization C,,, x2(M, N) = 1 decreases the number of independent pair concentrations to three. The symmetry requirement interactions, decreases this x,(1, -1) =&-I, l), in the case of ferromagnetic number to two. In the canonical ensemble, the equilibrium state of the system at a given temperature T is obtained through the minimization of the free energy
CLUSTER APPROACH
TO DILUTE MAGNETISM
433
where
fTT; {a>) = E({fll) - W{(+)) > and where E( { c}) is the total energy due to the magnetic lattice and S({a}) is the “configurational entropy”,
(1.4) interactions
in the
where g( { c}) is the “statistical weight”, which is the CVM approximation for the number of equivalent arrangements of up and down spins of the whole lattice consistent with the given set of concentrations {c}. In the case of pure systems, the system of equations that follows from the minimization of the CVM free energy is usually solved by either the NewtonRaphson (NR) method15) or by the natural iteration (NI) method”). The NI uses as minimization variables the concentrations of all configurations of the basic cluster. Constraints that appear among the variables due to normalization or “translational symmetry” requirements16) are introduced by means of Lagrange multipliers. On the other hand, in the NR procedure the constraints are used to explicitly reduce the minimization variables to a set of independent “correlation functions”, which are defined by means of a convenient “cluster algebra”15). A recent interesting alternative to both the NI and the NR method is the Desjardins-Steinsvoll (DS) algorithm”), which consists of converting the problem into a set of differential equations for the correlation functions with the inverse temperature as variable. The known uncorrelated results for infinite temperature are used as initial values for an iterative procedure that leads to the thermodynamic functions at all temperatures. In this work we propose the application of the DS algorithm to the site-dilute Ising model. The independent correlation functions to be used in the DS algorithm are defined, in the next section, by means of a generalization for spin 1 of the spin-i cluster algebra. The spin-l problem is then reduced to the site-dilute spin-i model through additional constraints that describe the random distribution of non-magnetic sites. The resulting set of equations is then solved by the usual Runge-Kutta method (instead of the procedure for solving these differential equations proposed by DS). The critical temperature and thermodynamic functions such as the susceptibility and the specific heat are directly calculated as functions of the degree of dilution. This scheme is illustrated for the case of the two-site cluster approximation in the fee lattice.
434
P.R.C.
HOLVORCEM
AND
R. OSORIO
2. Cluster algebra for spin 1 In this section
the Sanchez-de
ized for the spin-l
Ising model,
Fontaine
spin-4 cluster
of which the site-dilute
algebra15)
spin-i
is general-
model is a special
case. The spin-l model, usually known as the Blume-Emery-Griffiths model and in a particular case as the Blume-Cape1 model, is of course also of interest by itself and has been applied to certain magnetic systems”), 3He-4He mixtures’“), ternary fluids2”) and other systems that display tricritical behavior. The present cluster algebra permits a straightforward definition of the correlation functions, which are the independent variables for the minimization of the CVM free energy in a given approximation of a spin-l lattice system. It also gives directly the relation between the concentration of a given cluster configuration and the correlation functions. If we denote by Sj = 0, +l the possible states of site i and by T(M, S,), M = 0, +-1, the projection operator, which takes values 1, when S; = M, and 0, when Si # M, then it is simple to verify that
T(M, Si) = (1 - M2) + $MS; + (qM* - 1)s; Let us label the clusters involved in the let n(r) be the number of sites of the specified by a set of iz(~) spin values configurations by an index 1= 1,2, . configuration
of the rth cluster
(2.1)
calculation by an index Y = 1,2, . . . and rth cluster. Each of its configurations is {M, , M,, . , M,,,,}. If we label these s 3”“‘, then the concentration of the Ith
is
Xi’-)= (T(M, >4 MM,, S,>. . . T(M,,,, , S,,,,)) ,
(2.2)
with the ( . . . ) average taken over the whole lattice. In order to find the independent variables for the CVM, we substitute the projection operator (2.1) into eq. (2.2). Thus we obtain each concentration as a linear combination of correlation Y, =
functions
m, )Yq-
y, of the form
. . . &,)P”‘)
,
(2.3)
where t is a label for the set of values {a,, a2, . . . , CQ,,}, and (Y~= 0,l or 2. The coefficients of this linear combination depend on 1 = {M, , M,, . . . , M,,,,} . We have at most 3”“’ - 1 independent correlation functions (not counting y, = 1, where we define t = 0 as the label for the set of a-values (0, 0, . . . , O}). This number corresponds to the maximum number of independent configurations of the basic cluster r = T,,, considering the normalization C, x,(T~) = 1. However, this number is usually greatly reduced due to symmetry relations
CLUSTER APPROACH
435
TO DILUTE MAGNETISM
(e.g., t = (0, 1) = t’ = {l,O}). If we label the independent correlation functions that appear in the expressions for x1(r) by an index 4 = 1,2, . . . , q,, we find an expression of the form
(2.4) To give the general expression for the z+(r, q), let us define the functions A,(M)
= 1 - M2,
A,(M)
= $M,
Then it can be shown, from eqs. (2.1)-(2.5),
A,(M)
= ;M2 - 1.
(2.5)
that
(2.6) where the sum is over the indices t of the correlation functions yI that are equivalent by symmetry operations to the same independent yq. At this point, some of the y4 may have their values constrained by parameters of the problem that are fixed beforehand. This occurs, for instance, for a quenched random distribution of Si = 0, to be discussed in the next section in the context of a site-dilute Ising model.
3. The site-dilute
Ising model in the pair approximation
As a simple application of the above algebra, we consider here a case where the basic cluster is an nn pair. We study the site-dilute Ising model with quenched, randomly distributed vacancies and ferromagnetic interactions. The clusters involved are the site (r = 1) and the nn pair (r = 2). This level of the CVM hierarchy is equivalent to the Bethe-Peierls approximation and is expected to be exact in the central region of a Cayley tree. The configurational entropy (per site in units of k,) can be written as
s =
3/l
l$
P,(l)ax,(l)l+
3/2 ,$
PIwam)l
3
(3.1)
with Z(X) = x In (x). Here /3,(r ) is the degeneracy of the Zth configuration of the rth cluster, i.e., the number of configurations that are symmetry-equivalent to each independent one labeled as I, and the 7,. are geometrical factors introduced by Barker*l). For a coordination number z, we have x = -z/2 and
436
P.R.C.
HOLVORCEM
AND
R. OS6RIO
TABLE I
Configurations
of the site (Y = 1) and the pair (r = 2) clusters and their degeneracies P,(r),
M
I
P,(l)
(M> N) (1,‘) (1, -1)
0 _
1 1 1 _
_
_
_
_
1
1 2 3 4 5
-1
6
P,(2) :. 2
(1.0) (-1, -1) (-1.0) (07 0)
: 1
y1 = z - 1. The degeneracies of the three configurations of the site and of the six configurations of the pair are listed in table I. The independent spin-l correlations that occur in the expansion of eq. (2.4) for r=1,2 are
and the coefficients I+(Y, s) are listed in table II. (Note that in tables I and II we use M, N instead of M,, M2.)Since in the dilute model the vacancies are quenched with concentration 1 - p, we have the additional constraints
y2 =
(s:) =p,
so that there
5, =y,, From
Y,
are only three
independent
52= Y,l
eqs. (2.4),
variables
(3.3) left:
53=y4.
(3.2)-(3.4)
x,(l) = i(P + 5,)> Coefficients
= =P2>
and table
(3.4) II, one gets
x,(l) = t(P - t,>,
+(1)=1-p,
TABLE II v,(I, y) of the expansion (2.4) of the cluster configurations the correlation functions y, of eq. (3.2).
(3.5)
x!(r) in terms
4
v,( 1) 4)
GL 4)
0
1 - M’ Ml2 3M’i2 - 1 0 0 0
(1 - M’)(l - N’) (M/2)(1 - NZ) + (N/2)(1 - MZ) (3M’/2 - l)(lN’) + (3N’i2l)(l ~ M’) MN14 (M/2)(3N’/2 ~ 1) + (N/2)(3M’/2 - 1) (3ML/2 - 1)(3N’i2 - 1)
1 2 3 4 5
of
CLUSTER APPROACH
TO DILUTE MAGNETISM
il it.
437
(3.6)
5, P21
Let us call A the above matrix of coefficients Al,q+l = ~~(2, q) (with l= 1 ,*.., 6; q = 0, 1, . . . , 5), which is used in the differential equation method of solution of the next section.
4. Differential equation method of solution Our solution of the model follows the general method proposed by Desjardins and Steinsvollr7), which converts the problem of minimization of the free energy into an initial-value problem of ordinary differential equations. For a fixed spin concentration p, the correlation functions si depend on the temperature. The free energy (per site in units of the coupling constant J) can be written in terms of the l-variables as
(4.1)
where CY,= - BIJ, CY~= -z/2, CQ= 0,, 7 = k,TIJ, and s depends on ti through eqs. (3.1), (3.5) and (3.6). The condition of minimum free energy, df/dg,= 0, is written as (Yi/T = as/&&
(4.2)
)
which is differentiated
with respect to the parameter
t = 1 /r to yield
(4.3)
where ii = d$ldt
.
(4.4)
P.R.C.
438
This differential ki=C
equation
Ciiai,
HOLVORCEM
AND
R. OS6RIO
can be put in normal
form,
i-1,2,3,
(4.5)
j=l
where
Cij = (D-l),.
For our model
the D, are
(4.6) with
Knij = Pn(2)(‘;,An,+ ‘iAn
+ ‘j,An,)(‘j,An, + ‘jzA.4 + ‘j4.5) >
(4.7)
where S, is the Kronecker delta. The accessible region in the space of coordinates 5, can be determined by requiring that each concentration x,(l) and x,(2) be positive. The resulting solid is a triangular prism whose shape and size depend on p. Strictly speaking, the normal form (4.5) also requires that det D # 0. The set of points within the prism not satisfying this condition is a surface whose form also depends on p. At an infinite temperature (t = 0) and zero external field (CQ = 0), the initial conditions for eq. (4.5) are obviously & = 0. It is shown in appendix A that, as the temperature is lowered, the even correlation & increases, while the odd correlations 5, and t3 remain zero until the “critical surface” det D = 0 is reached at t = t,(p) (if p >p, = l/(z - 1)). Then the solution curve bifurcates and its two branches evolve (in the ferromagnetic phase) towards two vertices of the accessible prism as t - a. If an external field is added, the solution starts from non-zero values of <,, evolves without touching the critical surface, and, as t+ m, approaches one of the two mentioned vertices. Having a solution curve in t-space, we can compute response functions in a straightforward way. The specific heat (per site in units of kB) is
(4.8) The generalized susceptibility depending on a; and r:
x,~ = S<,/a~u, is obtained
o=&($)=S,-ri D,,$ I 1
k=l
when
we regard
5, as
(4.9) I
from which we have xi, = qjlr
.
(4.10)
CLUSTER APPROACH
Here
we are interested
particularly
439
TO DILUTE MAGNETISM
in the magnetic
susceptibility
x = -x,i.
The present formulation permits to find numerically the critical temperature in a simple way: Near t,, det D is a linear function of (t - fc), changing sign at t = t,, as shown in appendix t > t, as a thermodynamically thermodynamically
stable
A. The paramagnetic-phase solution continues for unstable branch. As discussed above, the two branches,
which
correspond
to the ferromagnetic-
phase solutions, originate at t = t, from a bifurcation of the paramagnetic-phase solution. Therefore, we need a new starting point at some t, > t,. The <;(t,,) are obtained by the usual natural iteration method (see appendix B). Then we solve the system (4.5) in two directions t+ t, and t+a, thus completing the solution curve. As a comparison following remarks
of the present scheme are worth mentioning:
with the NR and NI methods, the Both NR and NI are practical
methods for the minimization of the CVM free energy at a given temperature, needing an initial choice for the set of minimization variables at this temperature. The NR method only converges for adequate choices of the initial point. The NI method, on the other hand, although is expected to converge for any initial point”), does so in a usually much larger number of iterations, specially in the ordered phase and near the critical temperature. In both methods, phase diagrams and the temperature dependence of thermodynamic functions are obtained through independent solutions of the equations for a large number of temperature values. By contrast, the present procedure is not a method for solving the CVM equations at one given temperature, but a method for following the temperature dependence of correlation and response functions through the solution of the differential equation (4.5). Here we need only one initial point at each phase. At the disordered phase this is taken to be the known infinite temperature solution. At the ordered phase any temperature 0 < T < 7, is in principle a suitable starting point. At this temperature the solution can be obtained by any of the previous methods. The trajectory of the solution is again conveniently obtained form eq. (4.5), instead of by repetition of one of the previous methods for many temperatures.
5. Results and discussion We present in this section results for a dilute ferromagnet in the fee lattice (z = 12) in the entire dilution range 0 < p < 1. The system (4.5) was solved numerically by means of the second-order Runge-Kutta method22). The matrix inversion that is required to compute the right-hand side of eq. (4.5) was performed by the usual Gaussian elimination with a maximum-pivot strategy.
P.R.C.
440
HOLVORCEM
AND
a f
2
0
Fig. 1.
R. OS6RIO
CLUSTER
0
APPROACH
TO DILUTE
441
MAGNETISM
1
05
f
2
Fig. 1. Projections of the solution trajectories in the space of correlation functions on the planes (a) t3 = 0, (b) & = 0 and (c) 5, = 0, for several concentrations p. The dashed lines are the loci of the zero-temperature results.
In fig. 1 we show the projections of some solution curves on the sides of the unit cube in t-space. In fig. 2 we display the temperature dependence of 5,) &, t3. The disordered phase is represented in figs. l(a) and l(c) as a motion along the c,-axis, starting at the origin and ending at the critical surface at 5, = &(r,) = pl(z - 1) = p/11. For the relatively high concentrations (p s 0.2) shown in these figures, the order parameters 5, and 5, nearly attain their maximum allowed values p and p* repectively, as T+ 0. As the concentration is lowered towards p = p,, however, there is an increasing influence of finite clusters of aligned spins of both signs, relative to the infinite percolating cluster, on the odd correlations ,$I and t3. This causes the values of t1 and t3 at 7 =0 to vanish for pip,. For the short-range order parameter e,, all zero-temperature clusters contribute positively, thus yielding 5, = p2 at T = 0 for all p, as shown in fig. 3. In the present (pair) approximation, the percolation threshold is well known to be given by p, = l/(z - 1) = 1 /ll. In more sophisticated approximations, p, could be determined by extrapolation of TVagainst p without special difficulties
P.R.C.
442
HOLVORCEM
AND
R. OS6RIO
a
f,
\ O.!
\
’;
I/
10
kgT/J
b p=l
C 0
kgTIJ
Fig. 2. Correlation functions versus reduced temperature for several concentrations p. In (a) 5, is the reduced magnetization. In (b) the short-range order parameter tz = (S,S,) is represented by full lines while the long-range order parameter 5, = (S:S,) IS represented by dashed lines.
CLUSTER
APPROACH
TO DILUTE
443
MAGNETISM
T 1
0 08
0.04
0
0
1
kgT/J
2
0 03
T
2 0.02
001
\
0
0
1
0.024
-
0 016
-
kg T/J
2
4 3
0
1
2 kgT/J
Fig. 3. Low-temperature correlation threshold p, = 1111 = 0.0909.
functions
versus
reduced
temperature
near
the percolation
P.R.C.
444
such as the great number possible method
approach by Cahn
of iterations
is to apply
AND
R. OS6RIO
in the natural
iteration
a low-temperature
method.
expansion
Another
following
the
and Kikuchi’?).
Figs. 4 and 5 display at low concentrations diverges
HOLVORCEM
some curves of reduced and low temperatures.
at T = 0 and 7,. The first divergence
susceptibility and specific heat The susceptibility x = - C,, /r reflects
the paramagnetic
contri-
bution of isolated spins. The reduced specific heat has at 7c a finite anomaly due to the discontinuity in the derivative of &, which becomes less visible as p decreases towards p,, as can be verified in figs. 3 and 5. We have presented explicitly, as an illustration of the method, the solution for the site-dilute Ising model in the pair approximation. It is worthwhile to remark, in conclusion, that the present scheme, based on the spin-l cluster
a
x 5-
0 15
:
0
1
2
i
2 kBT/J
Fig. 4. Low-temperature P
reduced
susceptibility
versus
reduced
temperature
for (a) p > p, and
(b)
CLUSTER
APPROACH
TO DILUTE
kgT
Fig. 5. Low-temperature tions near p,.
reduced
specific
heat versus
reduced
445
MAGNETISM
/ J
temperature
for several
concentra-
algebra of section 2, the introduction of additional constraints due to dilution as presented in section 3 and the differential equation method of solution of section 4, constitutes a practical algorithm for obtaining the thermodynamics of site-dilute Ising models for any chosen basic cluster of the cluster-variation method. Calculations for larger clusters are now in progress both for site-dilute and for spin-l Ising models.
Acknowledgements We thank Belita Koiller, Maria A. Davidovich and Mark useful discussions. This work was supported by the Brazilian and CAPES.
Appendix
0. Robbins for agencies CNPq
A
Analytic solution in the disordered phase In the paramagnetic (3.6), the configuration
phase, 5, = E3 = 0. Let 51 concentrations reduce to
5,. Then,
by eqs.
(3.5)-
P.R.C.
446
_ r
P,
x,(l)
= x2(1)
= i
x,(2)
= t(p’
- 5),
The matrix
x,(2)
AND
=x4(2)
R. OS6RIO
= i(p’+
6) 3 (A.1)
x,(2)
D becomes,
[(z -
HOLVORCEM
= x,(2)
X6(2) = (1 - p)” .
= 1 P(1 - P),
from eq. (4.6),
l)P + -11
Z
0
P(1 -P)
P(1 ZP2
0
2(P2 + O(P”
P)
0
- 5)
z
_
0
P(1 -P)
-
4P
+ S)
PWp)(p2+5)
I. (A.4
Therefore,
det D = -
2(1-P)(P’+
To get the inverse
911 =
l)Sl
Z2[P - (z -
<)‘(P”
C, we compute
p)(p’+
02(P2
LBijof D,:
some minors
Z2P(P + 5) 2(1-
(A.3)
- 5) .
ia,, = 9*, = 0 )
- ‘5) ’
(A.4) 9
=
4P-(z--1)tl
22 P2(1-P)(P2+5). At zero external 8, =
field, the system
‘y*c,,,
$2
=
a2c22
3
(4.5)
i;
=
becomes
a2c23
.
Since at infinite temperature (t = 0) we have & = 0, initially 0 and C,, = s2,1det D = 0, and therefore 5, and t3 remain The equation of motion for E2 is i*=
g$
=
(P” - 52)(P2
which can be integrated 5, = p* tanh t .
+ WP’
>
(A.3 C,, = 9,,/det D = zero as t increases.
(A.6)
to (A.7)
CLUSTER
APPROACH
TO DILUTE
MAGNETISM
447
The specific heat is now
$ (p’+
c = -t2a2~* = ;
and the susceptibility
-@a,,
,$*)(p2-
We see that x-00
P4P+ 5) and det D-+0
t, = tanh-‘(l/[(z
(&I2
)
(A.81
is
= [p - (z - l)e]
’ = det
5,) = ;
pt( 1+ p tanh t) = [l - (z - 1)p tanh t] ’
(A.9)
as t-+ t,, with
- l)p]) .
(A. 10)
At the critical temperature 5, = t2(rC) = pl(z - 1). The critical temperature 7, goes to zero at the percolation threshold p, = l/(z - 1). Below p, the only divergence in the susceptibility occurs at T = 0.
Appendix Natural
B
iteration
method
The procedure of natural iterationi’) tions of the basic cluster configurations. system (3.6):
77i =
5i
=
[$I
4ilx,(2))
71
[?I
=
1,
772 =
usually has as variables the concentraIt is therefore convenient to invert the
PJ
i=l,2,3,
(cIilxlt2),
0
(y&) = i : 1
1
1
-2 0 00-l
1
-1 1
=p*
1
(B.la)
(B.lb)
where 122121 ($Ji[) = i 1 2 1 120100
773
01 )
-1 0 0 01 . 00
P.R.C.
448
A better
HOLVORCEM
form for the constraints
4 and the column
vector
(B.la)
R. OS6RIO
is obtained
by row-reducing
the matrix
n:
rl; = P2,
rl: = c &,x,(2), /=I
AND
77; = (1 -P>’
rl; = P(1 -PI,
,
(B.3)
where 120100 .
0
(B.4)
!
000001
We perform a constrained minimization variables x,(2), 1~ 1 s 6. Defining
with A,, A,, A, as Lagrange multipliers, respect to x,(2), we obtain fflhl
+
%ICIZn
-
r,45+,,[1
-y,p,(2)7[1+
of the free energy
and equating
+ln.q(l)l-
In .xn(2)] - i
of eq. (4.1) in the
to zero its derivative
ircI,,[l +lnx,(l)l)
+:,Ai =
0.
(B.6)
i=l
The key step in the NI procedure
x,(2) = C,(T>F,(TA,, A,,
with
is to isolate
x,(2)
in brackets
in eq. (B.6): (B.7)
,
AJ[x~(~)Ix~(~)]“,
with
n*1,
(B.8)
Dn =2Y*P,(2)’ Fn(7,A,, A,, A3)= exp Due
to the
simple
structure
+:n’i)
of the
matrix
.
c$‘, some
of the
F,, are
equal,
CLUSTER
APPROACH
Fl = F2 = F4 = g, , F3 = F5 =
gi = exp(-AJy,r)
TO DILUTE
MAGNETISM
449
g:‘2, F6 = g,, where 03.9)
.
We can now use the constraints to eliminate the A-dependence of eq. (B.7). Substituting (B.7) into the constraint equations (B.3) and solving for the gi we have g, = p2/(H, + 2H2 + Hz,) = p2/(H, + 2C, + Hz,) , g2
=
Ml
-
P)lW,
g, = (1 -p)‘/H,
+
(B.lO)
WI2 3
= (1 -p)*/C,
,
with H, = C,[~~(l)Ix~(l)]~~. W e h ave used the fact that D, = D, = 0. An iteration scheme based on eq. (B.7) was used to initialize the differential equation method of the main calculation. Given a set of initial values for the correlation functions, each iteration consists of computing successively x,(l), x2(1) from eq. (3.5), H,,, gi, F,,, x,(2) = H,,F,, and the new approximations &, &, S; from eq. (B.lb). Using a tolerance of lo-* requires usually lo* to lo3 iterations to convergence. To speed this up, we have used an Aitken extrapolation formula. Although this latter procedure destroys the convergence warranties of the NI, we have in this manner safely obtained convergence in 10 iterations.
References 1) For a review, see R.B. Stinchcombe, in: Phase Transitions and Critical Phenomena, vol. 7, C. Domb and J.L. Lebowitz, eds. (Academic, London, 1983), chap. 3. 2) H. Sato, A. Arrott and R. Kikuchi, J. Phys. Chem. Solids 10 (1959) 19. 3) J.S. Smart, J. Phys. Chem. Solids 16 (1960) 169. 4) R.J. Elliott, J. Phys. Chem. Solids 16 (1960) 165. 5) See, for instance, T. Oguchi and T. Obokata, J. Phys. Sot. Japan 27 (1969) 1111. K. Moorjani and S.K. Ghatak, J. Phys. Cl0 (1977) 1027. 6) T. Kaneyoshi, I.P. Fittipaldi and H. Bayer, Phys. Stat. Sol. (b) 102 (1980) 393. 7) G.B. Taggart, Physica A 116 (1982) 34. 8) O.F. de Alclntara Bonfim and I.P. Fittipaldi, Phys. Lett. A 98 (1983) 199. 9) T. Kaneyoshi, I. Tamura and R. Honmura, Phys. Rev. B 29 (1984) 2769. 10) R. Kikuchi, Phys. Rev. 81 (1951) 988. 11) R. Kikuchi, J. Physique Colloq. 38 (1977) C7-307. 12) For a small sample of recent applications of the CVM, see refs. 2-10 and 15-16 of R. Osorio and B. Koiller, Physica A 131 (1985) 263. 13) R. Osorio, B. Koiller and M.A. Davidovich, J. Magn. Magn. Mat. 54-57 (1986) 147; Phys. Rev. B 33 (1986) 1855. 14) J.L. Moran-Lopez and J.M. Sanchez, Phys. Rev. B 33 (1986) 5059.
450
15) 16) 17) 18) 19) 20) 21) 22) 23)
P.R.C.
HOLVORCEM
AND
R. OSGRIO
J.M. Sanchez and D. de Fontaine, Phys. Rev. B 17 (1978) 2926. R. Kikuchi, J. Chem. Phys. 65 (1976) 4545. J.S. Desjardins and 0. Steinsvoll, Phys. Scripta 28 (1983) 565. M. Blume, Phys. Rev. 141 (1966) 517. H.W. Capel, Physica 32 (1966) 966. M. Blume, V.J. Emery and R.B. Griffiths, Phys. Rev. A 4 (1971) 1071. D. Mukamel and M. Blume, Phys. Rev. A 10 (1974) 10. J.A. Barker, Proc. R. Sot. (London) A 216 (1953) 45. See, for instance, W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Recipes (Cambridge Univ. Press, Cambridge, 1986). J.W. Cahn and R. Kikuchi, Acta Metall. 27 (1979) 1329.
Numerical
CLUSTER
APPROACH
Paulo Departamento
R.C.
TO DILUTE
HOLVORCEM
and Roberto
MAGNETISM OS6RIO
de Fisica, Universidade de BrasrXa, 70910 Brasilia, DF, Brazil
Received 5 February manuscript received
Revised
1988 22 April
1988
A cluster algebra is developed for the definition of independent correlation functions in the cluster-variation method (CVM) for the spin-l Ising model. A scheme is then introduced for the study of site-dilute spin-i Ising models by means of the CVM. The procedure regards the site-dilute spin-i model as the spin-l model with additional constraints due to dilution. The Desjardins-Steinsvoll algorithm is used for the transformation of the CVM equations into a set of differential equations for the independent correlation functions with the inverse temperature as parameter. The evolution of the correlation functions with temperature and the behavior of response functions such as the specific heat and the susceptibility are then obtained for any degree of dilution. As an introduction to this scheme, its detailed application is presented here for the simple case of the pair approximation.
1. Introduction In spite of the success of series-expansion and renormalization-group methods for the treatment of critical behavior in dilute magnets’), closed form approximations for these systems are still an important subject of investigation due to their potential to lead to correct phase diagrams and to the general (non-critical) behavior of thermodynamic functions with little computational effort. We are concerned in this work with the problem of quenched site dilution in the Ising model, defined in a given crystal lattice by the following Hamiltonian: H = -J
2
UjUo;‘lirlj
-
(ii)
where
J is the coupling
external
rli =
field,
1, 0,
cj = tl
C uiirli
(1.1)
>
I
constant between nearest-neighbor (nn) spins, and vi is a random variable with values
with probability with probability
0378-4371/88/$03.50 0 (North-Holland Physics
B
p , 1- p .
Elsevier Science Publishers Publishing Division)
B is an
(1.2) B.V.
P.R.C. HOLVORCEM
432
The external formally tibility. Among
taken
field will be taken into
the closed
consideration
as zero
throughout
in the discussion
form approximations
methods were the first to lead pioneering works used a two-site
AND R. OS6RIO
this paper,
but it will be
of the zero-field
for site-dilute
Ising systems,
suscepcluster
to a non-zero percolation threshold. The cluster approach in the spirit of the Bethe-
Peierls2.‘) or the constant coupling4) approximations. Several extensions of these methods have been later carried ou?). More recently, effective field theories. which use approximations based on Callen’s identity, have been extensively applied to this class of problems”~‘). The results of such calculations, however, have the limitation of depending only on the coordination number of the lattice, not taking into account its dimension or other geometric properties adequately. Kikuchi’s cluster-variation method (CVM), on the other hand, is a hierarchy of closed form approximations for the configurational entropy of classical lattice systems that take into account the geometry of the problem in a more complete way “‘-I*). When applied to pure (non-dilute) systems, the CVM results are seen to converge to the exact solutions when the size of the basic cluster is increased. In its simplest forms the CVM reduces to the mean field (one-site cluster) and Bethe-P eier . 1s ( two-site cluster) approximations. Only very recently approximations beyond the two-site cluster have been carried out in the CVM for the study of dilute systems. Osorio et al.“) have discussed the bee antiferromagnet with nn and next-nearest-neighbor interactions within the tetrahedron approximation. Moran-Lopez and Sanchez”) have given a phase diagram for the fee dilute ferromagnet also within the tetrahedron approximation and discussed the results in the context of a model of ‘He-“He mixtures with a frozen chemical short-range order parameter. The starting point of the CVM is the definition of the basic cluster, whose size limits the description of correlations between sites and determines the accuracy of the approximation. We define {a} as the set of all independent concentrations of configurations of the basic cluster. For instance, in the case of the pair approximation, the different basic cluster configurations are given by {(M, N)}, where M, N = +1 for the spin-$ Ising model. In this case, { cr} is the set of independent pair concentrations x2(M, N). Notice that the site concentrations x,(M) are not independent variables, since x,(M) = C,,,x,(M, N). The normalization C,,, x2(M, N) = 1 decreases the number of independent pair concentrations to three. The symmetry requirement interactions, decreases this x,(1, -1) =&-I, l), in the case of ferromagnetic number to two. In the canonical ensemble, the equilibrium state of the system at a given temperature T is obtained through the minimization of the free energy
CLUSTER APPROACH
TO DILUTE MAGNETISM
433
where
fTT; {a>) = E({fll) - W{(+)) > and where E( { c}) is the total energy due to the magnetic lattice and S({a}) is the “configurational entropy”,
(1.4) interactions
in the
where g( { c}) is the “statistical weight”, which is the CVM approximation for the number of equivalent arrangements of up and down spins of the whole lattice consistent with the given set of concentrations {c}. In the case of pure systems, the system of equations that follows from the minimization of the CVM free energy is usually solved by either the NewtonRaphson (NR) method15) or by the natural iteration (NI) method”). The NI uses as minimization variables the concentrations of all configurations of the basic cluster. Constraints that appear among the variables due to normalization or “translational symmetry” requirements16) are introduced by means of Lagrange multipliers. On the other hand, in the NR procedure the constraints are used to explicitly reduce the minimization variables to a set of independent “correlation functions”, which are defined by means of a convenient “cluster algebra”15). A recent interesting alternative to both the NI and the NR method is the Desjardins-Steinsvoll (DS) algorithm”), which consists of converting the problem into a set of differential equations for the correlation functions with the inverse temperature as variable. The known uncorrelated results for infinite temperature are used as initial values for an iterative procedure that leads to the thermodynamic functions at all temperatures. In this work we propose the application of the DS algorithm to the site-dilute Ising model. The independent correlation functions to be used in the DS algorithm are defined, in the next section, by means of a generalization for spin 1 of the spin-i cluster algebra. The spin-l problem is then reduced to the site-dilute spin-i model through additional constraints that describe the random distribution of non-magnetic sites. The resulting set of equations is then solved by the usual Runge-Kutta method (instead of the procedure for solving these differential equations proposed by DS). The critical temperature and thermodynamic functions such as the susceptibility and the specific heat are directly calculated as functions of the degree of dilution. This scheme is illustrated for the case of the two-site cluster approximation in the fee lattice.
434
P.R.C.
HOLVORCEM
AND
R. OSORIO
2. Cluster algebra for spin 1 In this section
the Sanchez-de
ized for the spin-l
Ising model,
Fontaine
spin-4 cluster
of which the site-dilute
algebra15)
spin-i
is general-
model is a special
case. The spin-l model, usually known as the Blume-Emery-Griffiths model and in a particular case as the Blume-Cape1 model, is of course also of interest by itself and has been applied to certain magnetic systems”), 3He-4He mixtures’“), ternary fluids2”) and other systems that display tricritical behavior. The present cluster algebra permits a straightforward definition of the correlation functions, which are the independent variables for the minimization of the CVM free energy in a given approximation of a spin-l lattice system. It also gives directly the relation between the concentration of a given cluster configuration and the correlation functions. If we denote by Sj = 0, +l the possible states of site i and by T(M, S,), M = 0, +-1, the projection operator, which takes values 1, when S; = M, and 0, when Si # M, then it is simple to verify that
T(M, Si) = (1 - M2) + $MS; + (qM* - 1)s; Let us label the clusters involved in the let n(r) be the number of sites of the specified by a set of iz(~) spin values configurations by an index 1= 1,2, . configuration
of the rth cluster
(2.1)
calculation by an index Y = 1,2, . . . and rth cluster. Each of its configurations is {M, , M,, . , M,,,,}. If we label these s 3”“‘, then the concentration of the Ith
is
Xi’-)= (T(M, >4 MM,, S,>. . . T(M,,,, , S,,,,)) ,
(2.2)
with the ( . . . ) average taken over the whole lattice. In order to find the independent variables for the CVM, we substitute the projection operator (2.1) into eq. (2.2). Thus we obtain each concentration as a linear combination of correlation Y, =
functions
m, )Yq-
y, of the form
. . . &,)P”‘)
,
(2.3)
where t is a label for the set of values {a,, a2, . . . , CQ,,}, and (Y~= 0,l or 2. The coefficients of this linear combination depend on 1 = {M, , M,, . . . , M,,,,} . We have at most 3”“’ - 1 independent correlation functions (not counting y, = 1, where we define t = 0 as the label for the set of a-values (0, 0, . . . , O}). This number corresponds to the maximum number of independent configurations of the basic cluster r = T,,, considering the normalization C, x,(T~) = 1. However, this number is usually greatly reduced due to symmetry relations
CLUSTER APPROACH
435
TO DILUTE MAGNETISM
(e.g., t = (0, 1) = t’ = {l,O}). If we label the independent correlation functions that appear in the expressions for x1(r) by an index 4 = 1,2, . . . , q,, we find an expression of the form
(2.4) To give the general expression for the z+(r, q), let us define the functions A,(M)
= 1 - M2,
A,(M)
= $M,
Then it can be shown, from eqs. (2.1)-(2.5),
A,(M)
= ;M2 - 1.
(2.5)
that
(2.6) where the sum is over the indices t of the correlation functions yI that are equivalent by symmetry operations to the same independent yq. At this point, some of the y4 may have their values constrained by parameters of the problem that are fixed beforehand. This occurs, for instance, for a quenched random distribution of Si = 0, to be discussed in the next section in the context of a site-dilute Ising model.
3. The site-dilute
Ising model in the pair approximation
As a simple application of the above algebra, we consider here a case where the basic cluster is an nn pair. We study the site-dilute Ising model with quenched, randomly distributed vacancies and ferromagnetic interactions. The clusters involved are the site (r = 1) and the nn pair (r = 2). This level of the CVM hierarchy is equivalent to the Bethe-Peierls approximation and is expected to be exact in the central region of a Cayley tree. The configurational entropy (per site in units of k,) can be written as
s =
3/l
l$
P,(l)ax,(l)l+
3/2 ,$
PIwam)l
3
(3.1)
with Z(X) = x In (x). Here /3,(r ) is the degeneracy of the Zth configuration of the rth cluster, i.e., the number of configurations that are symmetry-equivalent to each independent one labeled as I, and the 7,. are geometrical factors introduced by Barker*l). For a coordination number z, we have x = -z/2 and
436
P.R.C.
HOLVORCEM
AND
R. OS6RIO
TABLE I
Configurations
of the site (Y = 1) and the pair (r = 2) clusters and their degeneracies P,(r),
M
I
P,(l)
(M> N) (1,‘) (1, -1)
0 _
1 1 1 _
_
_
_
_
1
1 2 3 4 5
-1
6
P,(2) :. 2
(1.0) (-1, -1) (-1.0) (07 0)
: 1
y1 = z - 1. The degeneracies of the three configurations of the site and of the six configurations of the pair are listed in table I. The independent spin-l correlations that occur in the expansion of eq. (2.4) for r=1,2 are
and the coefficients I+(Y, s) are listed in table II. (Note that in tables I and II we use M, N instead of M,, M2.)Since in the dilute model the vacancies are quenched with concentration 1 - p, we have the additional constraints
y2 =
(s:) =p,
so that there
5, =y,, From
Y,
are only three
independent
52= Y,l
eqs. (2.4),
variables
(3.3) left:
53=y4.
(3.2)-(3.4)
x,(l) = i(P + 5,)> Coefficients
=
and table
(3.4) II, one gets
x,(l) = t(P - t,>,
+(1)=1-p,
TABLE II v,(I, y) of the expansion (2.4) of the cluster configurations the correlation functions y, of eq. (3.2).
(3.5)
x!(r) in terms
4
v,( 1) 4)
GL 4)
0
1 - M’ Ml2 3M’i2 - 1 0 0 0
(1 - M’)(l - N’) (M/2)(1 - NZ) + (N/2)(1 - MZ) (3M’/2 - l)(lN’) + (3N’i2l)(l ~ M’) MN14 (M/2)(3N’/2 ~ 1) + (N/2)(3M’/2 - 1) (3ML/2 - 1)(3N’i2 - 1)
1 2 3 4 5
of
CLUSTER APPROACH
TO DILUTE MAGNETISM
il it.
437
(3.6)
5, P21
Let us call A the above matrix of coefficients Al,q+l = ~~(2, q) (with l= 1 ,*.., 6; q = 0, 1, . . . , 5), which is used in the differential equation method of solution of the next section.
4. Differential equation method of solution Our solution of the model follows the general method proposed by Desjardins and Steinsvollr7), which converts the problem of minimization of the free energy into an initial-value problem of ordinary differential equations. For a fixed spin concentration p, the correlation functions si depend on the temperature. The free energy (per site in units of the coupling constant J) can be written in terms of the l-variables as
(4.1)
where CY,= - BIJ, CY~= -z/2, CQ= 0,, 7 = k,TIJ, and s depends on ti through eqs. (3.1), (3.5) and (3.6). The condition of minimum free energy, df/dg,= 0, is written as (Yi/T = as/&&
(4.2)
)
which is differentiated
with respect to the parameter
t = 1 /r to yield
(4.3)
where ii = d$ldt
.
(4.4)
P.R.C.
438
This differential ki=C
equation
Ciiai,
HOLVORCEM
AND
R. OS6RIO
can be put in normal
form,
i-1,2,3,
(4.5)
j=l
where
Cij = (D-l),.
For our model
the D, are
(4.6) with
Knij = Pn(2)(‘;,An,+ ‘iAn
+ ‘j,An,)(‘j,An, + ‘jzA.4 + ‘j4.5) >
(4.7)
where S, is the Kronecker delta. The accessible region in the space of coordinates 5, can be determined by requiring that each concentration x,(l) and x,(2) be positive. The resulting solid is a triangular prism whose shape and size depend on p. Strictly speaking, the normal form (4.5) also requires that det D # 0. The set of points within the prism not satisfying this condition is a surface whose form also depends on p. At an infinite temperature (t = 0) and zero external field (CQ = 0), the initial conditions for eq. (4.5) are obviously & = 0. It is shown in appendix A that, as the temperature is lowered, the even correlation & increases, while the odd correlations 5, and t3 remain zero until the “critical surface” det D = 0 is reached at t = t,(p) (if p >p, = l/(z - 1)). Then the solution curve bifurcates and its two branches evolve (in the ferromagnetic phase) towards two vertices of the accessible prism as t - a. If an external field is added, the solution starts from non-zero values of <,, evolves without touching the critical surface, and, as t+ m, approaches one of the two mentioned vertices. Having a solution curve in t-space, we can compute response functions in a straightforward way. The specific heat (per site in units of kB) is
(4.8) The generalized susceptibility depending on a; and r:
x,~ = S<,/a~u, is obtained
o=&($)=S,-ri D,,$ I 1
k=l
when
we regard
5, as
(4.9) I
from which we have xi, = qjlr
.
(4.10)
CLUSTER APPROACH
Here
we are interested
particularly
439
TO DILUTE MAGNETISM
in the magnetic
susceptibility
x = -x,i.
The present formulation permits to find numerically the critical temperature in a simple way: Near t,, det D is a linear function of (t - fc), changing sign at t = t,, as shown in appendix t > t, as a thermodynamically thermodynamically
stable
A. The paramagnetic-phase solution continues for unstable branch. As discussed above, the two branches,
which
correspond
to the ferromagnetic-
phase solutions, originate at t = t, from a bifurcation of the paramagnetic-phase solution. Therefore, we need a new starting point at some t, > t,. The <;(t,,) are obtained by the usual natural iteration method (see appendix B). Then we solve the system (4.5) in two directions t+ t, and t+a, thus completing the solution curve. As a comparison following remarks
of the present scheme are worth mentioning:
with the NR and NI methods, the Both NR and NI are practical
methods for the minimization of the CVM free energy at a given temperature, needing an initial choice for the set of minimization variables at this temperature. The NR method only converges for adequate choices of the initial point. The NI method, on the other hand, although is expected to converge for any initial point”), does so in a usually much larger number of iterations, specially in the ordered phase and near the critical temperature. In both methods, phase diagrams and the temperature dependence of thermodynamic functions are obtained through independent solutions of the equations for a large number of temperature values. By contrast, the present procedure is not a method for solving the CVM equations at one given temperature, but a method for following the temperature dependence of correlation and response functions through the solution of the differential equation (4.5). Here we need only one initial point at each phase. At the disordered phase this is taken to be the known infinite temperature solution. At the ordered phase any temperature 0 < T < 7, is in principle a suitable starting point. At this temperature the solution can be obtained by any of the previous methods. The trajectory of the solution is again conveniently obtained form eq. (4.5), instead of by repetition of one of the previous methods for many temperatures.
5. Results and discussion We present in this section results for a dilute ferromagnet in the fee lattice (z = 12) in the entire dilution range 0 < p < 1. The system (4.5) was solved numerically by means of the second-order Runge-Kutta method22). The matrix inversion that is required to compute the right-hand side of eq. (4.5) was performed by the usual Gaussian elimination with a maximum-pivot strategy.
P.R.C.
440
HOLVORCEM
AND
a f
2
0
Fig. 1.
R. OS6RIO
CLUSTER
0
APPROACH
TO DILUTE
441
MAGNETISM
1
05
f
2
Fig. 1. Projections of the solution trajectories in the space of correlation functions on the planes (a) t3 = 0, (b) & = 0 and (c) 5, = 0, for several concentrations p. The dashed lines are the loci of the zero-temperature results.
In fig. 1 we show the projections of some solution curves on the sides of the unit cube in t-space. In fig. 2 we display the temperature dependence of 5,) &, t3. The disordered phase is represented in figs. l(a) and l(c) as a motion along the c,-axis, starting at the origin and ending at the critical surface at 5, = &(r,) = pl(z - 1) = p/11. For the relatively high concentrations (p s 0.2) shown in these figures, the order parameters 5, and 5, nearly attain their maximum allowed values p and p* repectively, as T+ 0. As the concentration is lowered towards p = p,, however, there is an increasing influence of finite clusters of aligned spins of both signs, relative to the infinite percolating cluster, on the odd correlations ,$I and t3. This causes the values of t1 and t3 at 7 =0 to vanish for pip,. For the short-range order parameter e,, all zero-temperature clusters contribute positively, thus yielding 5, = p2 at T = 0 for all p, as shown in fig. 3. In the present (pair) approximation, the percolation threshold is well known to be given by p, = l/(z - 1) = 1 /ll. In more sophisticated approximations, p, could be determined by extrapolation of TVagainst p without special difficulties
P.R.C.
442
HOLVORCEM
AND
R. OS6RIO
a
f,
\ O.!
\
’;
I/
10
kgT/J
b p=l
C 0
kgTIJ
Fig. 2. Correlation functions versus reduced temperature for several concentrations p. In (a) 5, is the reduced magnetization. In (b) the short-range order parameter tz = (S,S,) is represented by full lines while the long-range order parameter 5, = (S:S,) IS represented by dashed lines.
CLUSTER
APPROACH
TO DILUTE
443
MAGNETISM
T 1
0 08
0.04
0
0
1
kgT/J
2
0 03
T
2 0.02
001
\
0
0
1
0.024
-
0 016
-
kg T/J
2
4 3
0
1
2 kgT/J
Fig. 3. Low-temperature correlation threshold p, = 1111 = 0.0909.
functions
versus
reduced
temperature
near
the percolation
P.R.C.
444
such as the great number possible method
approach by Cahn
of iterations
is to apply
AND
R. OS6RIO
in the natural
iteration
a low-temperature
method.
expansion
Another
following
the
and Kikuchi’?).
Figs. 4 and 5 display at low concentrations diverges
HOLVORCEM
some curves of reduced and low temperatures.
at T = 0 and 7,. The first divergence
susceptibility and specific heat The susceptibility x = - C,, /r reflects
the paramagnetic
contri-
bution of isolated spins. The reduced specific heat has at 7c a finite anomaly due to the discontinuity in the derivative of &, which becomes less visible as p decreases towards p,, as can be verified in figs. 3 and 5. We have presented explicitly, as an illustration of the method, the solution for the site-dilute Ising model in the pair approximation. It is worthwhile to remark, in conclusion, that the present scheme, based on the spin-l cluster
a
x 5-
0 15
:
0
1
2
i
2 kBT/J
Fig. 4. Low-temperature P
reduced
susceptibility
versus
reduced
temperature
for (a) p > p, and
(b)
CLUSTER
APPROACH
TO DILUTE
kgT
Fig. 5. Low-temperature tions near p,.
reduced
specific
heat versus
reduced
445
MAGNETISM
/ J
temperature
for several
concentra-
algebra of section 2, the introduction of additional constraints due to dilution as presented in section 3 and the differential equation method of solution of section 4, constitutes a practical algorithm for obtaining the thermodynamics of site-dilute Ising models for any chosen basic cluster of the cluster-variation method. Calculations for larger clusters are now in progress both for site-dilute and for spin-l Ising models.
Acknowledgements We thank Belita Koiller, Maria A. Davidovich and Mark useful discussions. This work was supported by the Brazilian and CAPES.
Appendix
0. Robbins for agencies CNPq
A
Analytic solution in the disordered phase In the paramagnetic (3.6), the configuration
phase, 5, = E3 = 0. Let 51 concentrations reduce to
5,. Then,
by eqs.
(3.5)-
P.R.C.
446
_ r
P,
x,(l)
= x2(1)
= i
x,(2)
= t(p’
- 5),
The matrix
x,(2)
AND
=x4(2)
R. OS6RIO
= i(p’+
6) 3 (A.1)
x,(2)
D becomes,
[(z -
HOLVORCEM
= x,(2)
X6(2) = (1 - p)” .
= 1 P(1 - P),
from eq. (4.6),
l)P + -11
Z
0
P(1 -P)
P(1 ZP2
0
2(P2 + O(P”
P)
0
- 5)
z
_
0
P(1 -P)
-
4P
+ S)
PWp)(p2+5)
I. (A.4
Therefore,
det D = -
2(1-P)(P’+
To get the inverse
911 =
l)Sl
Z2[P - (z -
<)‘(P”
C, we compute
p)(p’+
02(P2
LBijof D,:
some minors
Z2P(P + 5) 2(1-
(A.3)
- 5) .
ia,, = 9*, = 0 )
- ‘5) ’
(A.4) 9
=
4P-(z--1)tl
22 P2(1-P)(P2+5). At zero external 8, =
field, the system
‘y*c,,,
$2
=
a2c22
3
(4.5)
i;
=
becomes
a2c23
.
Since at infinite temperature (t = 0) we have & = 0, initially 0 and C,, = s2,1det D = 0, and therefore 5, and t3 remain The equation of motion for E2 is i*=
g$
=
(P” - 52)(P2
which can be integrated 5, = p* tanh t .
+ WP’
>
(A.3 C,, = 9,,/det D = zero as t increases.
(A.6)
to (A.7)
CLUSTER
APPROACH
TO DILUTE
MAGNETISM
447
The specific heat is now
$ (p’+
c = -t2a2~* = ;
and the susceptibility
-@a,,
,$*)(p2-
We see that x-00
P4P+ 5) and det D-+0
t, = tanh-‘(l/[(z
(&I2
)
(A.81
is
= [p - (z - l)e]
’ = det
5,) = ;
pt( 1+ p tanh t) = [l - (z - 1)p tanh t] ’
(A.9)
as t-+ t,, with
- l)p]) .
(A. 10)
At the critical temperature 5, = t2(rC) = pl(z - 1). The critical temperature 7, goes to zero at the percolation threshold p, = l/(z - 1). Below p, the only divergence in the susceptibility occurs at T = 0.
Appendix Natural
B
iteration
method
The procedure of natural iterationi’) tions of the basic cluster configurations. system (3.6):
77i =
5i
=
[$I
4ilx,(2))
71
[?I
=
1,
772 =
usually has as variables the concentraIt is therefore convenient to invert the
PJ
i=l,2,3,
(cIilxlt2),
0
(y&) = i : 1
1
1
-2 0 00-l
1
-1 1
=p*
1
(B.la)
(B.lb)
where 122121 ($Ji[) = i 1 2 1 120100
773
01 )
-1 0 0 01 . 00
P.R.C.
448
A better
HOLVORCEM
form for the constraints
4 and the column
vector
(B.la)
R. OS6RIO
is obtained
by row-reducing
the matrix
n:
rl; = P2,
rl: = c &,x,(2), /=I
AND
77; = (1 -P>’
rl; = P(1 -PI,
,
(B.3)
where 120100 .
0
(B.4)
!
000001
We perform a constrained minimization variables x,(2), 1~ 1 s 6. Defining
with A,, A,, A, as Lagrange multipliers, respect to x,(2), we obtain fflhl
+
%ICIZn
-
r,45+,,[1
-y,p,(2)7[1+
of the free energy
and equating
+ln.q(l)l-
In .xn(2)] - i
of eq. (4.1) in the
to zero its derivative
ircI,,[l +lnx,(l)l)
+:,Ai =
0.
(B.6)
i=l
The key step in the NI procedure
x,(2) = C,(T>F,(TA,, A,,
with
is to isolate
x,(2)
in brackets
in eq. (B.6): (B.7)
,
AJ[x~(~)Ix~(~)]“,
with
n*1,
(B.8)
Dn =2Y*P,(2)’ Fn(7,A,, A,, A3)= exp Due
to the
simple
structure
+:n’i)
of the
matrix
.
c$‘, some
of the
F,, are
equal,
CLUSTER
APPROACH
Fl = F2 = F4 = g, , F3 = F5 =
gi = exp(-AJy,r)
TO DILUTE
MAGNETISM
449
g:‘2, F6 = g,, where 03.9)
.
We can now use the constraints to eliminate the A-dependence of eq. (B.7). Substituting (B.7) into the constraint equations (B.3) and solving for the gi we have g, = p2/(H, + 2H2 + Hz,) = p2/(H, + 2C, + Hz,) , g2
=
Ml
-
P)lW,
g, = (1 -p)‘/H,
+
(B.lO)
WI2 3
= (1 -p)*/C,
,
with H, = C,[~~(l)Ix~(l)]~~. W e h ave used the fact that D, = D, = 0. An iteration scheme based on eq. (B.7) was used to initialize the differential equation method of the main calculation. Given a set of initial values for the correlation functions, each iteration consists of computing successively x,(l), x2(1) from eq. (3.5), H,,, gi, F,,, x,(2) = H,,F,, and the new approximations &, &, S; from eq. (B.lb). Using a tolerance of lo-* requires usually lo* to lo3 iterations to convergence. To speed this up, we have used an Aitken extrapolation formula. Although this latter procedure destroys the convergence warranties of the NI, we have in this manner safely obtained convergence in 10 iterations.
References 1) For a review, see R.B. Stinchcombe, in: Phase Transitions and Critical Phenomena, vol. 7, C. Domb and J.L. Lebowitz, eds. (Academic, London, 1983), chap. 3. 2) H. Sato, A. Arrott and R. Kikuchi, J. Phys. Chem. Solids 10 (1959) 19. 3) J.S. Smart, J. Phys. Chem. Solids 16 (1960) 169. 4) R.J. Elliott, J. Phys. Chem. Solids 16 (1960) 165. 5) See, for instance, T. Oguchi and T. Obokata, J. Phys. Sot. Japan 27 (1969) 1111. K. Moorjani and S.K. Ghatak, J. Phys. Cl0 (1977) 1027. 6) T. Kaneyoshi, I.P. Fittipaldi and H. Bayer, Phys. Stat. Sol. (b) 102 (1980) 393. 7) G.B. Taggart, Physica A 116 (1982) 34. 8) O.F. de Alclntara Bonfim and I.P. Fittipaldi, Phys. Lett. A 98 (1983) 199. 9) T. Kaneyoshi, I. Tamura and R. Honmura, Phys. Rev. B 29 (1984) 2769. 10) R. Kikuchi, Phys. Rev. 81 (1951) 988. 11) R. Kikuchi, J. Physique Colloq. 38 (1977) C7-307. 12) For a small sample of recent applications of the CVM, see refs. 2-10 and 15-16 of R. Osorio and B. Koiller, Physica A 131 (1985) 263. 13) R. Osorio, B. Koiller and M.A. Davidovich, J. Magn. Magn. Mat. 54-57 (1986) 147; Phys. Rev. B 33 (1986) 1855. 14) J.L. Moran-Lopez and J.M. Sanchez, Phys. Rev. B 33 (1986) 5059.
450
15) 16) 17) 18) 19) 20) 21) 22) 23)
P.R.C.
HOLVORCEM
AND
R. OSGRIO
J.M. Sanchez and D. de Fontaine, Phys. Rev. B 17 (1978) 2926. R. Kikuchi, J. Chem. Phys. 65 (1976) 4545. J.S. Desjardins and 0. Steinsvoll, Phys. Scripta 28 (1983) 565. M. Blume, Phys. Rev. 141 (1966) 517. H.W. Capel, Physica 32 (1966) 966. M. Blume, V.J. Emery and R.B. Griffiths, Phys. Rev. A 4 (1971) 1071. D. Mukamel and M. Blume, Phys. Rev. A 10 (1974) 10. J.A. Barker, Proc. R. Sot. (London) A 216 (1953) 45. See, for instance, W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Recipes (Cambridge Univ. Press, Cambridge, 1986). J.W. Cahn and R. Kikuchi, Acta Metall. 27 (1979) 1329.
Numerical