Phase diagram of the repulsive Blume–Emery–Griffiths model in the presence of an external magnetic field on a complete graph

The Blume–Emery–Griffiths (BEG) model is the simplest model that incorporates biquadratic interactions [1]. The presence of a biquadratic exchange interaction is known to be relevant to understanding the properties of rare-earth compounds. Biquadratic exchange was first suggested by Kittel as part of the theory of magnetoelastic effects in NiAs-type structures [2], and by Anderson in the superexchange interaction of iron group oxides and fluorides [3]. In rare-earth compounds, the unpaired 4f electrons lie deep inside the 5d and 5s orbitals, so these electrons do not experience the strong crystal field generated by other ions in the crystal. Therefore, their spherically symmetric potential is not completely destroyed and, as a result, their orbital angular momentum is not entirely quenched. The superexchange between these unquenched orbital momenta gives rise to a biquadratic exchange interaction term in the Hamiltonian [4]. Other interactions, such as phonon exchange between ions [5] and Schrödinger's spin-one exchange operator [6], can also result in the inclusion of such an interaction. Both attractive and repulsive biquadratic interactions are of interest. The requirement for a small repulsive exchange interaction in a Hamiltonian was first mentioned by Harris and Owen [7] and Rodbell et al [8] in order to explain the paramagnetic resonance of the Mn ion pairs that are present as an impurity in the crystals of MgO.

Biquadratic exchange interaction is represented by a fourth-order term in spin operators. The spin-one BEG model has been shown to successfully capture the physics of these higher-order interactions and has been widely studied. It incorporates a uniform crystal field (Δ) and a biquadratic exchange interaction (K) along with the bilinear exchange interaction term. This model was first introduced in order to explain the phase separation and superfluidity of a ³He − ⁴He mixture [1]. Apart from this example, many other physical systems, such as metamagnets, liquid crystals, semiconducting alloys, microemulsions, etc. can also be mapped onto the BEG model. This model has a rich phase diagram that depends on the sign and magnitude of the biquadratic term. The special case, K = 0 is known as the Blume–Capel model. The Blume–Capel model was first studied by Blume [9] and Capel [10] in order to explain the first-order transition in UO₂. Another extreme case of zero bilinear exchange was studied by Griffiths [11].

The attractive biquadratic exchange interaction (K > 0) BEG model has been extensively studied. Its phase diagram changes with the value of K. For small K, there is a transition from a ferromagnetic to a paramagnetic phase in the (T − Δ) plane. This transition line changes from a continuous to a first-order transition line at a tricritical point (TCP). As K increases further, another paramagnetic state emerges and the two paramagnetic states are separated by another first-order line. The two first-order lines meet at a triple point. For larger values of K, the continuous transition line terminates on the first-order line at a critical end point (CEP), and the TCP disappears. The phase diagram has been extensively studied using various techniques such as the mean-field theory [1, 12–15], cluster variation [16], the Bethe lattice [17], high-temperature series expansion [18], etc. Apart from the mean field, the phase diagram has also been studied in finite dimensions using renormalization group [19, 20] and Monte-Carlo simulations in two and three dimensions [21–23]. However, the simulations have mostly been performed using the Blume–Capel model to study the continuous transition line and the TCP. The other multicritical points have not been studied, as they are hard to locate in simulations.

For a repulsive biquadratic interaction (K < 0), competition between the biquadratic and the bilinear interactions gives rise to a very different behavior from the behavior for positive K. The negative K term prefers the non-magnetic spins to the magnetic ones. This creates competition between the magnetic and non-magnetic spins. In a recent study on a complete graph [24] in the (T − Δ) plane, it was shown that as K becomes more negative, the TCP changes to a quadrupolar point at K = −0.0828 and K = −0.1838 in the microcanonical and canonical ensembles, respectively. The authors of that paper studied the system for small negative values of K up to K = −0.4. In this paper, we study the BEG model using a complete graph for the entire range of K, with an emphasis on the large negative K regime in the (T − Δ − H) space. Since the multicritical points occur in systems described by three or more thermodynamic fields, it is useful to study them in (T − Δ − H) space. Compared to the attractive BEG model, the repulsive BEG has been less well studied. Though there have been some studies in the canonical ensemble on bipartite lattices using the mean-field method [24, 25], renormalization group [26–28], Monte Carlo simulations in two and three dimensions [29–33], cluster-variation methods [34], and hierarchical models [35].

We solve the phase diagram in both the microcanonical and canonical ensembles. We observe four topologies of the phase diagram, depending on the values of K. We find that for low negative values of K there is a TCP at which three critical lines (λ, λ_±) meet (see figure 1(a)). The λ line is the line of continuous transition between the ferromagnetic phase (with magnetization m ≠ 0) and the paramagnetic phase (with magnetization m = 0) in the H = 0 plane. The λ_± lines are the lines of continuous transition in the ±H planes, respectively. The λ_± lines enclose two first-order surfaces which meet along a triple line in the H = 0 plane. Since these first-order surfaces appear symmetrically in the phase diagram like the wings of a bird, they are referred to as 'wings' [1]. As we decrease K, TCP becomes a quadrupolar point at K = −0.0828 and K = −0.1838 in the microcanonical and canonical ensembles, respectively. On reducing K further, the λ_± lines move inside the ordered region and meet at a new multicritical point, the bicritical end point (BEP), and the λ line is truncated at the first-order line at a CEP, as shown in figure 1(b). Previously, this BEP has been reported to be an ordered critical point (CP) [24–34]. By introducing the external field, we realize that this point is the junction of two critical lines (λ_±) and hence a BEP. We find that the width of the wings in temperature shrink as K approaches K = −1. At exactly K = −1, the BEP and the CEP move to T = Δ = 0 and the wings vanish (see figure 1(c)). With a further reduction in K, we see no transition in the finite H plane. There is only a transition from a ferromagnetic state to a paramagnetic state in the H = 0 plane (figure 1(d)). The area under the λ line shrinks as K becomes more and more negative. At K → −∞, only the paramagnetic state survives and there is no transition.

Zoom In Zoom Out Reset image size

The competition introduced by the repulsive biquadratic interaction makes this a very interesting model to study. In fact, we find that the phase diagram for the repulsive (−1 < K ⩽ 0) BEG model is similar to the topology of the phase diagram for the Blume–Capel model with a random crystal field recently studied [36] for intermediate and weak disorder. In this paper, we examine ensemble inequivalence by not only looking at the first-order line in the (T − Δ) plane, but also by computing the critical lines (λ_±) in the H ≠ 0 plane. We find that these two critical lines and also the multicritical points are different in the two ensembles in general. Another interesting observation we have is that for K ⩽ −1, the two ensembles are equivalent. The attractive BEG model has been comprehensively studied, and we find results similar to those reported in earlier studies in the (T − Δ) plane [1, 12–21]. However, in the (T − Δ − H) space, we found non-monotonic behavior of the wings in terms of temperature as K becomes greater than K = 1. This, as far as we know, has not been previously reported.

