Effect of biquadratic exchange on phase transitions of a planar classical Heisenberg ferromagnet
M.Žukovič, T.Idogaki and K.Takeda
Institute of Environmental Systems, Kyushu University
Department of Applied Quantum Physics, Kyushu University
Abstract. Effect of biquadratic exchange on phase transitions of a planar classical
Heisenberg (or XY) ferromagnet on a stacked triangular lattice is investigated by Standard Monte Carlo
and Histogram Monte Carlo simulations in the region of a bilinear to biquadratic exchange interaction
ratio . The biquadratic exchange is found to cause separate secondorder phase
transitions in a strong biquadratic exchange limit, followed by simultaneous dipole and quadrupole
ordering, which is of first order for an intermediate range of the exchange ratio and changes to a
secondorder one again as is further increased. Thus, a phase diagram featuring both
triple and tricritical points is obtained. Furthermore, a finitesize scaling analysis is used to
calculate the critical indices for both dipole and quadrupole kinds of ordering.
: 75.10.Hk; 75.30.Kz; 75.40.Cx; 75.40.Mg.
: Planar Heisenberg ferromagnet; Biquadratic exchange; Phase transition; Quadrupole ordering;
Multicritical point; Histogram Monte Carlo simulation;
Corresponding author.
Permanent address: Institute of Environmental Systems, Faculty of Engineering, Kyushu University,
Fukuoka 8128581, Japan
Tel.: +81926423811; Fax: +81926336958
Email: milantap@mbox.nc.kyushuu.ac.jp
1.Introduction
It has long been recognized that there are materials in which also biquadratic (or generally
multipolar) interactions may play an important role as in some cases they are comparable or even much
larger than the bilinear ones. The problem of higherorder interactions in systems with Heisenberg
symmetry has been tackled in several mean field approximation (MFA) studies
sivar72 ()chenlevy1 (), by hightemperature series expansion (HTSE) calculations
chenlevy2 (), as well as within a framework of some other approximative schemes
micnas (); chaddha99 (). Those studies have shown that such interactions can induce various interesting
properties such as tricritical and triple points, quadrupole ordering, separate dipole and quadrupole
phase transitions etc. Much less attention, however, has been paid to this problem on systems with XY
spin symmetry. Chen chenetal1 (); chenetal2 () calculated transition temperatures and the
susceptibility critical indices for a planar Heisenberg ferromagnet with biquadratic exchange on cubic
lattices by the HTSE method. However, they only investigated the region of ratio for which
MFA predicts a secondorder transition and hence, they failed to provide a complete phase diagram with
all possible phase transitions and their nature for this model. For example, since MFA predicts
firstorder transitions for dipole ordering if , only quadrupole phase transitions,
which are expected to be of second order carmesin (), were investigated. Here we note that, in spite
of all the previous conclusions based on various approaches, the rigorous proof of the existence of
dipole longrange order (DLRO), corresponding to the ferromagnetic directional arrangement of spins, and
quadrupole longrange order (QLRO), representing an axially ordered state in which spins can point in
either direction along the axis of ordering, at finite temperature on the classical bilinearbiquadratic
exchange model has only recently been provided independently by Tanaka and Idogaki tanaka (), and
Campbell and Chayes campbell ().
We have recently performed MC simulations on the XY model with the bilinearbiquadratic exchange
Hamiltonian on a simple cubic lattice in order to establish phase boundaries and describe critical
behaviour of the system nagata (). Unfortunately, we did not succeed in finding a conclusive answer
to the problem of the order of the DLRO transition in the region of the separate DLRO and QLRO
transitions. In the present paper we chose the lattice with a hexagonal (stacked triangular) symmetry,
increased the lattice size and performed a careful finitesize scaling (FSS) analysis in order to obtain
more reliable data and deliver a complete picture of a longrange ordering and its nature for the
considered system via Standard Monte Carlo (SMC) and Histogram Monte Carlo (HMC) simulations. The phase
diagram obtained captures all important features induced by the biquadratic exchange such as separate
ordering of dipoles and quadrupoles, the appearance of triple and tricritical points, shortrange order
above the quadrupole ordering temperature, etc. Furthermore, we perform a FSS analysis to determine the
order of a transition from the energy cumulant scaling, and to calculate the critical indices for both
DLRO and QLRO transitions.
2.Model and simulation technique
We are concerned with the planar classical Heisenberg model with bilinear and biquadratic
interactions, described by the Hamiltonian
(1) 
where is a twodimensional unit vector at the th lattice
site , denotes the sum over nearest neighbors, and are the
bilinear and biquadratic exchange interaction constants, respectively.
We first perform SMC simulations on systems of spins, where = 12, 15, 18, 24, 30,
assuming periodic boundary condition throughout. For a fixed value of the exchange ratio
we run simulations starting at low (high) temperatures from a ferromagnetic (random) initial
configuration and gradually raise (lower) temperature. These heatingcooling loops serve to check
possible hysteresis, accompanying firstorder transitions. As we move in
space, we use the last spin configuration as an input for calculation at the next point. We sweep
through the spins in sequence and updating follows a Metropolis dynamics. In the updating process, the
new direction of spin in the spin flip is selected completely at random, without any limitations by a
maximum angle of spin rotation or allowed discrete set of resulting angle values. Thermal averages are
calculated using at most Monte Carlo steps per spin (MCS/s) after equilibrating over
another MCS/s. We calculate the system internal energy and some other physical
quantities defined as follows:
the specific heat per site
(2) 
the DLRO parameter (magnetization) ,
(3) 
the QLRO parameter ,
(4) 
and the following quantities which are functions of the parameter () :
the susceptibility per site
(5) 
the logarithmic derivatives of and with respect to
(6) 
(7) 
the fourthorder longrange order (LRO) cumulants and
(8) 
and the fourthorder energy cumulant
(9) 
These quantities are used to estimate the nature of a transition as well as its position. Firstorder
transitions usually manifest themselves by discontinuities in the order parameter and energy, and
hysteresis when cooling and heating. If transition is second order, it can be approximately localized by
the peak position and also doublechecked by the position of the fourthorder LRO and energy
cumulants curves intersection for various lattice sizes.
In the next stage we perform HMC calculations, developed by Ferrenberg and Swendsen
ferrswen1 (); ferrswen2 (), at the transition temperatures estimated from the SMC calculations for
each lattice size. Here, MCS are used for calculating averages after discarding another
MCS for thermalization. We calculate the energy histogram , the order parameters
histograms , as well as the physical quantities (2)(9). Using data
from the histograms one can calculate physical quantities at neighboring temperatures, and thus
determine the values of extrema of various quantities and their locations with high precision for each
lattice size. In such a way we can obtain quality data for FSS analysis which determines the order of
the transition and, in the case of a secondorder transition, it also allows us to extract critical
indices. For example, the energy cumulant exhibits a minimum near critical temperature ,
which achieves the value in the limit for a secondorder
transition, while is expected for a firstorder transition
ferrswen1 (); ferrswen2 (). The extrema of a variety of thermodynamic quantities at a secondorder
transition are known to scale with a lattice size as, for example:
(10) 
(11) 
(12) 
where and represent the correlation length and susceptibility critical
exponents, respectively. In the case of a firstorder transition (except for the order parameters), they
display a volumedependent scaling, . The simulations were performed on the vector
supercomputer FUJITSU VPP700/56.
3.FSS analysis and phase diagram
We started our calculations in a region where the biquadratic exchange is equal to the bilinear
one, i.e. , and gradually lowered the ratio down to 0. For the ratio greater than
we found only one phase transition corresponding to a dipole longrange
ordering which is of second order. Hence, as far as phase transitions are concerned, in this range of
the exchange ratio the presence of the biquadratic exchange is only reflected in the position of the
transition temperature (as well as critical indices to some extent) but causes no qualitative changes
compared with the case of guillou (); janke (). The secondorder character of a transition in
this region was also confirmed by HMC calculations. Neither doublepeaked energy or order parameter
histogram nor volumedependent scaling of the calculated quantities were observed. Fig.1 shows the FSS
calculation of from Eqs. (11) and (12), and from Eqn.
(10) in a lnln plot, for . The best fit value
As the exchange ratio is lowered slightly below , the order of the
transition changes. Still no separate quadrupole ordering is observed, however, the transition becomes
apparently first order. This is clearly seen in Fig.2, where bimodal nature of energy distribution
histograms for various lattice sizes at is obvious. With increasing lattice size the
barrier between the two energy states is increasing, indicating that the energy is discontinuous and
hence, the transition of first order. However, as can be seen from Fig.3, for the case of
the bimodal distribution can only be observed at sufficiently large , indicating
vicinity of the multicritical point. We note that similar behaviour was also observed in both order
parameters distribution histograms, although we do not show it here.
The bimodal distribution vanishes as the ratio drops below the value .
Moreover, below quadrupoles start ordering separately at temperatures higher
than those for dipole ordering. Thus the phase boundary branches and a new middle phase of axial
quadrupole longrange order (QLRO) without magnetic dipole ordering opens between the paramagnetic and
DLRO phases. This phase broadens as decreases, since the QLRO branch is little sensitive
to the ratio variation and levels off, while the DLRO branch turns down approaching the
point (. This means that the ground state is always magnetic as long as
there is a finite dipole exchange interaction. In Fig.4 we present the temperature variation of both
DLRO and QLRO parameters and , respectively, at . We can see that quadrupoles
order before dipoles, forming a fairly broad region of QLRO without DLRO. The question is of what order
are these transitions. MFA predicts a firstorder transition to DLRO phase in this region
chenlevy2 (); chenetal1 (), although some other MFA results allow possibility of a secondorder DLRO
transition in the high biquadratic exchange limit sivar72 (). Unfortunately, no other calculations
by more reliable techniques are available. In fact, determination of the order of a transition here
turned out to be quite tricky. Although no discontinuities in the internal energy temperature variation
were observed, extremely sharp magnetization slopes and consequently, extremely sharp susceptibility
peaks at the transition (Fig.5) would rather suggest a firstorder transition. Moreover, the
fourthorder cumulant at the transition falls to negative values displaying a minimum, i.e.
displays behaviour seemingly typical for a firstorder transition. This, however, must be taken with
precaution since, here, with increasing temperature the system does not enter a paramagnetic phase but a
different ordered phase. Indeed, further investigation showed that the minimum did not scale with the
lattice size as it should in the case of a firstorder transition vollmayr () and, hence, that in
this case it is of different origin and can not be seen as an indicator of a firstorder transition.
Therefore, in order to resolve the ambiguities in the phenomena observed from SMC calculations we
performed a finitesize scaling analysis from HMC data. As mentioned in Introduction, a firstorder
transition can be reliably detected via the fourthorder energy cumulant scaling, which is shown in
Fig.6 for . As we can see, a minimum near extrapolated to the limit
, , reaches the value very close to , as it should be in the
case of a secondorder transition. Also the other quantities clearly did not scale with volume. Instead,
as shown in Fig.7, the obtained scaling indices and match quite well the values expected at a secondorder transition for the given quantities and
universality class (See comments on the universality class in question in the Discussion and concluding
remarks section). Although the transition at appears to be undoubtedly of second
order, it might change to a firstorder one at higher ratio but before the DLRO and QLRO
boundaries merge, as suggested in Ref. sivar72 (). Therefore, we checked for such a possibility but
did not find it to be so. Scaling analysis, similar to that for the case of , performed
for the ratio , which is close to the limiting value ,
produced , a value which deviates slightly more from than in the previous
case, but is still too close to be considered as a sign of a firstorder transition. Also the critical
indices and confirm the secondorder nature
of the transition. Therefore, we conclude that the DLRO transition remains second order in the whole
region up to .
On the other hand, the nature of the separate QLRO transition is much more obvious. All
calculated observables behave typically like those for a secondorder transition and so does the
scaling. This is in agreement with the work of Carmesin carmesin (), who performed mapping of
quadrupolar interactions to dipolar ones for planar systems and, hence, showed that a transition in
these systems should be second order. Critical indices describing a QLRO transition at are extracted from the scaling depicted in Fig.8. Their magnitudes and
are in a good agreement with the values obtained for critical indices
connected with a DLRO transition. These values are also in a good agreement with the values obtained
from the hightemperature series expansion calculations for cubic lattices chenetal2 (); nagata ().
Observing the specific heat temperature dependence we noticed a shoulderlike broadening just
above the quadrupoleordering transition temperature (Fig.9). We ascribe this phenomenon to shortrange
order, which was already predicted to exist by Chen and Levy chenlevy1 (), however, apparently
could not be seen within their MFA scheme. Here, we show a snapshot of the spin arrangement in a XY
plane taken at (inset). In the snapshot we can see that the quadrupoles are not totally
disordered but tend to align locally in regions equally partitioned among the crystalographically
equivalent three directions. At the quadrupoleordering temperature the quadrupoles choose one ordering
direction and the entire crystal forms a single domain.
The resulting phase diagram is drawn in Fig.10 and critical indices calculated for secondorder
transitions are summarized in Table 1.
5.Discussion and conclusions
In this work we studied effects of the biquadratic exchange on the phase diagram of the
classical XY ferromagnet on a stacked triangular lattice (STL). We believe that our investigations,
which are the first of this kind for a system with hexagonal lattice symmetry, covered all significant
phenomena induced by the presence of the biquadratic exchange and present a fairly compact picture of
the role of this higherorder exchange interaction on the critical behaviour of the system considered.
We obtained the phase diagram with two ordered phases, featuring some interesting phenomena such as the
appearance of tricritical and triple points. We found that in the region where the bilinear exchange is
dominant there is only one phase transition to the DLRO phase, which remains second order until the
exchange ratio reaches the value . Upon further lowering of the ratio, the
transition changes to a firstorder one at the tricritical point and remains so down to . Below this value the phase boundary splits into the QLRO transition line at higher
temperatures and the DLRO transition line at lower temperatures, which are both of second order. Hence,
there is a triple point at , where two secondorder and one firstorder
transition boundaries meet.
We speculate that the mechanism driving the system through the mentioned transitions could be as
follows. At () the transition is known to be second order. If
we introduce the biquadratic exchange and reduce , i.e. reduce the bilinear couplings, this
might induce in a sense a kind of tension (competition) between the two exchange interactions. Namely,
while the decreasing bilinear exchange drives the transition temperature down to the lower values, the
biquadratic exchange does not follow this tendency and rather prevents the ordering temperature from
rapid decrease. This tendency is clearly seen from the phase diagram both in the region of separate
transitions, where does not vary much with decreasing , as well as in the region of
simultaneous ordering, where the transition temperature is apparently enhanced by the presence of the
biquadratic exchange (the case of absent biquadratic exchange is represented by the dashdot straight
line in () parameter space). Put differently, quadrupoles would prefer ordering at
higher temperatures but as long as there is a single transition they are prevented to do so by too low
bilinear exchange, and order occurs only if the temperature is lowered still further. This frustration
is enhanced as decreases, and below it results in a firstorder
transition when the strength of the quadrupole ordering prevails and frustrated quadrupoles order
abruptly along with dipoles. However, below the frustration is too high for
the two kinds of ordering to occur simultaneously  they separate, the tension is released and the
transitions become second order again.
The calculated critical indices in the region for QLRO and for DLRO transitions were found to belong to the standard threedimensional XY
universality class, although they showed a slight variation with changing , as had also
been observed in the HTSE calculations chenetal1 (); chenetal2 (). The values of the indices for the
DLRO transition in the region appear to be in a reasonable agreement with
those just mentioned above, and one might tend to include them into the same XY universality class. We
believe, however, that they rather belong to the Ising universality class (Note that the critical
indices of both universality classes take fairly similar values: , and , ferrlandau () for ). The reason
for our claim is that in this region there is first axial quadrupole ordering and only then directional
dipole ordering within the given axis. Therefore, the only difference from the Ising case is that
dipoles can order along any of the three axes, not only the zaxis. The Isinglike nature of the
ordering would also explain the remarkably sharp behaviour of the physical quantities at the transition
as, for example, shown in Fig.5.
We further intend to perform similar simulations on STL with antiferromagnetic bilinear and
antiferromagnetic/ferromagnetic biquadratic exchange. Such a spin system would present a geometrically
frustrated system with additional frustration arising from competition between the antiferromagnetic
bilinear and ferromagnetic biquadratic, or interplane antiferromagnetic bilinear and antiferromagnetic
biquadratic exchange interactions. Among the main points of interest of such investigations will be the
nature of the groundstate order, the universality class of the antiquadrupolar ordering, etc.
(Ref.chenetal2 ())  
0  0.667 0.008  1.333 0.018  1.32 0.03 
0.2  0.660 0.008  1.305 0.011  1.30 0.03 
0.1  0.642 0.006  1.241 0.011   
0.3  0.654 0.009  1.216 0.006   
0.8  0.677 0.008  1.177 0.015   

Footnotes
 The errors for and are calculated from standard errors of the respective slopes in the linear regression .
References
 J. Sivardiere, Phys. Rev. B 6, 4284, (1972).
 J. Sivardiere, A.N. Berker and M. Wortis, Phys. Rev. B 7, 343, (1973).
 H.H. Chen and P. Levy, Phys. Rev. Lett. 27, 1383, (1971); Phys. Rev. B 7, 4267, (1973).
 H.H. Chen and P. Levy, Phys. Rev. B 7, 4284 , (1973).
 R. Micnas, J. Phys. C: Solid St. Phys. 9, 3307, (1976).
 G.S. Chaddha and A. Sharma, J. Magn. Magn. Mater. 191, 373, (1999).
 K.G. Chen, H.H. Chen, C.S. Hsue and F.Y. Wu, Physica 87A, 629, (1977).
 K.G. Chen, H.H. Chen and C.S. Hsue, Physica 93A, 526, (1978).
 H.O. Carmesin, Phys. Lett. A 125, 294, (1987).
 A. Tanaka and T. Idogaki, J. Phys. Soc. Japan 67, 604, (1998).
 M. Campbell and L. Chayes, J. Phys. A 32, 8881, (1999).
 H. Nagata, M.Žukovič and T.Idogaki  unpublished.
 A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61, 2635 , (1988).
 A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63, 1195 , (1989).
 J.C. Le Guillou and J. ZinnJustin, Phys. Rev. Lett. 39, 95 , (1977).
 W. Janke, Phys. Lett. 148, 306 , (1990).
 G.A.T. Allan and D.D. Betts, Proc. Phys. Soc. 91, 341 , (1967).
 K. Vollmayr, J.D. Reger, M. Scheucher and K. Binder, Z. Phys. B 91, 113 , (1993).
 A.M. Ferrenberg and D.P. Landau, Phys. Rev. B 44, 5081, (1991).