- Email: [email protected]
Chiral order is a property of a directed bond
Volume 150, number 5,6,7
PHYSICS LETTERS A
12 November 1990
Chiral order is a property of a directed bond T. Morita and S. Fujiki Department of EngineeringScience, Facultyof Engineering, Tohoku University, Sendai 980, Japan Received 26 July 1989; accepted for publication 12 September 1990 Communicated by A.A. Maradudin
It is argued that the chiral order in a magnetic system is a property of a directed bond on a lattice, and the chiral order is discussed in terms of the distribution function for a pair of lattice sites. This fact is illustrated by the calculation of the chiral order in the triangle approximation in the cluster variation method, for the six- and twelve-clockmodels on the triangular lattice.
The chiral order in a magnetic system is a property characteristic to a plaquette and hence properties o f a plaquette are required in order to discuss the chiral order [ 1 ]. The purpose o f the present note is to call attention to the fact that the chiral order can also be regarded as a property o f a bond. We consider a lattice on which lattice sites are connected by bonds in such a way that any two o f the bonds do not intersect with each other and the total n u m b e r o f bonds connected to a lattice site is even. The space in the lattice is divided into plaquettes by the bonds. We now consider such a state in the chiral order for this lattice, that every plaquette has chirality up or down and the chiralities o f the neighbouring plaquettes which have a c o m m o n b o n d on the boundaries are opposite. In order to describe this state, we draw an arrow on every bond in such a way that the bonds surrounding a plaquette o f chirality up are directed counterclockwise and those with chirality down clockwise. We call a plaquette counterclockwise (clockwise) if the bonds surrounding it are directed counterclockwise (clockwise). We now label the spins on the sites surrounding a plaquette by Sl, s2, ..., st along the direction on the bonds, where l is the total n u m b e r o f sites or bonds surrounding the plaquette. The chirality of the plaquette is defined, up to a normalization constant, by
Cpl ~
~pl
Si × $i + i
!
,
(1)
where st+~ is understood to denote s~, the brackets denote the thermodynamic average, and ypj is 1 ( - 1 ) if the plaquette is counterclockwise (clockwise). If all the spins s~ are restricted within the xy-plane, the chirality is up (down) if cp~ is along the positive (negative) z-axis. Expression ( 1 ) shows that the chirality can be calculated if we know (s~×si+~) for the l bonds. We now define the chirality of a directed bond ij with a direction from site i to j by
eij= (si Xsj ) .
(2)
In order to calculate this we need to know the pair distribution function for the directed bond. In the typical case that the lattice is regular and all the c~j are equal, we only have to calculate this quantity for a directed bond. As an example we take the n-clock model on a triangular lattice. In this model every spin has magnitude one and can have n directions within the xyplane. The Hamiltonian is given by
H=J ~ cos(0,-0j), (~j)
(3)
where 0~ is the spin variable taking values 0, 2n/n, ..., 2 ( n - 1 )n/n, and J i s assumed to be positive, i.e., we consider here the antiferromagnetic case. The summation in (3) is taken over the pairs o f nearest neighbour sites. In order to discuss the chiral order, we put an arrow representing the direction to every b o n d connecting a pair o f nearest neighbour lattice
0375-9601/90/$ 03.50 © 1990 - Elsevier Science Publishers B.V. (North-Holland)
299
Volume 150, number 5,6,7
PHYSICS LETTERS A
sites, in such a way that every upward plaquette is counterclockwise and every downward one clockwise, as shown in fig. I. We can calculate the distribution function for a directed bond, employing the cluster variation method, which is a systematic approximation method of calculating the distribution functions. In order to determine the distribution function for a directed bond in the present case, we have to use the triangle approximation, and write the distribution functions for a site, for a directed bond and for a triangle. The lattice is divided into three sublattices, and we have distribution functions for the clusters shown in fig. 2. They are denoted by p} t ) (0) and po, 0' ) for ~(2)(0, i= 1, 2, 3, and p~(O, 0', 0"). Here and in the following, i, j and k denote l, 2 and 3 in cyclic order. The calculation is made by using two different methods. In the first method, the distribution functions in the present approximation are expressed in the following form [ 2 ]:
p(~')( O) =exp [flF,,(') +A(~t)( O) ] , i= 1, 2, 3,
(4a)
12 November 1990
i= 1, 2, 3 ,
(4b)
p123~(0, 0 ' , / 9 " ) = e x p { f l F ( 3 ) - k
+ cos(0' - 0" ) +cos(0" - 0) ]
+AI~(O, 0', 0" )},
(4c)
where fl= 1/kBT, kBis the Boltzmann constant, T is the temperature, and
A}t)(O)=3262)(O)+32}2)(O),
(5a)
+,~,(~) (o) + ~j~2)(o' ) + 2 ~ 3) (0, o' ), i = 1 , 2 , 3,
(5b)
Al3~(O, 0', 0" ) =211 ) (0) +2~)(0' ) + ~ ) ( 0 " ) +21~ )(0) + ~ ) ( 0 '
) +~?,)(0" )
+~)(o',o")+~,)(o",o)+2l~)(o,o')
.
(5c)
The constants F} t), F/~2) and F (3) are determined by the normalization conditions of the distribution functions, and 2~z)(O), 2~)(0) and 2~3)(0, 0') for i=1, 2, 3 and for 0, 0'=0, 2x/n, ..., 2 ( n - 1 ) x / n , by the reducibility conditions: 0'
=exp[flF~ 2) + f l J c o s ( 0 - 0' ) +A~2)(0, 0') l ,
i= 1, 2, 3,
Ab2)(O, 0') = 2 [2~z) (0) +2j(,2)(0') l
p,(')(0)= Y'.p~Z)(O,O'),
pb~)(o,o ' )
flJ[cos( O-O ' )
p ~ ' ) ( 0 ' ) = Zp~2)(O,O'),
i=1,2,3, i=1,2,3,
(6a) (6b)
0
pl~)(O, 0 ' ) = ZPI~(O, 0', 0"),
(7a)
p~)(O', 0") = ~'. p}~(O, 0', 0"),
(7b)
p~)(O", 0)= ~pt]](O, 0', 0").
(7c)
0"
0
O'
The magnetization of a sublattice i, m~, is given by
mi= 2~ cosOp!~)(O) , 0
Fig. 1. Triangular lattice with directed bonds.
(8)
and the chirality of a bond ij, c~., that is the z-component of e~j, is calculated by .___. 2
.___. 3
c~= 2~ Y. sin(O-O')p~:)(O,O').
3
2
1
(9)
0 O'
2
The free energy per site f is given by 1
2
3
1
2
Fig. 2. Clusters which are incorporated in the triangle approximation.
300
3
f= ~ i=1
~.3--i/1 17(1) ----01~(2))+ 2 F ( 3 )
(10)
Volume 150, n u m b e r 5,6,7
P H Y S I C S LETTERS A
In the low-temperature limit, the solution should have the three-sublattice antiferromagnetic symmetry [3]:
p~)(O)=p)l)(OT2n/3), p~2) ( 0 , 0 ' ) = p j ~ ) ( O ~
i = 1 , 2 , 3,
12 N o v e m b e r 1990
q.O sin[a-ta)
(lla)
......
m,
2n/3, fiT- 2 n / 3 ) ,
i=1,2,3,
(lib)
and
pl~(o, o', o") =pt~(O"T-2n/3, 0~ 2~/3, 0'-T-2~/3),
where, if 0 + 2 g / 3 > 2n or 0 - 2 ~ / 3 <0, the variable 0 denotes O-2n or 0+2n. ( 6 a ) - ( 7 c ) were solved, assuming this symmetry, by iteration of the following 2n+ n 2 quantities: g, (0) =NI exp[21zZ)(0) ] ,
I
o.o
(llc)
0.0
I
0.1
I
I
0.3
1 0.5
11
I 07
k~T/J
Fig. 3. Magnetization of a site m~ and ehiral order of a directed bond c~ plotted against temperature for the antiferromagnetic nclock models with n = 6 and 12. Dashed lines at low temperature show the asymptotes of the curves for the plane rotator model (n = ~ ) in the same approximation.
g'~ (0) =N~ exp [2t32) (0) ] ,
g2(O, 0') =Nz exp[212)(0) +2~2) (O' ) +2t23)(0, 0 ' ) ] , where N~, NI and Nz denote the normalization constants ofgl (O),g'~(0) and gz(O, O'), respectively. We choose initial values for the quantities, put them in the right-hand sides of (6a), (6b) and (7a), and calculate the new values and their changes from the previous values. We repeat the procedure by using the new values for the initial values, until the sum of the absolute values of the changes is less than 10 -6 , in the calculation for n = 6. In the calculation at a temperature in the ordered and the paramagnetic phase, the final values in a lower and a higher temperature, respectively, were used as the initial values. For the system for n = 6, the ordered-phase solution persists up to kB T/J= 0.6 122, but the free energy for the paramagnetic-phase solution becomes lower at kB T/J =0.61 13, and hence this is the temperature of the first-order transition. The obtained m~ and c~ are plotted in fig. 3 for n = 6 and also for n = 12. The accuracy of the results for n = 12 is not as good as for n = 6. The asymptotes at the low-temperature limit of the curves for n = oo are also shown. In the second method, the free energy in the present approximation per site f i s expressed as follows:
3
f=J E Z E cos(O-o')pb~)( o, o') i=1 0 O' 3
- T Z (]S~ 1) -S~2)) -2TS(3) ,
(12)
i= 1
where S~ ~J = -kB ~ p~l)(0) l n p ~ l ) ( 0 ) ,
(13a)
8
S~2)=-kBY. ~pbZ)(O,O')lnpb2)(O,O'),
(13b)
0 a'
S~3)=-kn E Z Zpl~(O,O',O") 0 0' O"
×lnp123](0, 0', 0 " ) .
(13c)
We regard (12) as a functional ofpl3~(0, 0', 0") by determining p ~ ) ( 0 ) and p~2)(0, 0') by ( 6 a ) - ( 7 c ) from it. By the original variational principle of the cluster variation method, we get the distribution functions and the free energy f by minimizing this functional with respect to the variations of p l~(O, 0', O"), under the normalization condition. As stated in ref. [4], we can exploit the knowledge that thepl,3~(0, 0', 0" ) to be obtained takes the form (4c) with (5c), so we put
301
Volume 150, number 5,6,7
PHYSICS LETTERS A
pt~(o,o',o") =exp [flF~3)+2~ (0) +22(0' ) + 2 3 ( 0 " ) + ~ 12(0, 0' ) + 2 2 3 ( 0 ' , 0" ) +~31 (0", 0) ] ,
(14)
where F ~3) is the normalization constant. Assuming the symmetry ( 11 c) and expressing 2 j (0) as a Fourier series, [n/2l
;t,(O)=½ao+
Y'. amcosmO
[(n-l)/2]
+
~
bmsinmO,
(15)
m=l
and 2 t2 (0, 0' ) in a similar fashion, as a double Fourier series with respect to 0and 0', we have n~r-n unknown parameters to be determined in expression (14). Here [ n / 2 ] denotes the greatest integer not exceeding n/2. By determining the unknowns so as to make expression (12) minimum, we confirmed the above results for n = 6. The more accurate value of the first-order-transition temperature is 0.6113119, and the values of rn~ and c,~ at the transition temperature are 0.2623 and 0.1107, respectively. In order to confirm that expression (12) is minimum, the Hessian matrix of this quantity as a function of the n 2_ n parameters is calculated and all its eigenvalues are confirmed to be positive. As the temperature is lowered for the paramagnetic-phase solution, a doubly degenerate eigenvalue, which corresponds to the three-sublattice antiferromagnetic order, becomes
302
12 November 1990
negative below 0.6088808. At further lowering the temperature, one more eigenvalue, corresponding to the chiral order without magnetization, becomes negative below 0.5552625. In conclusion, we see that the chirality can be calculated when knowing the correlation function for a directed bond, provided that the directions are put on bonds in the lattice in such a way that all the bonds surrounding every plaquette are directed either counterclockwise or clockwise. Two methods of calculation in the triangle approximation of the cluster variation method are sketched for a system on the triangular lattice. In the cluster variation method, the lowest approximation that enables us to discuss the chiral order, which is a property determined by the distribution function of a directed bond, is the triangle approximation given in the present note. References [ 1 ] J. Villain, J. Phys. C 10 (1977) 1717, 4793; S. Miyashita and H. Shiba, J. Phys. Soc. Japan 53 (1984) 1145; S. Fujiki and D.D. Betts, Prog. Theor. Phys. Suppl. No. 87 (1986) 268; F. Matsubara and S. Inawashiro, Solid State Commun. 67 (1988) 229; N. Kawashima and M. Suzuki, preprint (1989). [2] T. Morita, J. Math. Phys. 13 (1972) 115. [3] S. Katsura, T. Ide and T. Morita, J. Star. Phys. 42 (1986) 381, [4] T. Morita, Phys. Lett. A 138 (1989) 485.
PHYSICS LETTERS A
12 November 1990
Chiral order is a property of a directed bond T. Morita and S. Fujiki Department of EngineeringScience, Facultyof Engineering, Tohoku University, Sendai 980, Japan Received 26 July 1989; accepted for publication 12 September 1990 Communicated by A.A. Maradudin
It is argued that the chiral order in a magnetic system is a property of a directed bond on a lattice, and the chiral order is discussed in terms of the distribution function for a pair of lattice sites. This fact is illustrated by the calculation of the chiral order in the triangle approximation in the cluster variation method, for the six- and twelve-clockmodels on the triangular lattice.
The chiral order in a magnetic system is a property characteristic to a plaquette and hence properties o f a plaquette are required in order to discuss the chiral order [ 1 ]. The purpose o f the present note is to call attention to the fact that the chiral order can also be regarded as a property o f a bond. We consider a lattice on which lattice sites are connected by bonds in such a way that any two o f the bonds do not intersect with each other and the total n u m b e r o f bonds connected to a lattice site is even. The space in the lattice is divided into plaquettes by the bonds. We now consider such a state in the chiral order for this lattice, that every plaquette has chirality up or down and the chiralities o f the neighbouring plaquettes which have a c o m m o n b o n d on the boundaries are opposite. In order to describe this state, we draw an arrow on every bond in such a way that the bonds surrounding a plaquette o f chirality up are directed counterclockwise and those with chirality down clockwise. We call a plaquette counterclockwise (clockwise) if the bonds surrounding it are directed counterclockwise (clockwise). We now label the spins on the sites surrounding a plaquette by Sl, s2, ..., st along the direction on the bonds, where l is the total n u m b e r o f sites or bonds surrounding the plaquette. The chirality of the plaquette is defined, up to a normalization constant, by
Cpl ~
~pl
Si × $i + i
!
,
(1)
where st+~ is understood to denote s~, the brackets denote the thermodynamic average, and ypj is 1 ( - 1 ) if the plaquette is counterclockwise (clockwise). If all the spins s~ are restricted within the xy-plane, the chirality is up (down) if cp~ is along the positive (negative) z-axis. Expression ( 1 ) shows that the chirality can be calculated if we know (s~×si+~) for the l bonds. We now define the chirality of a directed bond ij with a direction from site i to j by
eij= (si Xsj ) .
(2)
In order to calculate this we need to know the pair distribution function for the directed bond. In the typical case that the lattice is regular and all the c~j are equal, we only have to calculate this quantity for a directed bond. As an example we take the n-clock model on a triangular lattice. In this model every spin has magnitude one and can have n directions within the xyplane. The Hamiltonian is given by
H=J ~ cos(0,-0j), (~j)
(3)
where 0~ is the spin variable taking values 0, 2n/n, ..., 2 ( n - 1 )n/n, and J i s assumed to be positive, i.e., we consider here the antiferromagnetic case. The summation in (3) is taken over the pairs o f nearest neighbour sites. In order to discuss the chiral order, we put an arrow representing the direction to every b o n d connecting a pair o f nearest neighbour lattice
0375-9601/90/$ 03.50 © 1990 - Elsevier Science Publishers B.V. (North-Holland)
299
Volume 150, number 5,6,7
PHYSICS LETTERS A
sites, in such a way that every upward plaquette is counterclockwise and every downward one clockwise, as shown in fig. I. We can calculate the distribution function for a directed bond, employing the cluster variation method, which is a systematic approximation method of calculating the distribution functions. In order to determine the distribution function for a directed bond in the present case, we have to use the triangle approximation, and write the distribution functions for a site, for a directed bond and for a triangle. The lattice is divided into three sublattices, and we have distribution functions for the clusters shown in fig. 2. They are denoted by p} t ) (0) and po, 0' ) for ~(2)(0, i= 1, 2, 3, and p~(O, 0', 0"). Here and in the following, i, j and k denote l, 2 and 3 in cyclic order. The calculation is made by using two different methods. In the first method, the distribution functions in the present approximation are expressed in the following form [ 2 ]:
p(~')( O) =exp [flF,,(') +A(~t)( O) ] , i= 1, 2, 3,
(4a)
12 November 1990
i= 1, 2, 3 ,
(4b)
p123~(0, 0 ' , / 9 " ) = e x p { f l F ( 3 ) - k
+ cos(0' - 0" ) +cos(0" - 0) ]
+AI~(O, 0', 0" )},
(4c)
where fl= 1/kBT, kBis the Boltzmann constant, T is the temperature, and
A}t)(O)=3262)(O)+32}2)(O),
(5a)
+,~,(~) (o) + ~j~2)(o' ) + 2 ~ 3) (0, o' ), i = 1 , 2 , 3,
(5b)
Al3~(O, 0', 0" ) =211 ) (0) +2~)(0' ) + ~ ) ( 0 " ) +21~ )(0) + ~ ) ( 0 '
) +~?,)(0" )
+~)(o',o")+~,)(o",o)+2l~)(o,o')
.
(5c)
The constants F} t), F/~2) and F (3) are determined by the normalization conditions of the distribution functions, and 2~z)(O), 2~)(0) and 2~3)(0, 0') for i=1, 2, 3 and for 0, 0'=0, 2x/n, ..., 2 ( n - 1 ) x / n , by the reducibility conditions: 0'
=exp[flF~ 2) + f l J c o s ( 0 - 0' ) +A~2)(0, 0') l ,
i= 1, 2, 3,
Ab2)(O, 0') = 2 [2~z) (0) +2j(,2)(0') l
p,(')(0)= Y'.p~Z)(O,O'),
pb~)(o,o ' )
flJ[cos( O-O ' )
p ~ ' ) ( 0 ' ) = Zp~2)(O,O'),
i=1,2,3, i=1,2,3,
(6a) (6b)
0
pl~)(O, 0 ' ) = ZPI~(O, 0', 0"),
(7a)
p~)(O', 0") = ~'. p}~(O, 0', 0"),
(7b)
p~)(O", 0)= ~pt]](O, 0', 0").
(7c)
0"
0
O'
The magnetization of a sublattice i, m~, is given by
mi= 2~ cosOp!~)(O) , 0
Fig. 1. Triangular lattice with directed bonds.
(8)
and the chirality of a bond ij, c~., that is the z-component of e~j, is calculated by .___. 2
.___. 3
c~= 2~ Y. sin(O-O')p~:)(O,O').
3
2
1
(9)
0 O'
2
The free energy per site f is given by 1
2
3
1
2
Fig. 2. Clusters which are incorporated in the triangle approximation.
300
3
f= ~ i=1
~.3--i/1 17(1) ----01~(2))+ 2 F ( 3 )
(10)
Volume 150, n u m b e r 5,6,7
P H Y S I C S LETTERS A
In the low-temperature limit, the solution should have the three-sublattice antiferromagnetic symmetry [3]:
p~)(O)=p)l)(OT2n/3), p~2) ( 0 , 0 ' ) = p j ~ ) ( O ~
i = 1 , 2 , 3,
12 N o v e m b e r 1990
q.O sin[a-ta)
(lla)
......
m,
2n/3, fiT- 2 n / 3 ) ,
i=1,2,3,
(lib)
and
pl~(o, o', o") =pt~(O"T-2n/3, 0~ 2~/3, 0'-T-2~/3),
where, if 0 + 2 g / 3 > 2n or 0 - 2 ~ / 3 <0, the variable 0 denotes O-2n or 0+2n. ( 6 a ) - ( 7 c ) were solved, assuming this symmetry, by iteration of the following 2n+ n 2 quantities: g, (0) =NI exp[21zZ)(0) ] ,
I
o.o
(llc)
0.0
I
0.1
I
I
0.3
1 0.5
11
I 07
k~T/J
Fig. 3. Magnetization of a site m~ and ehiral order of a directed bond c~ plotted against temperature for the antiferromagnetic nclock models with n = 6 and 12. Dashed lines at low temperature show the asymptotes of the curves for the plane rotator model (n = ~ ) in the same approximation.
g'~ (0) =N~ exp [2t32) (0) ] ,
g2(O, 0') =Nz exp[212)(0) +2~2) (O' ) +2t23)(0, 0 ' ) ] , where N~, NI and Nz denote the normalization constants ofgl (O),g'~(0) and gz(O, O'), respectively. We choose initial values for the quantities, put them in the right-hand sides of (6a), (6b) and (7a), and calculate the new values and their changes from the previous values. We repeat the procedure by using the new values for the initial values, until the sum of the absolute values of the changes is less than 10 -6 , in the calculation for n = 6. In the calculation at a temperature in the ordered and the paramagnetic phase, the final values in a lower and a higher temperature, respectively, were used as the initial values. For the system for n = 6, the ordered-phase solution persists up to kB T/J= 0.6 122, but the free energy for the paramagnetic-phase solution becomes lower at kB T/J =0.61 13, and hence this is the temperature of the first-order transition. The obtained m~ and c~ are plotted in fig. 3 for n = 6 and also for n = 12. The accuracy of the results for n = 12 is not as good as for n = 6. The asymptotes at the low-temperature limit of the curves for n = oo are also shown. In the second method, the free energy in the present approximation per site f i s expressed as follows:
3
f=J E Z E cos(O-o')pb~)( o, o') i=1 0 O' 3
- T Z (]S~ 1) -S~2)) -2TS(3) ,
(12)
i= 1
where S~ ~J = -kB ~ p~l)(0) l n p ~ l ) ( 0 ) ,
(13a)
8
S~2)=-kBY. ~pbZ)(O,O')lnpb2)(O,O'),
(13b)
0 a'
S~3)=-kn E Z Zpl~(O,O',O") 0 0' O"
×lnp123](0, 0', 0 " ) .
(13c)
We regard (12) as a functional ofpl3~(0, 0', 0") by determining p ~ ) ( 0 ) and p~2)(0, 0') by ( 6 a ) - ( 7 c ) from it. By the original variational principle of the cluster variation method, we get the distribution functions and the free energy f by minimizing this functional with respect to the variations of p l~(O, 0', O"), under the normalization condition. As stated in ref. [4], we can exploit the knowledge that thepl,3~(0, 0', 0" ) to be obtained takes the form (4c) with (5c), so we put
301
Volume 150, number 5,6,7
PHYSICS LETTERS A
pt~(o,o',o") =exp [flF~3)+2~ (0) +22(0' ) + 2 3 ( 0 " ) + ~ 12(0, 0' ) + 2 2 3 ( 0 ' , 0" ) +~31 (0", 0) ] ,
(14)
where F ~3) is the normalization constant. Assuming the symmetry ( 11 c) and expressing 2 j (0) as a Fourier series, [n/2l
;t,(O)=½ao+
Y'. amcosmO
[(n-l)/2]
+
~
bmsinmO,
(15)
m=l
and 2 t2 (0, 0' ) in a similar fashion, as a double Fourier series with respect to 0and 0', we have n~r-n unknown parameters to be determined in expression (14). Here [ n / 2 ] denotes the greatest integer not exceeding n/2. By determining the unknowns so as to make expression (12) minimum, we confirmed the above results for n = 6. The more accurate value of the first-order-transition temperature is 0.6113119, and the values of rn~ and c,~ at the transition temperature are 0.2623 and 0.1107, respectively. In order to confirm that expression (12) is minimum, the Hessian matrix of this quantity as a function of the n 2_ n parameters is calculated and all its eigenvalues are confirmed to be positive. As the temperature is lowered for the paramagnetic-phase solution, a doubly degenerate eigenvalue, which corresponds to the three-sublattice antiferromagnetic order, becomes
302
12 November 1990
negative below 0.6088808. At further lowering the temperature, one more eigenvalue, corresponding to the chiral order without magnetization, becomes negative below 0.5552625. In conclusion, we see that the chirality can be calculated when knowing the correlation function for a directed bond, provided that the directions are put on bonds in the lattice in such a way that all the bonds surrounding every plaquette are directed either counterclockwise or clockwise. Two methods of calculation in the triangle approximation of the cluster variation method are sketched for a system on the triangular lattice. In the cluster variation method, the lowest approximation that enables us to discuss the chiral order, which is a property determined by the distribution function of a directed bond, is the triangle approximation given in the present note. References [ 1 ] J. Villain, J. Phys. C 10 (1977) 1717, 4793; S. Miyashita and H. Shiba, J. Phys. Soc. Japan 53 (1984) 1145; S. Fujiki and D.D. Betts, Prog. Theor. Phys. Suppl. No. 87 (1986) 268; F. Matsubara and S. Inawashiro, Solid State Commun. 67 (1988) 229; N. Kawashima and M. Suzuki, preprint (1989). [2] T. Morita, J. Math. Phys. 13 (1972) 115. [3] S. Katsura, T. Ide and T. Morita, J. Star. Phys. 42 (1986) 381, [4] T. Morita, Phys. Lett. A 138 (1989) 485.