The structure of this paper is as follows: in section 2, we introduce the BEG model and discuss its zero-temperature phase diagram. In sections 3 and 4, we derive the equations of the critical lines in the (T − Δ − H) space for both repulsive and attractive BEG models for the canonical ensemble and the microcanonical ensemble, respectively. In section 5 we discuss the ensemble inequivalence in detail. We conclude in section 6.

The Hamiltonian of the BEG model on a complete graph in the presence of an external magnetic field is given by:

$$\begin{equation}\mathcal{H}=-\frac{1}{2N}{\left(\sum\limits _{i}{S}_{i}\right)}^{2}-\frac{K}{2N}{\left(\sum\limits _{i}{S}_{i}^{2}\right)}^{2}+{\triangle}\sum\limits _{i}{S}_{i}^{2}-H\sum\limits _{i}{S}_{i}\end{equation} \tag{ 1 }$$

where S_i can take three values ±1 and 0, H is a constant external field coupled with the order parameter, △ is the crystal field, and K is the biquadratic interaction coefficient. The two order parameters are magnetization, ${x}_{1}={\sum }_{i}\frac{{S}_{i}}{N}$ and the density of the ±1 spins, ${x}_{2}={\sum }_{i}\frac{{S}_{i}^{2}}{N}$. For any finite K, as Δ → −∞, this model becomes equivalent to the Ising model, as the spins take only ±1 values. As Δ increases, the number of vacancies in the system increases. For negative K, spins are more likely to take the value of zero. At finite temperatures, when K < 0, both the biquadratic term and the crystal field term prefer zero spins. Hence, the λ transition occurs at lower values of Δ as K decreases. On the other hand, for positive K the magnetic spins are more likely to be chosen. Hence, when K > 0 the biquadratic and crystal field terms in the Hamiltonian compete, and the λ transition occurs at a higher Δ for positive K.

First, let us examine the zero-temperature phase diagram of the system. The energy per particle can be written as (from the Hamiltonian equation (1)): ${\epsilon}=-\frac{1}{2}\left({x}_{1}^{2}+K{x}_{2}^{2}\right)+{\Delta}{x}_{2}-H{x}_{1}$. When all the spins are zero, the energy is = 0. Apart from this paramagnetic phase, there are other states which are possible, depending on the parameter values. For −1 ⩽ K ⩽ ∞, the ferromagnetic state, x₁ = ±1 and x₂ = 1 dominates. The energy of this state is ${\epsilon}=-\frac{1}{2}\left(1+K\right)+{\Delta}$. If 2Δ > 1 + K, then the phase is paramagnetic, and for 2Δ < 1 + K, the phase is ferromagnetic. At exactly 2Δ = 1 + K there is a first-order phase transition. For K < −1, the term $-\frac{1}{2}\left(1+K\right)$ in the energy contributes a positive value. Therefore, for any Δ ⩾ 0, the paramagnetic phase is the stable state. As Δ becomes negative, there is another ferromagnetic state where |x₁| = x₂ < 1, which becomes stable when $\vert {\Delta}\vert {< }-\frac{1}{2}\left(1+K\right)$. For $\vert {\Delta}\vert { >}-\frac{1}{2}\left(1+K\right)$, the state with |x₁| = x₂ = 1 becomes stable and there is now a first-order transition between these two ferromagnetic states at ${\Delta}=\frac{1+K}{2}$.

Given the Hamiltonian (equation (1)), the probability of the spin configuration (C_N = {S_i }) for N spins can be expressed as:

$$\begin{equation}P\left({C}_{N}\right)=\frac{{\text{e}}^{\beta \left(\frac{1}{2N}{\left({\sum }_{i}{S}_{i}\right)}^{2}+\frac{K}{2N}{\left({\sum }_{i}{S}_{i}^{2}\right)}^{2}-{\triangle}{\sum }_{i}{S}_{i}^{2}+H{\sum }_{i}{S}_{i}\right)}}{{Z}_{N}}.\end{equation} \tag{ 2 }$$

Here, Z_N is the partition function and $\beta =\frac{1}{T}$. The free energy of the mean-field models can be calculated in many ways. We calculated the free energy of the system using the large deviation principle (LDP) [37], which states that the probability of a spin configuration with a magnetization of x₁ and a density of x₂ for N → ∞ can be expressed as:

$$\begin{equation}P\left({C}_{N}\right)\asymp {\mathrm{e}}^{-NI\left({x}_{1},{x}_{2}\right)},\end{equation} \tag{ 3 }$$

where I(x₁, x₂) is called the rate function. For a system to satisfy LDP, the following limit should hold:

$$\begin{equation}I\left({x}_{1},{x}_{2}\right)\asymp -\underset{N\to \infty }{\mathrm{lim}}\enspace \frac{1}{N}\enspace \mathrm{ln}P\left({C}_{N}\right).\end{equation} \tag{ 4 }$$

The rate function I(x₁, x₂) can be seen as the full Landau free-energy functional for the system. Minimizing the rate function w.r.t. x₁ and x₂ gives the free energy of the system. The detailed calculation of the rate function I(x₁, x₂) is shown in appendix A. The minimization of the rate function w.r.t. x₁ and x₂ gives the two following coupled equations for the two order parameters:

$$\begin{equation}m=\frac{2\enspace {\text{e}}^{\beta \left(Kq-{\Delta}\right)}\enspace \mathrm{sinh}\enspace \beta \left(m+H\right)}{1+2\enspace {\text{e}}^{\beta \left(Kq-{\Delta}\right)}\enspace \mathrm{cosh}\enspace \beta \left(m+H\right)}\end{equation} \tag{ 5 }$$

$$\begin{equation}q=\frac{2\enspace {\text{e}}^{\beta \left(Kq-{\Delta}\right)}\enspace \mathrm{cosh}\enspace \beta \left(m+H\right)}{1+2\enspace {\text{e}}^{\beta \left(Kq-{\Delta}\right)}\enspace \mathrm{cosh}\enspace \beta \left(m+H\right)},\end{equation} \tag{ 6 }$$

where m and q are the extremums of x₁ and x₂. For m ≠ 0, the two fixed-point equations are connected via:

$$\begin{equation}q=m\enspace \mathrm{coth}\enspace \beta \left(m+H\right)\end{equation} \tag{ 7 }$$

and the free energy at the fixed point can be written as (substituting equation (7) into equation (A10));

$$\begin{align}\hfill f\left(m\right)& =\frac{\beta {m}^{2}}{2}+\frac{\beta K{m}^{2}\enspace {\mathrm{coth}}^{2}\enspace \beta m}{2}+\mathrm{log}\left(1+2\enspace {\text{e}}^{-\beta {\Delta}}\enspace \mathrm{cosh}\enspace \beta H\right)\hfill \\ \hfill & \quad -\mathrm{log}\left(1+2\enspace {\text{e}}^{\beta \left(Km\enspace \mathrm{coth}\enspace \beta m-{\Delta}\right)}\enspace \mathrm{cosh}\enspace \beta \left(m+H\right)\right).\hfill \end{align} \tag{ 8 }$$

