References

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 second-order 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 second-order one again as is further increased. Thus, a phase diagram featuring both triple and tricritical points is obtained. Furthermore, a finite-size 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 812-8581, Japan
Tel.: +81-92-642-3811; Fax: +81-92-633-6958
E-mail: milantap@mbox.nc.kyushu-u.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 higher-order interactions in systems with Heisenberg symmetry has been tackled in several mean field approximation (MFA) studies sivar72 ()-chen-levy1 (), by high-temperature series expansion (HTSE) calculations chen-levy2 (), 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 chen-etal1 (); chen-etal2 () 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 second-order 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 first-order 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 long-range order (DLRO), corresponding to the ferromagnetic directional arrangement of spins, and quadrupole long-range 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 bilinear-biquadratic 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 bilinear-biquadratic 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 finite-size scaling (FSS) analysis in order to obtain more reliable data and deliver a complete picture of a long-range 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, short-range 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

 H=−J1∑⟨i,j⟩\boldmathSi⋅\boldmathSj−J2∑⟨i,j⟩(\boldmathSi⋅% \boldmathSj)2 , (1)

where is a two-dimensional 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 heating-cooling loops serve to check possible hysteresis, accompanying first-order 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

 c=(⟨E2⟩−⟨E⟩2)NkBT2 , (2)

the DLRO parameter (magnetization) ,

 m=⟨M⟩/N, where M=⎡⎣(∑iSix)2+(∑iSiy)2⎤⎦12 , (3)

the QLRO parameter ,

 q=⟨Q⟩/N, where Q=⎡⎣(∑i((Six)2−(Siy)2))2+(∑i2SixSiy)2⎤⎦12 , (4)

and the following quantities which are functions of the parameter () :
the susceptibility per site

 χO=(⟨O2⟩−⟨O⟩2)NkBT , (5)

the logarithmic derivatives of and with respect to

 D1O=∂∂Kln⟨O⟩=⟨OE⟩⟨O⟩−⟨E⟩ , (6)
 D2O=∂∂Kln⟨O2⟩=⟨O2E⟩⟨O2⟩−⟨E⟩ , (7)

the fourth-order long-range order (LRO) cumulants and

 U1=1−⟨O4⟩3⟨O2⟩2 ,  U2=2−⟨O4⟩⟨O2⟩2 , (8)

and the fourth-order energy cumulant

 V=1−⟨E4⟩3⟨E2⟩2 . (9)

These quantities are used to estimate the nature of a transition as well as its position. First-order 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 double-checked by the position of the fourth-order LRO and energy cumulants curves intersection for various lattice sizes.
In the next stage we perform HMC calculations, developed by Ferrenberg and Swendsen ferr-swen1 (); ferr-swen2 (), 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 second-order 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 second-order transition, while is expected for a first-order transition ferr-swen1 (); ferr-swen2 (). The extrema of a variety of thermodynamic quantities at a second-order transition are known to scale with a lattice size as, for example:

 χO,max(L)∝LγO/νO , (10)
 D1O,max(L)∝L1/νO , (11)
 D2O,max(L)∝L1/νO , (12)

where and represent the correlation length and susceptibility critical exponents, respectively. In the case of a first-order transition (except for the order parameters), they display a volume-dependent 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 long-range 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 second-order character of a transition in this region was also confirmed by HMC calculations. Neither double-peaked energy or order parameter histogram nor volume-dependent 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 ln-ln plot, for . The best fit value 1 agrees very well with the values obtained from other calculations for , but is apparently lower than the values obtained for ( and from Refs. guillou (); janke ()). However, the critical indices are expected to change by the presence of the biquadratic exchange allan () and, as shown by the series-expansion calculations for chen-etal1 (), is continuously lowered by increasing biquadratic exchange.
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 long-range 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 first-order transition to DLRO phase in this region chen-levy2 (); chen-etal1 (), although some other MFA results allow possibility of a second-order 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 first-order transition. Moreover, the fourth-order cumulant at the transition falls to negative values displaying a minimum, i.e. displays behaviour seemingly typical for a first-order 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 first-order transition vollmayr () and, hence, that in this case it is of different origin and can not be seen as an indicator of a first-order transition. Therefore, in order to resolve the ambiguities in the phenomena observed from SMC calculations we performed a finite-size scaling analysis from HMC data. As mentioned in Introduction, a first-order transition can be reliably detected via the fourth-order 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 second-order 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 second-order 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 first-order 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 first-order transition. Also the critical indices and confirm the second-order 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 second-order 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 high-temperature series expansion calculations for cubic lattices chen-etal2 (); nagata ().
Observing the specific heat temperature dependence we noticed a shoulder-like broadening just above the quadrupole-ordering transition temperature (Fig.9). We ascribe this phenomenon to short-range order, which was already predicted to exist by Chen and Levy chen-levy1 (), 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 quadrupole-ordering 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 second-order 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 higher-order 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 first-order 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 second-order and one first-order 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 dash-dot 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 first-order 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 three-dimensional XY universality class, although they showed a slight variation with changing , as had also been observed in the HTSE calculations chen-etal1 (); chen-etal2 (). 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 , ferr-landau () 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 z-axis. The Ising-like 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 inter-plane antiferromagnetic bilinear and antiferromagnetic biquadratic exchange interactions. Among the main points of interest of such investigations will be the nature of the ground-state order, the universality class of the antiquadrupolar ordering, etc.

### Footnotes

1. The errors for and are calculated from standard errors of the respective slopes in the linear regression .

### References

1. J. Sivardiere, Phys. Rev. B 6, 4284, (1972).
2. J. Sivardiere, A.N. Berker and M. Wortis, Phys. Rev. B 7, 343, (1973).
3. H.H. Chen and P. Levy, Phys. Rev. Lett. 27, 1383, (1971); Phys. Rev. B 7, 4267, (1973).
4. H.H. Chen and P. Levy, Phys. Rev. B 7, 4284 , (1973).
5. R. Micnas, J. Phys. C: Solid St. Phys. 9, 3307, (1976).
6. G.S. Chaddha and A. Sharma, J. Magn. Magn. Mater. 191, 373, (1999).
7. K.G. Chen, H.H. Chen, C.S. Hsue and F.Y. Wu, Physica 87A, 629, (1977).
8. K.G. Chen, H.H. Chen and C.S. Hsue, Physica 93A, 526, (1978).
9. H.-O. Carmesin, Phys. Lett. A 125, 294, (1987).
10. A. Tanaka and T. Idogaki, J. Phys. Soc. Japan 67, 604, (1998).
11. M. Campbell and L. Chayes, J. Phys. A 32, 8881, (1999).
12. H. Nagata, M.Žukovič and T.Idogaki - unpublished.
13. A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61, 2635 , (1988).
14. A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63, 1195 , (1989).
15. J.C. Le Guillou and J. Zinn-Justin, Phys. Rev. Lett. 39, 95 , (1977).
16. W. Janke, Phys. Lett. 148, 306 , (1990).
17. G.A.T. Allan and D.D. Betts, Proc. Phys. Soc. 91, 341 , (1967).
18. K. Vollmayr, J.D. Reger, M. Scheucher and K. Binder, Z. Phys. B 91, 113 , (1993).
19. A.M. Ferrenberg and D.P. Landau, Phys. Rev. B 44, 5081, (1991).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters