PACS Numbers: 05.50.+q
PACS Numbers: 05.50.+q
PACS Numbers: 05.50.+q
4 AUGUST 2003
S. Ranjbar∗
I. INTRODUCTION
As mentioned in a previous paper [1], the first realistic model to be solved exactly
and which exhibited a phase transition was the two-dimensional Ising model with nearest-
neighbor interactions in a zero external field. The exact expression for the partition function
of this model was first derived by Onsager [2] in one of the most celebrated and much quoted
papers in mathematical physics. The surprising feature of Onsager’s solution at the time
was that the model exhibited a phase transition with a symmetric logarithmic divergence
of the specific heat, which was, of course, at variance with the predictions of the classical
theories. It was known from the work of Peierls [3], that the model had a phase transition,
but prior to Onsager, various approximation schemes due to Bragg and Williams [4,5], Bethe
[6], Peierls [3,7], and others all predicted classical critical behavior. One of the purposes of
trying to find exact solutions in equilibrium statistical mechanics is, of course, to provide
a testing ground for approximate methods; in so doing one hopes that some insight will be
gained into the intricacies of phase transitions.
Since Onsager’s solution of the two-dimensional Ising model in a zero field, numerous
man-hours of effort have been spent and many a good researcher has fallen foul to the
dreaded ‘Ising’s disease’, a failure to yield an exact solution to the two-dimensional model
in non-zero field and the three-dimensional model, even with a zero field.
In the previous work [1], although we could not find an analytical expression for a
square lattice in the presence of the magnetic field, yet we were able to produce useful
results. In this work, our aim is the investigation of the two-dimensional triangular Ising
model in the absence and presence of a magnetic field, using the transfer matrix method. In
this model, which is solved exactly in the absence of a magnetic field [8], each lattice point
interacts with six of its neighbors. The triangular lattice has an order-disorder transition
of a type very similar to that of the rectangular lattice. The reader can refer to the first
reference for more information about the Ising model and its applications.
The original idea of the transfer matrix method was due to Kramers and Wannier [9]
and it, in fact, formed the basis for Onsager’s solution of the two-dimensional model. It
has since proved to be useful for formulating and analyzing a number of model systems in
statistical mechanics.
In a general setting let us suppose that we have a lattice composed of l ‘layers’, which
could, for example, consist of points, lines, or planes for one-, two-, and three dimensional
lattices, respectively. If we suppose that the ‘state’ of the ith layer is specified by a quantity
σi , which can take only a finite number of discrete values, and that interactions occur only
within layers and between nearest neighboring layers, the interaction energy E{σ} in a given
configuration {σ} = (σ1 , σ2 , . . . , σl ) or state of the lattice can, in general, be expressed in
the form [10]
l
X
E{σ} = U (σi ) + V (σi , σi+1 ) , (1)
i=1
where U (σi ) represents the interaction energy within a layer in state σ i and V (σi , σi+1 )
represents the interaction energy between nearest neighboring layers in states σ i and σi+1 .
For simplicity we have assumed periodic boundary conditions, that is, we assume that lth
layer is coupled to the first, or equivalently, that σ l+1 = σ1 .
The partition function for a lattice model whose interaction energy can be expressed
in the form of Eq. (1) is given by
l
!
X X X
Zl = exp(−βE(σ)) = exp −β (U (σi ) + V (σi , σi+1 )) ,
{σ} (σ1 ,...,σl ) i=1
(2)
X l
Y
= T (σi , σi+1 ) ,
(σ1 ,...,σl ) i=1
VOL. 41 S. RANJBAR 385
where T2 (σ1 , σ3 ) denotes the (σ1 , σ3 ) component of the matrix T2 . It can also be noted
that the matrix T in effect transfers us from one layer to the next, hence the term transfer
matrix is applicable.
By repeating the above process and summing successively over the layer states
σ2 , . . . , σl in Eq. (2) we then obtain
X
Zl = Tl (σ1 , σ1 ) = Tr(Tl ) , (5)
σ
where Tr denotes the trace of the matrix, i.e. the sum of the diagonal entries.
Using now the fact that the trace of a matrix is also the sum of its eigenvalues and
that the eigenvalues of Tl are λli , i = 1, 2, . . . , M , where the λi s are eigenvalues of T, we
can write Eq. (5) as follows:
M
X
Zl = λli . (6)
i=1
Since the elements T (σ, σ 0 ) are all positive, it follows from the Perron-Frobenius [11]
theorem, that the maximum eigenvalue is both non-degenerate and an analytic function of
its arguments and, in particular, from Eq [3], an analytic function of β or the temperature
(T 6= 0). Without loss of generality, we can then order the eigenvalues such that
taking λmax to be an analytic function of β. If the number of layers goes to infinity, then
from Eqs. (6) and (7), the free energy can be obtained as
M
" l #
l
X λi
A = −kT ln Zl = −kT ln λmax (1 + ) = −lkT ln λmax . (8)
λmax
i=2
In the following section, we illustrate the transfer matrix method by evaluating the
limiting free energy for the one-dimensional Ising model with nearest-neighbor interactions
in an external magnetic field.
386 INVESTIGATION OF TWO-DIMENSIONAL TRIANGULAR . . . VOL. 41
where the cB term represents the magnetic dipole energy in the applied field B and J is
the nearest-neighbor coupling constant. As in Eq. (1), we have imposed periodic boundary
conditions by coupling the last spin to the first, or equivalently, setting σ N +1 ≡ σ1 .
In the notation of the previous section, the layers are simply lattice sites on the chain;
the number of states per layer (M ) is two, corresponding to spin up σ = +1 or spin down
σ = −1. The interaction energy between the nearest-neighbor layers is given by
The transfer matrix T is then a 2 × 2 matrix with entries given from Eq. (3) by
where
m = ej , h = eb . (14)
IV. SOLVING THE ISING MODEL FOR THE TRIANGULAR PLANAR LATTICE
FIG. 2: The simplest triangular Ising model with three rows and an infinite number of columns.
FIG. 3: The modified triangular lattice model with three rows. Each line connecting two spins
shows the interacting nearest neighbor spins.
TABLE I: All Possible Configurations for Each Layer of FIG. 3 Along With Their Configurational
Energies.
State A1 A2 A3 A4 A5 A6 A7 A8
u e u u u e e e
J
J
J
J
J
J
J
J
u
J u u
J u e
J u u
J e e
J e u
J e e
J u e
J e
Layer
By solving the secular determinant, |T 3 − λI| = 0, in the absence of magnetic field, we may
3/N
obtain the maximum eigenvalue , λmax = Z3 = z33 , as
where
s √
3 1 − 3x + 1 − 6x − 3x2
f3 (x) = (20)
2
and
m4 − 1
x= , (21)
(m4 + 1)2
in which z3 is the partition function for each spin in the 3-row model.
Now, we can extend our model by adding a row to the model which is shown in FIG.
2. Again, we face the problem that the number of nearest neighbor spins in such a 4-row
model is different, depending on the location of the spin. To eliminate such differences, and
to have a model in which each spin has six nearest neighbors, the model must be modified
[1]. For this model, by including the periodic boundary condition, the secular determinant
has an order of 16. In the absence of any magnetic field, the maximum eigenvalue, z 4 , can
be obtained as
where
s p
4 a1 + a 2 +(a1 + a2 )2 − 32x4
f4 (x) = ,
4
(23)
a1 = 1 − 4x + 2x2 ,
p
a2 = (1 − 4x + 2x2 )2 + 16x2 (2x − 1) .
It is evident that by adding two or three rows to FIG. 2 and doing the same modifi-
cation as previously mentioned, we can obtain models with five and six rows, respectively.
The exact partition functions for each spin in the absence of any magnetic field in five and
six row models can be obtained as:
where
s p
5 b1 + b 2 +(b1 + b2 )2 − 80x4
f5 (x) = ,
4
(25)
b1 = 1 − 5x + 5x2 ,
q
b2 = b21 − 20x2 (1 − x)2 + 40x3 ,
and
where
s p
6 c1 +c21 − 8x6
f6 (x) = ,
2
1 − 6x + 9x2 − 2x3 + (1 − 2x)c2 + c3
c1 = , (27)
p 4
c2 = x4 + 8x3 + 2x2 − 8x + 1 ,
p
c3 = 2(4x2 + 4x − 1)(7x4 − 10x2 + 8x − 1 − (x2 − 4x + 1)c2 ) .
It is possible to extend the model with successive increments in the number of rows.
However, for n = 7, 8, and 9 the secular determinant, in the absence of any magnetic field,
gives a polynomial with an order of 8, 18, and 23, respectively. Since such polynomials
cannot be solved analytically, we are not able to obtain a mathematical expression for
the partition function. However, the polynomials and implicit derivatives can be solved
numerically to obtain the exact partition functions, free energies (−A/N kT = ln z), internal
energies (−E/N J = ∂ ln z/∂j) and heat capacities (C/N k = j 2 ∂ 2 ln z/∂j 2 ).
All calculations were carried out on a personal computer, using the well-known soft-
ware program “Maple”. Unfortunately, we were not able to solve for the secular determinant
VOL. 41 S. RANJBAR 391
FIG. 4: The exact reduced free energy, −A/N kT , for the models with n = 3 − 11 compared to
that for the square lattice (n → ∞). The top curve (dashed) is for n = 3, the bottom curve is for
n → ∞, and the other curves are for n = 4 − 11 from top to the bottom, respectively.
FIG. 5: Same as in FIG. 3, but for the reduced internal energy (the results for the model with
n = 10 and 11 are not included).
392 INVESTIGATION OF TWO-DIMENSIONAL TRIANGULAR . . . VOL. 41
for n > 9 analytically, mainly due to the very long computational time required. However,
the determinant was solved numerically only for n = 10 and 11, from which we calculated
the exact free energy. The results of such calculations are given in FIGs. 4-6 for the free
energy, the internal energy, and heat capacity, respectively. The exact thermodynamic
properties for the triangular planar lattice (n → ∞) are also shown in these figures [12,13].
By adding one row to our model, the order of the transfer matrix will increase two-
fold. Therefore, when n increases, the solution of the secular determinant becomes more
complicated. This means that the model cannot be solved for the triangular planar lattice
(n → ∞) because of the computer limitations.
However, we may refer to results for which the analytical expression for the free
energy of the model has been obtained. The z n is given in Eqs. (19), (22), (24), and (26)
for n = 3, 4, 5, and 6, respectively. These equations have a common term, 2e j cosh 2j, from
which we may expect that zn will include such a term regardless of the value of n. In
addition, we can see in FIG. 5 that all the curves, irrespective of the number of the rows
at a specific temperature, Jc = 1/4 ln 3, pass through a fixed point −E/N J = 2, because
∂ ln f (x) ∂x e−j ∂ ln f
= = 2 (1 − 2 tanh 2j) . (28)
∂f ∂f cosh 2j ∂f
This is the zero at jc , the critical temperature of the triangular lattice. Thus, it is
VOL. 41 S. RANJBAR 393
possible to give the following general form for the free energy of the model:
Now, we have to obtain a mathematical expression for f (n, x). For this purpose, we
expand the functions f (3, x) to f (6, x). By keeping in mind that the value of x is between
zero and 0.125, it is expected that the expansions of these functions will rapidly converge.
The expansions of these functions as a power series in x, can be obtained as
20 3 82 4
f3 = 1 − x − 2x2 − x − x + O(x5 ) , (30)
3 3
2 3 127 4 641 5
f4 = 1 − x − 2x − 7x − x − x + O(x6 ) , (31)
4 4
829 5 4631 6
f5 = 1 − x − 2x2 − 7x3 − 32x4 − x − x + O(x7 ) , (32)
5 5
5597 6 33181 7
f6 = 1 − x − 2x2 − 7x3 − 32x4 − 166x5 − x − x + O(x8 ) . (33)
6 6
According to Eqs. (30-33), the following expressions are proposed for f 7 and f8
38765 7
f7 = 1 − x − 2x2 − 7x3 − 32x4 − 166x5 − 933x6 − x + O(x8 ) , (34)
7
f8 = 1 − x − 2x2 − 7x3 − 32x4 − 166x5 − 933x6 − 5538x7 + O(x8 ) . (35)
Thus, according to Eqs. (30-35), we suggest a polynomial for the model with any number
of rows:
∞
X
f (n, x) = 1 − tk,n xk , (36)
k=1
where the exact values of the coefficients t k,n in the power series can be determined from
the coefficient of xk at the k + 1 row model. Therefore, for n > 8 the coefficients x to x 7 in
Eq. (36) are the same as the coefficients x to x 7 in Eq. (35). Now, a question that arises is:
how can we determine the other coefficients for n > 8? If a mathematical relation among
the tk ’s existed then, we could show them as an analytical function of k. Unfortunately,
such a relation does not exist, so we must approximate them.
As mentioned previously, our aim is to obtain a mathematical expression for the
partition function with any number of rows. In this way, our studies show that the following
expression is appropriate for the f (n, x):
sn
X 8k−1
f (n, x) ' 1 − xk . (37)
k2
k=1
The coefficient sn can be evaluated by fitting the exact values of f (n, x) into Eq. (37) for
different values of n. The next step is to obtain the value of the s n from Eq. (37) for the
394 INVESTIGATION OF TWO-DIMENSIONAL TRIANGULAR . . . VOL. 41
FIG. 7: Comparison of the heat capacity calculated from Eq. (39) (dotted curve) with the exact
value (solid curve) for n = 4 and 9 and n → ∞.
square lattice (n → ∞). To do so, we have fitted the value s n versus n. The approximated
mathematical expression for the sn (except for n = 3) is found to be
Now we can use Eqs. (29), (37), and (38) to calculate the thermodynamic properties of a
triangular lattice with any number of rows. Using these equations, we have
sn
!
1X 8k k
ln zn = j + ln 2 cosh 2j + ln 1 − x . (39)
8 k2
k=1
Eq. (39) can be used to calculate the internal energy (−E/N J = ∂ ln z/∂j) and the
heat capacity (C/N k = j 2 ∂ 2 ln z/∂j 2 ). In order to investigate the accuracy of the proposed
free energy given by Eq. (39), we have compared its heat capacity with the exact value in
FIG. 7 for n = 4, 9 and n → ∞. As shown, we may conclude that the heat capacity, and
hence the proposed free energy given by Eq. (39), becomes more accurate when n becomes
larger, the case which we are most interested in.
Using Eq. (39), in FIG. 8 we consider the behavior of the heat capacity for n = 50.
As shown, the singularity of the heat capacity can be observed only for n → ∞ (Onsager
transition [14]). In other words, for a model with a limited number of rows we may expect
to observe only a diffuse transition [14] or a continuous transition [15].
VOL. 41 S. RANJBAR 395
FIG. 8: The reduced heat capacity versus the inverse of the reduced temperature, j/kT , calculated
from Eq. (39) for n = 50 (dotted curve) and n → ∞ (solid curve) at B = 0.
In the presence of a magnetic field, the Ising model cannot be solved exactly in two-
dimensions. Here, also we do not claim to solve a triangular Ising model in the presence of
a magnetic field. Instead, the approach presented for the two-dimensional Ising model [1]
is carried out step by step for the triangular model.
The exact approach given in the previous sections for solving the model in the absence
of any external magnetic field can be used to obtain the exact thermodynamic properties
of a model in the presence of a magnetic field. For example, Eq. (18) can be used to obtain
the exact partition function for a model with n = 3. As explained before, for the case that
B = 0, we can set up a transfer matrix with larger n values, as well. Unfortunately, due
to algebraic limitations, irrespective of n = 3, we cannot derive the maximum eigenvalues
analytically and, subsequently, we cannot introduce an analytical equation for the partition
function such as Eq. (39). Thus, for the model with n = 4, 5, and 6, we calculate the
partition function by the numerical method. The calculated partition function is used
to obtain the thermodynamic properties of the model. For instance, the result for the
heat capacities are shown in FIG. 9 for a given number of rows at specified magnetic field
strengths, where b = cB/kT . We have also used the exact partition function to calculate
the magnetization per spin, I,
∂ ln z
I= . (40)
∂b T
396 INVESTIGATION OF TWO-DIMENSIONAL TRIANGULAR . . . VOL. 41
FIG. 9: The exact reduced heat capacity versus the inverse of the reduced temperature for a given
magnetic field for n = 3 (solid curve), n = 4 (dotted curve), n = 5 (dashed curve), n = 6 (dashed
dotted curve).
The results of such calculations are shown in FIG. 10 for some given number of rows at
specified magnetic field strengths.
VII. CONCLUSION
We have used a lattice model with a limited number of rows, n, to calculate the
exact thermodynamic properties of the triangular Ising lattice by using the transfer matrix
method. The exact free energy was calculated only for the case where n ≤ 11, while
the exact internal energy and heat capacity were calculated for n ≤ 9. The calculated
free energy was fitted into a polynomial expression. We assumed that the free energy of
the model is given by a polynomial whose degree depends on the number of rows in the
model (see Eq (39)). This approach gives more accurate results when n becomes larger.
Expression Eq. (39), has been used to calculate the thermodynamic properties of the planar
lattice (n → ∞) (see FIG. 7). As shown in the figure, our calculated heat capacity (derived
from Eq. (39)) is in very good agreement with the exact one.
According to Eq. (39), a singular point for the heat capacity is expected to be observed
only for a model with an unlimited number of rows (n → ∞).
We have shown how to set up the transfer matrix in the presence of an external
magnetic field (see Eq. (18) for n = 3). The solution to such a matrix equation gives
the exact solution for the partition function of the model in the presence of a magnetic
VOL. 41 S. RANJBAR 397
FIG. 10: Same as FIG. 9, for the magnetization for n = 4 (solid curve), n = 5 (dotted curve), n = 6
(dashed curve).
field. The matrix can be written for the model with any number of rows. However, due
to the mentioned limitations, it can be solved only for the model with a limited number of
rows. In fact, the expression obtained in this case is more complicated than that in which
B = 0, for any value of n. For this reason, unlike the previous case, we have not been able
to give a general analytical expression for the partition function when B 6= 0. However,
the exact heat capacity and magnetization given in FIG.’s 9 and 10 reveal an important
conclusion: the calculated results for n = 3, 4, 5, and 6 are almost the same when b ' j or
b > j. To emphasize this point more clearly, we have plotted the exact values of the heat
capacity for n = 6 against that for n = 5 in FIG. 11, when b = j. From such results, we
conclude that our exact calculated values, for the model with n = 4, 5, 6, or 7, are almost
the same as those for the actual model (n → ∞), when b ' j or b > j. For the case where
B = 0, the main difference between the thermodynamic properties of a model with different
values of n appears around the singular point of the heat capacity (see FIG. 6). For the
case when B 6= 0, in which there is no singularity, we may expect that such differences
become insignificant, especially for stronger magnetic fields. This is in accordance with our
conclusion cited above.
Acknowledgment
I would like to thank the Razi University, Kermanshah-Iran for financial support.
398 INVESTIGATION OF TWO-DIMENSIONAL TRIANGULAR . . . VOL. 41
FIG. 11: The exact reduced heat capacity for n = 6 versus that for n = 5 when b = j.
References
∗
Electronic address: [email protected]
[1] S. Ranjbar and G. A. Parsafar, J. Phys. Chem. B. 103, 7514 (1999).
[2] L. Onsager, Phys. Rev. 65, 117 (1944).
[3] R. Peierls, Proc. Camb. Philos. Soc. 32, 477 (1963).
[4] W. L. Bragg and E. J. Williams, Proc. R. Soc. London A 145, 699 (1934).
[5] W. L. Bragg and E. J. Williams, Proc. R. Soc. London A 151, 540 (1935).
[6] Bethe, H. A. Proc. R. Soc. London A 150, 552 (1935).
[7] R. Peierls, Proc. R. Soc. London A 154, 207 (1936).
[8] C. Domb, The Critical Point, (Taylor and Francis, London, 1996).
[9] H. A. Kramers and G. H. Wannier, I-II Phys. Rev. 60, 252 (1941).
[10] C. J. Thompson, Classical Equilibrium Statistical Mechanics, (Clarendon Press, Oxford, 1988).
[11] F. R. Gantmacher, The Theory of Matrices, (Chelsea, New York, 1964),Vol 2.
[12] G. F. Newell, Phys. Rev. 79, 876 (1950).
[13] G. H. Wannier, Phys. Rev. 79, 357 (1950).
[14] A. Múnster, Statistical Thermodynamics (Springer-Verlag, Berlin, 1969),Vol.I.
[15] L. E. Reichl, A modern course in statistical physics, (Wiley, New York, 1998).