For H = 0, the system has a line of continuous transition (λ line) in the (T − Δ) plane. The equation of this line can be obtained by linearizing equation (7). On linearizing, we get the equation of the λ line, which is

$$\begin{equation}2\left(\beta -1\right)={\text{e}}^{\beta {\Delta}-K}.\end{equation} \tag{ 9 }$$

For H ≠ 0, the system has a line of continuous transition in the H₊ and H₋ planes, separating the two magnetic states. To calculate these lines, we take f'(m) = f''(m) = f‴(m) = 0 and f‴(m) > 0, to get the loci of the critical points. This gives the following equations:

$$\begin{equation}{f}_{1}\equiv m-\frac{az}{C}=0\end{equation} \tag{ 10 }$$

$$\begin{equation}{f}_{2}\equiv \frac{1}{za\beta }-\frac{az{y}^{2}+\sqrt{{y}^{2}+1}y+yK\left(\sqrt{{y}^{2}+1}-\beta m{y}^{2}\right)}{{C}^{2}}=0\end{equation} \tag{ 11 }$$

$$\begin{align}\hfill {f}_{3}& \equiv \frac{{\beta }^{3}z}{{C}^{3}}\left(az{K}^{2}\sqrt[\frac{3}{2}]{{y}^{2}+1}-K\left({y}^{2}+1\right)y\left[\left(2+K\right)+2az\left(1+K\right)my\beta \right]\right.\hfill \\ \hfill & \quad +\sqrt{{y}^{2}+1}\left[az\left(1+2K\right)+2\beta {K}^{2}m{y}^{3}+{\beta }^{2}az{K}^{2}{m}^{2}{y}^{4}\right]\hfill \\ \hfill & \quad -\left.y\left[1-2{a}^{2}{z}^{2}+2\beta Kazym-2K{y}^{2}-2\beta azKm{y}^{3}+{\beta }^{2}{K}^{2}{m}^{2}{y}^{4}\right]\right)=0\hfill \end{align} \tag{ 12 }$$

where a = e^−βΔ, y = cosech β(H + m), z = 2 e^{βKm coth β(H+m)}, and $C=\sqrt{{y}^{2}+1}za+y$. For K = 0, the above equations reduce to the following equations for the pure Blume–Capel model:

$$\begin{equation}m={\pm}\frac{2a\sqrt{{x}^{2}-1}}{\left(2ax+1\right)}\end{equation} \tag{ 13 }$$

$$\begin{equation}\frac{2a+x}{{\left(2ax+1\right)}^{2}}=\frac{1}{2a\beta }\end{equation} \tag{ 14 }$$

$$\begin{equation}\frac{8{a}^{2}+2ax-1}{{\left(2ax+1\right)}^{3}}=0\end{equation} \tag{ 15 }$$

where x = cosh β(m + H). The solutions of equations (14) and (15) are hence given by:

$$\begin{equation}x=\mathrm{cosh}\enspace \beta \left(m+H\right)=\frac{\beta -2}{\sqrt{4-\beta }}\end{equation} \tag{ 16 }$$

$$\begin{equation}a={\mathrm{e}}^{-\beta {\Delta}}=\frac{\sqrt{4-\beta }}{4}\end{equation} \tag{ 17 }$$

and the critical lines for the H ≠ 0 plane are the following:

$$\begin{equation}m={\pm}\sqrt{\frac{\beta -3}{\beta }}\end{equation} \tag{ 18 }$$

$$\begin{equation}H={\pm}\frac{1}{\beta }\enspace \mathrm{log}\left(\frac{\beta -2+\sqrt{{\beta }^{2}-3\beta }}{\sqrt{4-\beta }}\right)-m.\end{equation} \tag{ 19 }$$

These lines are the λ_± lines [1] (depending on the sign of H). These lines enclose two first-order surfaces in the H ≠ 0 plane, called the wings.

For K ≠ 0, solving equations (10)–(12) is not analytically possible. Hence, we use graphical methods to get the coordinates of the critical points in the (T − Δ − H) space for a given K. We plot f₁, f₂, and f₃ in the (β − Δ) plane for m, fixing K and varying the value of H. The value of m for a fixed value of H and K at which the three equations meet gives the co-ordinates of the critical point. If we now take H = 0 in equations (10)–(12), we get the coordinates of the point of intersection of the λ_± lines in the (T − Δ) plane. We can now use these to locate the multicritical points (TCP and BEP) in the (T − Δ) plane. We use these to obtain the phase diagram for various values of K. For example: figure 2(a) is a contour plot of f₁, f₂, and f₃ in the (β − Δ) plane at K = −0.4 and H = 0. The intersection of the three functions gives the co-ordinates of the critical point. We find that for K > −0.1838 the functions only intersect for m = 0, which is the point where the λ_± lines meet the λ line. Hence, this point is the TCP. For the range −0.1838 > K > −1, we find that the intersection occurs for m ≠ 0. This m ≠ 0 solution gives the locus of the BEP where the λ_± lines meet the H = 0 plane in this regime of K. Interestingly, we find that for K < −1 the three functions never intersect at the same point for any m. For example, in figure 2(b), we plot three functions for K = −2. We will discuss these results in detail in the next section (section 3.1).

Zoom In Zoom Out Reset image size

3.1. The repulsive Blume–Emery–Griffiths model

In this section, we analyze the results of the repulsive BEG model. We find that for 0 ⩾ K ⩾ −0.1838, two critical lines (λ_±) at H ≠ 0 meet the λ line (at H = 0) at the TCP. The temperature decreases exponentially and Δ increases linearly with increasing H along the λ_± lines, as shown in figure 3. As K becomes more negative (for −0.1838 > K > −1), the (λ_±) lines no longer meet the λ line, instead they enter the ordered region and meet at the first-order surface (H = 0) at a BEP. In this case, the λ line terminates at the first-order line at a CEP.

Zoom In Zoom Out Reset image size

As K approaches −1, the wing width in temperature (which is the difference between the temperature at the BEP (T_BEP) and the saturation value of the temperature (T_sat) at which both the Δ and H → ∞) denoted by T_w starts to shrink. At exactly K = −1 the BEP, CEP, and T_w reach zero. As we decrease K further, we find that there are no transitions in the H plane. This is also supported by the fact that now there are no multicritical points in the (T − Δ) plane. Hence, we conclude that for K < −1, the wing surfaces completely disappear. The phase diagram just consists of a continuous transition line (λ line) from the ferromagnetic phase to the paramagnetic phase in the H = 0 plane. For large negative K, the area enclosed by the λ line in the (T − Δ) plane shrinks. At K → −∞, there is no phase transition, and only the S = 0 state dominates. The decreasing width of the wings with decreasing K is shown in figure 4. We also observe that T_sat can be approximated numerically as T_sat ≃ (K + 1)/4. This will be discussed in more detail in section 4.1, where we obtain similar results in the microcanonical ensemble. The values of T_w, T_sat, and the co-ordinates of the multicritical points (TCP, BEP) for the repulsive BEG model are listed in table 1.

Zoom In Zoom Out Reset image size

Table 1. Co-ordinates of the TCP and BEP for different values of K. T_sat is the saturation value of the temperature at which both the Δ and H → ∞. T_w is the width of the wing lines for different K.

Canonical: −1 ⩽ K ⩽ 0
	TCP/BEP		Tsat	Tw
K	T	Δ	≃(K + 1)/4	(=TTCP/BEP − Tsat)
0	0.333 33	0.462 098	0.25	0.083 3333
−0.05	0.310 3448	0.447 41	0.237 501	0.072 843 316
−0.1	0.285 7142	0.431 268	0.225 0124	0.060 7018
−0.2	0.231 2737	0.391 831	0.2	0.031 273 76
−0.3	0.187 5335	0.346 377	0.174 999 56	0.012 533 94
−0.4	0.154 462 22	0.298 727	0.149 999 25	0.004 462 97
−0.6	0.100 192 068	0.199 958	0.1	0.000 192 068
−0.7	0.075 009 189	0.149 998	0.075 000 18	0.000 009 009
−0.8	≃0.050 000 25	≃0.1	≃0.049 999 875	≃0.000 000 375
−0.9	≃0.024 999 969	≃0.05	≃0.024 999 968	≃0.000 000 001
−0.999	≃0.000 25	≃0.0005	≃0.000 25	≃0.00

The absence of a phase transition for K < −1 in the H plane can also be seen by looking at the magnetization and susceptibility. We find that the magnetic susceptibility diverges around the expected critical point for K > −1. On the other hand, for K ⩽ −1, magnetic susceptibility is finite over the entire H plane. We plotted the magnetization and the susceptibility for H = 0.5 at K = −0.6 and K = −1.2. In figure 5(a) we show them as a function of Δ for K = −0.6, by fixing T = 0.1. The susceptibility shows singular behavior at Δ = 0.7. The point of divergence matches the co-ordinates of the transition obtained from equations (10)–(12). On the other hand, for K = −1.2, we find no such divergence. In figure 5(b) for K = −1.2, by fixing T = 0.025, we plot the magnetization and find that it changes continuously along Δ and the susceptibility shows a cusp but does not diverge. Though we only plotted the results for a fixed T, we have checked the entire plane by changing the value of T. The magnetic susceptibility has no divergence for any T.

Zoom In Zoom Out Reset image size

3.2. The attractive Blume–Emery–Griffiths model

The attractive BEG model has already been extensively studied by various authors [1–14] and the topology of the phase diagram is known as a function of K in the (T − Δ) plane. We observe a similar phase diagram topology. We study the (T − Δ − H) phase diagram and find that the topology of the phase diagram for different K's are similar to [14]. To recap, we find that for 0 < K ⩽ 2.78, the phase diagram is similar to what we find for 0 ⩾ K ⩾ −0.1838. The λ_± lines meet at the TCP. For 2.78 < K < 3, a new first-order surface appears separating two paramagnetic states: P2 (m = 0, q₋ < 0.5) and P1 (m = 0, q₊ > 0.5). This surface meets the first-order line (at H = 0) at a triple point. This new first-order surface terminates on a line of critical points (at H ≠ 0 plane). As K changes from K = 2.78, this line of critical points in the paramagnetic region moves to higher temperatures, and at exactly K = 3 it intersects the λ_± lines and then extends to infinity. For 3 < K ⩽ 3.8, the λ_± lines terminate at the first-order surface that separates the P1 and P2 phase, and become finite. For K > 3.8, the λ line terminates at a CEP, and thus the wings vanish.

We observe that for K > 1, the wings show non-monotonic behavior in temperature, in contrast to what happens in the range of −1 < K < 1 (figure 3). For small values of H, the λ_± lines move towards lower values of T and higher values of Δ. As H becomes larger, these lines start moving towards higher Δ and higher T as shown in figure 6. This non-monotonic behavior observed in the wings for all K ⩾ 1 values can be interpreted as follows: for any positive K, there are two possible solutions for q for a fixed value of energy, q_± (for more detail, see section 4) with q₊ > q₋. When K ⩾ 1, the first term in the energy (mentioned in section 2) dominates over the crystal field term and the density of the S = ±1 spin increases. For smaller values of H the q₋ solution dominates in the system and the wings show monotonic behavior, as before. When H increases further, the q₊ solution becomes favorable, which in turn lowers the energy. Thus the energy–entropy balance occurs at a higher temperature.

Zoom In Zoom Out Reset image size

In order to analyze the system in the microcanonical ensemble, we need to express the energy in terms of the number of particles with spin ±1 and 0. Let us assume the number of particles with ±1 spin are N_± and the number of particles with zero spin are N₀, such that N = N₊ + N₋ + N₀, where N is the total number of particles in the system. The energy of the system can thus be written as,

$$\begin{equation}E={\Delta}Q-\frac{1}{2N}{M}^{2}-\frac{K}{2N}{Q}^{2}-HM\end{equation} \tag{ 20 }$$

where M = N₊ − N₋ is the total magnetization and Q = N₊ + N₋ is the spin density of the system. In terms of m (=M/N) and q (=Q/N), the expression for the energy is,

$$\begin{equation}{\epsilon}={\Delta}q-\frac{1}{2}{m}^{2}-\frac{K}{2}{q}^{2}-Hm\end{equation} \tag{ 21 }$$

where ${\epsilon}=\frac{E}{N}$ is the energy per particle, and m and q are the single site magnetization and density (as mentioned in section 2). The total number of microstates of the system can be written in terms of N, N₊, N₋ and N₀ as,

$$\begin{equation}{\Omega}=\frac{N!}{{N}_{+}!{N}_{-}!{N}_{0}!}.\end{equation} \tag{ 22 }$$

In the limit when N₊, N₋, N₀ are large, the expression for entropy, i.e. S = k_B ln(Ω) can be written using the Stirling approximation as,

$$\begin{equation}s=\frac{S}{{k}_{\mathrm{B}}N}=q\enspace \mathrm{ln}\left(2\right)-\left(1-q\right)\enspace \mathrm{ln}\left(1-q\right)-\frac{1}{2}\left(q+m\right)\enspace \mathrm{ln}\left(q+m\right)-\frac{1}{2}\left(q-m\right)\enspace \mathrm{ln}\left(q-m\right)\end{equation} \tag{ 23 }$$

where s is the entropy per particle of the system. The equilibrium entropy can be obtained by maximizing the entropy of equation (23) with respect to m and q. We can express q in terms of m and the other variables as,

$$\begin{equation}{q}_{{\pm}}=\frac{{\Delta}}{K}{\pm}{\gamma }^{1/2}\end{equation} \tag{ 24 }$$

where $\gamma ={\left(\frac{{\Delta}}{K}\right)}^{2}-\frac{2{\epsilon}}{K}-\frac{{m}^{2}}{K}-\frac{2Hm}{K}$. For K = 0, the expression has a much simpler form, $q=\frac{1}{{\Delta}}\left({\epsilon}+\frac{1}{2}{m}^{2}+Hm\right)$.

Since there are two values of q, the one that is in the range [0, 1] is accepted. There is also a possibility that both q values are in the range [0, 1]; in this case, the equilibrium entropy will be the one with the maximum value at its corresponding equilibrium m. We find that for K < 0, only q₋ is acceptable, however, for K > 0, both the q_± solutions are acceptable.

We now aim to find the second-order transition line in the (T − Δ − H) space. In the H = 0 plane, the value of the magnetization m on the line of continuous transition is zero, however, for any nonzero H, the magnetization m has a nonzero value on the continuous transition line. In order to obtain this continuous transition line, we need to equate the first three derivatives of s (with respect to m) to zero, with the constraint that the fourth derivative is negative. The first four derivatives of the entropy s are:

$$\begin{equation}\frac{\partial s}{\partial m}={q}^{\prime }\enspace \mathrm{ln}\left\{\frac{2\left(1-q\right)}{\sqrt{{q}^{2}-{m}^{2}}}\right\}-\mathrm{ln}\sqrt{\frac{q+m}{q-m}}\end{equation} \tag{ 25 }$$

$$\begin{equation}\frac{{\partial }^{2}s}{\partial {m}^{2}}={q}^{{\prime\prime}}\enspace \mathrm{ln}\left\{\frac{2\left(1-q\right)}{\sqrt{{q}^{2}-{m}^{2}}}\right\}-\frac{{q}^{\prime 2}}{1-q}-\frac{1}{2}\left\{\frac{{\left({q}^{\prime }+1\right)}^{2}}{q+m}+\frac{{\left({q}^{\prime }-1\right)}^{2}}{q-m}\right\}\end{equation} \tag{ 26 }$$

$$\begin{align}\hfill \frac{{\partial }^{3}s}{\partial {m}^{3}}& ={q}^{\prime \prime \prime }\enspace \mathrm{ln}\left\{\frac{2\left(1-q\right)}{\sqrt{{q}^{2}-{m}^{2}}}\right\}-\frac{{q}^{\prime 3}}{{\left(1-q\right)}^{2}}+\frac{1}{2}\left\{\frac{{\left({q}^{\prime }+1\right)}^{3}}{{\left(q+m\right)}^{2}}+\frac{{\left({q}^{\prime }-1\right)}^{3}}{{\left(q-m\right)}^{2}}\right\}\hfill \\ \hfill & \quad -\frac{3}{2}\left\{\frac{2{q}^{\prime }{q}^{{\prime\prime}}}{1-q}+\frac{{q}^{{\prime\prime}}\left({q}^{\prime }+1\right)}{q+m}+\frac{{q}^{{\prime\prime}}\left({q}^{\prime }-1\right)}{q-m}\right\}\hfill \end{align} \tag{ 27 }$$

$$\begin{align}\hfill \frac{{\partial }^{4}s}{\partial {m}^{4}}& ={q}^{\prime \prime \prime \prime }\enspace \mathrm{ln}\left\{\frac{2\left(1-q\right)}{\sqrt{{q}^{2}-{m}^{2}}}\right\}-2{q}^{\prime \prime \prime }\left\{\frac{{q}^{\prime }+1}{q+m}+\frac{{q}^{\prime }-1}{q-m}+\frac{2{q}^{\prime }}{1-q}\right\}\hfill \\ \hfill & \quad +3{q}^{{\prime\prime}}\left\{\frac{{\left({q}^{\prime }+1\right)}^{2}}{{\left(q+m\right)}^{2}}+\frac{{\left({q}^{\prime }-1\right)}^{2}}{{\left(q+m\right)}^{2}}-\frac{2{q}^{\prime 2}}{{\left(1-q\right)}^{2}}\right\}\hfill \\ \hfill & \quad -\frac{3}{2}{q}^{\prime \prime 2}\left\{\frac{2}{1-q}+\frac{1}{q+m}+\frac{1}{q-m}\right\}-\frac{{\left({q}^{\prime }+1\right)}^{4}}{{\left(q+m\right)}^{3}}-\frac{{\left({q}^{\prime }-1\right)}^{4}}{{\left(q-m\right)}^{3}}-\frac{2{q}^{\prime 4}}{{\left(1-q\right)}^{3}}\hfill \end{align} \tag{ 28 }$$

where q', q'' ...... are partial derivatives of q w.r.t. m. We solve the above first three equations numerically and obtain a set of physical solutions (Δ, , m), such that the fourth derivative is negative. We then calculate the temperature, using the relation $\beta =\frac{\partial s}{\partial {\epsilon}}$,

$$\begin{equation}\beta =\mp \frac{1}{K{\gamma }^{1/2}}\enspace \mathrm{ln}\left\{\frac{2\left(1-{q}_{{\pm}}\right)}{\sqrt{{q}_{{\pm}}^{2}-{m}^{2}}}\right\}\end{equation} \tag{ 29 }$$

which gives the equivalent phase diagram in the (T − Δ − H) space.

4.1. The repulsive Blume–Emery–Griffiths model

In this section, we show our results for the repulsive BEG model in the microcanonical ensemble in the (T − Δ − H) space. In the absence of a magnetic field, this model has been recently studied in [12, 24]. We find that for −0.0828 ⩽ K ⩽ 0, the phase diagram consists of a TCP where the λ_± lines meet the λ line in the H = 0 plane. As the K decreases further, for −1 < K < −0.0828, it was previously reported in [24] that a CP appears in the ordered region of the system along with a CEP. In this topology, as we switch on the field H, we note that the λ_± lines meet at the proposed CP. Thus, the CP is actually a BEP. We show our results for K = −0.4 in the (Δ, H), (, H) and (T, H) planes in figure 7. Here, we show the behavior of the λ₊ line for positive H. We note that the value of Δ on the λ₊ line increases almost linearly with H in the large H limit. The values of and T decrease with H and saturate for large H. We note that the variation of Δ in the large H limit is of the type: Δ ≃ (K + 1)/2 + H. Also, the saturation values are _sat ≃ (K + 1)/8 and T_sat ≃ (K + 1)/4. The values of (Δ, , T) for the BEP and the saturation values of and T are listed in table 2.

Zoom In Zoom Out Reset image size

Table 2. Coordinates of the multicritical points (TCP, BEP), saturation values of and T and the width of the wings for −1 < K ⩽ 0.

Microcanonical: −1 ⩽ K ⩽ 0
	TCP/BEP			sat	Tsat	w	Tw
K	Δ		T	≃(K + 1)/8	≃(K + 1)/4	(=TCP/BEP − sat)	(=TTCP/BEP − Tsat)
0	0.462 40	0.152 75	0.330 33	0.125 02	0.250 07	0.027 73	0.080 26
−0.05	0.447 41	0.141 25	0.310 32	0.118 75	0.237 50	0.0225	0.072 82
−0.1	0.430 79	0.125 56	0.279 64	0.112 50	0.225 00	0.013 06	0.054 64
−0.2	0.391 00	0.103 43	0.224 54	0.100 00	0.200 00	0.003 43	0.024 54
−0.3	0.346 24	0.088 37	0.185 64	0.087 50	0.175 00	0.000 87	0.010 64
−0.4	0.298 71	0.075 19	0.154 01	0.075 00	0.150 00	0.000 19	0.004 01
−0.5	0.249 68	0.062 529	0.126 145	0.062 50	0.125 00	0.000 029	0.001 145
−0.6	0.199 96	0.050 0025	0.100 1879	0.050 000 00	0.100 0000	0.000 0025	0.000 1879
−0.7	0.15	0.0375 000 615	0.075 009 205	00.037 500 00	0.075 0000	0.000 000 061	0.000 0092
−0.75	0.125	0.031 250 0035	0.062 500 78	0.031 250 000	0.062 5000	0.000 000 0035	0.000 000 78
−0.8	≃0.1	≃0.025 000 001 03	≃0.050 000 0192	≃0.025	≃0.05	≃0.000 000 001 03	≃0.000 000 0192

On the λ_± line, the variation of Δ and (or T) in the limit H → ∞ can be explained in a simple way. In the limit H → ∞, we can safely assume that there are no particles with spin −1, or in other words, N₋ = 0. Thus, q is equal to m. In this limit, the entropy of the system (per particle) can be written as s = −(1 − m)ln(1 − m) − m ln(m), which has a maximum at m = 1/2. The energy per particle, in this limit, turns out to be → {(Δ − H)/2 − (K + 1)/8}. In order for the energy (per particle) to be finite on the transition line, Δ should also increase linearly with H. We indeed see a linear variation of Δ with H on the λ_± line. If we use the variation of Δ as approximated numerically, i.e. Δ ≃ (K + 1)/2 + H, we can estimate the saturation value of → _sat ≃ (K + 1)/8. Using these values in the expression for calculating the temperature (equation (29)), it can easily be shown that the saturation of T is T_sat ≃ (K + 1)/4. Hence, the saturation values of and T will become zero for K = −1.

We measure the widths of the wings of energy (and temperature), denoted by _w (and T_w). We also list the saturation values of and T and the widths of the wings (_w and T_w) in table 2. We plot the BEP and the width of the wings (_w and T_w) in figure 8. We note that the BEP tends to = T = 0 as K → −1. The widths of the wings are also found to decrease exponentially and tend to zero as K tends to −1. From all the above observations, it is clear that at K = −1, the width of the wings vanishes and the BEP reaches = T = 0. Thus, for K ⩽ −1, there is no phase transition in the non-zero H plane for a finite T.

Zoom In Zoom Out Reset image size

4.2. The attractive Blume–Emery–Griffiths model

The attractive BEG model has been previously studied for the microcanonical ensemble in [12], in the (T − Δ) plane. The full phase diagram (T − Δ − H) has not been studied before for the microcanonical ensemble, to the best of our knowledge. In this section, we present the results for the attractive BEG model in the (T − Δ − H) space. We find that in the range 0 < K < 3, the phase diagram is similar to the case of 0 > K ⩾ −0.0828. For K > 3, the λ line is truncated at the first-order line at a CEP. The first-order line continues to exist in the paramagnetic region and becomes a surface in the (T − Δ − H) space which separates the two paramagnetic phases P1 and P2 (as discussed previously in section 3.2) and the wings no longer exist.

For small positive values of K, the variation of is monotonic with H on the λ_± lines, similarly to the case of negative K. For large positive K (⩾1), however, the variation in is non-monotonic on the transition line as shown in figure 9(a). This can be understood by separating the expression for into two parts, = ₁ + ₂, where ${{\epsilon}}_{1}={\Delta}q-\frac{1}{2}{m}^{2}-Hm$ and ${{\epsilon}}_{2}=-\frac{K}{2}{q}^{2}$. We note that the variation of ₁ remains similar for small as well as large K, however, the variation of ₂ is different for small and large K. It decreases with H for small K but increases with H for large K (see figure 9(b)). The variation in ₂ is mainly due to the variable q, which itself shows this behavior. In ₁, we also have the variable q, but it appears with other terms. ₁ does not change its qualitative behavior when we change K. For small K, since both ₁ and ₂ decrease with H, the sum also decreases with H. For large K, there is competition between ₁ and ₂. In the small H regime, the variation in ₂ dominates, which gives rise to an increase in with H. For large H, the variation in ₁ starts to dominate and decreases with H. For very large H, ₁, _2, and all finally saturate. The saturation values of follow a similar relationship with K to that obtained for negative K. The variation of T also shows a similar non-monotonic behavior in the same range of K. We made similar observations for the canonical ensemble in section 3.2.

Zoom In Zoom Out Reset image size

The inequivalence of different ensembles in the BEG model has been reported earlier in [12, 24] in the absence of a magnetic field. In the (T − Δ) plane, while the λ line equation is same in both the ensembles, the first-order line and the multicritical points are known to be located differently [38–41]. It has been reported that the TCP and other multicritical points are different for canonical and microcanonical ensembles for a given value of K [12, 24].

In this work, we have examined all three continuous transition lines (λ, λ₊, λ₋) and the first-order surfaces. We found that not only the multicritical points but the continuous transition lines λ₊ and λ₋ were also different in the two ensembles. In fact, the ensemble inequivalence of the two ensembles can be seen as a consequence of this inequivalence. For K = 0, which corresponds to the Blume–Capel model, in figure 10, we plot the locus of the λ₊ line in two ensembles and one can see that they are different (the λ₋ line also behaves in a similar way). We plot the product of βΔ on the λ₊ line as a function of H for both ensembles, and note that for H → 0, these lines meet at different points, which is the TCP of their corresponding ensembles. For the canonical ensemble, these λ_± lines meet at (βΔ)_TCP ≃ 1.3863, while for the microcanonical ensemble, these lines meet at (βΔ)_TCP ≃ 1.3998 (see figure 10). We also note that the λ₊ lines for the two ensembles become close to each other for large H. We plot the difference in the value of βΔ for the two ensembles for a given K, and plot it as a function of H in the inset of figure 10. We note that this value decreases exponentially to zero as H becomes large.

Zoom In Zoom Out Reset image size

For non-zero K, the λ_± lines meet the λ line at the TCP. This topology persists for 0 ⩾ K ⩾ −0.1838 in the canonical ensemble, whereas for the microcanonical ensemble, this topology occurs for 0 ⩾ K ⩾ −0.0828. As K decreases further (for the canonical ensemble, −0.1838 < K < −1 and for microcanonical ensemble, −0.0828 < K < −1), the λ_± lines move inside the ordered region and meet at the BEP in the H = 0 plane. Interestingly, we find that the difference in the positions of the BEP and the CEP in the two ensembles decreases with decreasing K and for K = −1, the two ensembles become equivalent. In figure 11(a), we plot the value of βΔ at the BEP for both ensembles. We note that the value of β_BEPΔ_BEP for the two ensembles becomes closer as K → −1. In the inset of figure 11(a), we also plot the difference in the value of β_BEPΔ_BEP for the microcanonical and canonical ensembles, and note that this difference decreases exponentially as K → −1. Thus, for K ⩽ −1, we find that there is no ensemble inequivalence in the H = 0 plane.

Zoom In Zoom Out Reset image size

We have shown in sections 3.1 and 4.1 that for K ⩽ −1, there is no phase transition for a finite magnetic field in either of the ensembles and hence there are no wings. Thus, there is also no inequivalence in the H ≠ 0 plane. For K > −1, however, we do have wings, and the continuous transition lines λ₊ and λ₋ meet at their corresponding TCPs or BEPs for the canonical and microcanonical ensembles in the limit H → 0. Thus, the critical lines in the H plane are different for the two ensembles for K > −1. In figure 11(b), we plot the value of βΔ on the continuous transition line for K = −0.3 as a function of H for both ensembles. We note that the two lines are different for small H, however, these lines tend to meet each other for large H. We measure the difference between the values of βΔ for the two ensembles for a given K, and plot it as a function of H in the inset in figure 11(b). We note that this difference almost exponentially reaches zero as H increases. Thus, in the limit H → ∞, these critical lines for both ensembles become equivalent.

From the above discussion, it is clear that ensemble inequivalence is observed for K > −1 with small H values, however, for K ⩽ −1, the phase diagrams in the two ensembles become equivalent. In previous literature [39, 42–44], where ensemble inequivalence with long-range interactions were studied, it was found that whenever the two ensembles (either micro-canonical/canonical or canonical/grand-canonical) have a continuous transition, the transition occurs at the same point. The phase diagrams of the two ensembles can, however, be different from each other when the phase transition becomes first-order in one of the ensembles. This kind of behavior is observed in many systems such as the spin-1 BEG model [42, 43], the ABC model [39, 44], etc.

However, this is not always true. Even the continuous transition points can differ in the two ensembles. For example, in a generalized ABC model in [45], the canonical and grand-canonical ensembles are found to exhibit a second-order phase transition at different points in the phase space. In [38], a general statement is provided to check the possibility of ensemble inequivalence for continuous transition using Landau theory. The transition is observed for a system undergoing a phase transition governed by some order parameter, 'μ₁' (for example), in a given ensemble. This parameter can be the average magnetization in the case of a magnetic transition, or the difference in the density of the two phases for a liquid–gas phase transition. The model is then considered within a 'higher' ensemble, where a certain thermodynamic variable, denoted by 'μ₂', is allowed to fluctuate (within the 'lower' ensemble, 'μ₂' was kept at a fixed value). In the case where μ₂ denotes energy, the two ensembles would correspond to the canonical and micro-canonical ensembles, while in the case where μ₂ is the particle density, they would correspond to the grand-canonical and canonical ensembles. The system is thus described by the Landau free energy denoted by f(μ₁, μ₂). In [38], the authors of that paper found that if f(μ₁, μ₂) = f(−μ₁, μ₂), the two ensembles were equivalent when any of them showed a continuous transition. If, on the other hand, f(μ₁, μ₂) = f(−μ₁, −μ₂), the system showed ensemble inequivalence, even for continuous transition.

In our case, μ₁ is the magnetization m, and μ₂ is the energy of the system. The lower ensemble in our case is thus the microcanonical ensemble and the higher one is canonical. If we check the above symmetries, we note that neither of the conditions studied in [38] is satisfied. When we add a magnetic-field term, the symmetry of the problem is broken and we find that we have ensemble inequivalence, even when the two ensembles show a second-order transition.

The repulsive and attractive BEG model has previously been studied in the canonical and in microcanonical ensembles. This model is known to exhibit many multicritical points along with first- and second-order transition lines. The model was previously studied in the (T − Δ) plane [1, 12–21, 24–34] and ensemble inequivalence was reported [12, 24]. The full phase diagram in the (T − Δ − H) space was only studied for the attractive BEG model in the presence of an external field in the canonical ensemble [14]. We revisited the model in order to study the full phase diagram in the (T − Δ − H) space in both ensembles on a complete graph. Although we explored the phase diagram for the entire range of K, we mainly focused on the repulsive BEG model. We found that for small negative values of K, the model exhibits a TCP where the λ_± lines meet. For the canonical ensemble, the range of K for such a topology was 0 ⩾ K ⩾ −0.1838, whereas for the microcanonical ensemble it was 0 ⩾ K ⩾ −0.0828. As K decreases further, the wings meet inside the ordered phase at a BEP (for the canonical ensemble, −0.1838 > K > −1 and for the microcanonical ensemble, −0.0828 > K > −1). This point was identified as an ordered CP in an earlier study [24]. We also observed that as K → −1, the width of the wings decreases. At exactly K = −1, the wing width in temperature becomes zero and the CEP and the BEP reach T = Δ = 0 in both the canonical and microcanonical ensembles. For K > −1, we observe that the λ_± lines are different in the two ensembles and that they meet at different multicritical points in the H = 0 plane. For K ⩽ −1, we find that there is no phase transition in the H plane for both ensembles.

The absence of a transition in the H plane for K ⩽ −1 can be argued by considering the energy. The energy of the system in terms of the order parameters can be written as: ${\epsilon}=-\frac{1}{2}\left({m}^{2}+K{q}^{2}\right)+{\Delta}q-Hm$. For low temperatures, we can take m ≈ q. In that case, only for K > −1 can the first term lower the energy and take over the entropy at sufficiently small temperatures. Hence, the transition in the H plane is only likely for K > −1.

Disorder, in general, is known to smooth the first-order transition, and has been known to convert a TCP into a BEP [36]. Here, we studied a pure BEG model. We found that the competition induced by negative K affected the phase diagram in a similar manner. It would be interesting to see if there is a similarity in the phase diagram of the two problems, even in finite dimensions.

Using the three-dimensional Blume–Capel model, it was shown in [22] by constraining the number of vacancies that such a constraint modifies the behavior near the tricriticality and thus some of the critical exponents are renormalized, although the universality class remains unchanged. This gives rise to a discrepancy between the constrained and unconstrained systems which acts as an ensemble inequivalence. It might also be interesting to investigate whether such discrepancies exist in the full phase diagram ((T − Δ − H) space).

Also, earlier work on frustrated BEG on bipartite lattices [25, 28] showed that two new ordered phases, namely, the antiquadrupolar and ferrimagnetic phases, occur for repulsive BEG. These phases are not possible on a complete graph. Hence, it would be interesting to study the effects of large negative K on these two phases, especially in the mean field limit on a bipartite lattice.

Raj Kumar Sadhu acknowledges National Institute of Science Education Research, Bhubaneswar (India) for funding the visiting research programme during the initiation of this project.

In order to solve the Hamiltonian (equation (1)), we take the non-interacting Hamiltonian: $H{\sum }_{i}{S}_{i}-{\triangle}{\sum }_{i}{S}_{i}^{2}$, with the probability measure

$$\begin{equation}P\left(1\right)=\frac{{\text{e}}^{\beta \left(H-{\Delta}\right)}}{\left.1+2\enspace {\mathrm{e}}^{-\beta {\Delta}}\enspace \mathrm{cosh}\enspace \beta H\right)}\end{equation} \tag{ A1 }$$

$$\begin{equation}P\left(-1\right)=\frac{{\text{e}}^{-\beta \left(H+{\Delta}\right)}}{\left.1+2\enspace {\mathrm{e}}^{-\beta {\Delta}}\enspace \mathrm{cosh}\enspace \beta H\right)}\end{equation} \tag{ A2 }$$

$$\begin{equation}P\left(0\right)=\frac{1}{\left.1+2\enspace {\mathrm{e}}^{-\beta {\Delta}}\enspace \mathrm{cosh}\enspace \beta H\right)}.\end{equation} \tag{ A3 }$$

The scaled cumulant generating function is:

$$\begin{align}\hfill {\Lambda}\left({k}_{1},{k}_{2}\right)& =\underset{N\to \infty }{\mathrm{lim}}\enspace \frac{1}{N}\enspace \mathrm{log}{{< }{\mathrm{e}}^{N\left({x}_{1}{k}_{1}+{x}_{2}{k}_{2}\right)}{ >}}_{P}\hfill \\ \hfill & =\mathrm{log}\left(1+2\enspace {\mathrm{e}}^{{k}_{2}-\beta {\triangle}}\enspace \mathrm{cosh}\left({k}_{1}+\beta H\right)\right)-\mathrm{log}\left(1+2\enspace {\mathrm{e}}^{-\beta {\triangle}}\enspace \mathrm{cosh}\enspace \beta H\right).\hfill \end{align} \tag{ A4 }$$

The rate function R for the non-interacting Hamiltonian can then be evaluated using the Gärtner–Ellis theorem [37] and is given by:

$$\begin{align}\hfill R\left({x}_{1},{x}_{2}\right)& =\underset{{k}_{1},{k}_{2}}{\mathrm{sup}}\left[{x}_{1}{k}_{1}+{x}_{2}{k}_{2}-{\Lambda}\left({k}_{1},{k}_{2}\right)\right]\hfill \\ \hfill & =\underset{{k}_{1},{k}_{2}}{\mathrm{sup}}\left[{x}_{1}{k}_{1}+{x}_{2}{k}_{2}-\mathrm{log}\left(1+2\enspace {\mathrm{e}}^{{k}_{2}-\beta {\triangle}}\enspace \mathrm{cosh}\left({k}_{1}+\beta H\right)\right)\right]\hfill \\ \hfill & \quad +\left.\mathrm{log}\left(1+2\enspace {\mathrm{e}}^{-\beta {\triangle}}\enspace \mathrm{cosh}\enspace \beta H\right)\right].\hfill \end{align} \tag{ A5 }$$

Minimizing the above equation w.r.t. k₁ and k₂ gives the following relations:

$$\begin{equation}{\Rightarrow}{x}_{1}=\frac{2\enspace {\mathrm{e}}^{{k}_{2}^{{\ast}}-\beta {\Delta}}\enspace \mathrm{sinh}\enspace {k}_{1}^{{\ast}}}{1+2\enspace {\mathrm{e}}^{{k}_{2}^{{\ast}}-\beta {\triangle}}\enspace \mathrm{cosh}\enspace {k}_{1}^{{\ast}}}\end{equation} \tag{ A6 }$$

$$\begin{equation}{\Rightarrow}{x}_{1}=\frac{2\enspace {\mathrm{e}}^{{k}_{2}^{{\ast}}-\beta {\Delta}}\enspace \mathrm{cosh}\enspace {k}_{1}^{{\ast}}}{1+2\enspace {\mathrm{e}}^{{k}_{2}^{{\ast}}-\beta {\triangle}}\enspace \mathrm{cosh}\enspace {k}_{1}^{{\ast}}}\end{equation} \tag{ A7 }$$

where ${k}_{1}^{{\ast}}$ and ${k}_{2}^{{\ast}}$ are the minimums of k₁ and k₂. This implies

$$\begin{equation}\frac{{x}_{1}}{{x}_{2}}=\mathrm{tanh}\enspace {k}_{1}^{{\ast}}.\end{equation} \tag{ A8 }$$

The interacting part of the Hamiltonian is $-\frac{1}{2}{x}_{1}^{2}-\frac{K}{2}{x}_{2}^{2}$. The full rate function of the total Hamiltonian can now be obtained by making use of the tilted LDP. The tilted LDP allows us to generate a new LDP from an old LDP by a change of measure.

Let W_n be a sequence of random variables taking their values from $\mathcal{H}$ and a subset A of $\mathcal{H}$. We define the probability measures

$$\begin{equation*}{Q}_{n,{\Phi}}=\frac{1}{{Z}_{n}}{\int }_{A}{\mathrm{e}}^{\left[n{\Phi}\left(x\right)\right]}{P}_{n}\left\{{W}_{n}\in \mathrm{d}x\right\}\end{equation*}$$

where Z_n denotes the normalizing constant. Here, P_n is the probability measure on $\mathcal{H}$ which satisfies the LDP with a rate function of I and Φ is a continuous function mapping $\mathcal{H}$ into $\mathcal{R}$ which is bounded from above. According to the tilted LDP, the sequence of probability measures {Q_n,Φ, n ∈ N} now satisfies the LDP on $\mathcal{H}$ with the rate function

$$\begin{equation*}{I}_{{\Phi}}\left(x\right)=\left[I\left(x\right)-{\Phi}\left(x\right)\right]-\underset{y\in \mathcal{H}}{\mathrm{inf}}\left\{I\left(y\right)-{\Phi}\left(y\right)\right\}.\end{equation*}$$

In our system, the rate function R(x₁, x₂) can hence be tilted to give the full rate function as

$$\begin{equation}I\left({x}_{1},{x}_{2}\right)={x}_{1}{k}_{1}^{{\ast}}+{x}_{2}{k}_{2}^{{\ast}}-{\Lambda}\left({k}_{1}^{{\ast}},{k}_{2}^{{\ast}}\right)-\frac{\beta {x}_{1}^{2}}{2}-\frac{\beta K{x}_{2}^{2}}{2}-\underset{{k}_{1},{k}_{2}}{\mathrm{inf}}\left[R\left({k}_{1},{k}_{2}\right)-\frac{\beta {k}_{1}^{2}}{2}-\frac{\beta K{k}_{2}^{2}}{2}\right],\end{equation} \tag{ A9 }$$

the minimization of which w.r.t. x₁ and x₂ gives the following free-energy functional:

$$\begin{align}\hfill \tilde {f}\left(m,q\right)& =\frac{\beta {m}^{2}}{2}+\frac{\beta K{q}^{2}}{2}-\mathrm{log}\left(1+2\enspace {\text{e}}^{\beta \left(Kq-{\Delta}\right)}\enspace \mathrm{cosh}\enspace \beta \left(m+H\right)\right)\hfill \\ \hfill & \quad +\mathrm{log}\left(1+2\enspace {\text{e}}^{-\beta {\Delta}}\enspace \mathrm{cosh}\enspace \beta H\right)\hfill \end{align} \tag{ A10 }$$

where the minimums are denoted by ${x}_{1}^{{\ast}}=m$, ${x}_{2}^{{\ast}}=q$, which give equation (8) at the fixed points in section 3.