Phase diagram of the frustrated FCC antiferromagnet from effective-field theory

The phase diagram of a face-centred cubic (FCC) antiferromagnet is calculated from the effective field theory (EFT) of Honmura and Kaneyoshi taking into account not only the effect of interaction with nearest neighbours, J1, but also the effect of second neighbours, J2. The phase diagram for the nearest neighbour case away from the triple point, which in our calculations is predicted to be at H = 4 and T = 0, is close to cluster variation method (CVM) and Monte Carlo (MC) results. Similar to MC and CVM predictions, we observe that the increasing second neighbours interaction pushes the triple point towards zero field. Our calculations also show that for α = −J2/J1 = 0.3, the triple point merges with the transition point of the L10 phase, one of the ground states, at H = 0 and changes the nature of phase transition from first- to second-order, in full agreement with Monte Carlo predictions. The phase diagram with the effect of second neighbours calculated for several values of α are in good agreement with available MC and CVM predictions.


Introduction
Frustration, a characteristics of antiferromagnets, significantly changes the collective ordering behaviour of a system. For example, the ferromagnetic triangular Ising model with nearest neighbour interaction, J, undergoes a continuous phase transition at finite temperature, T c ≈ 3.4J (in this article we assume k B = 1), while the transition temperature for its antiferromagnetic counterpart is T N = 0. This dramatic shift of transition temperature is due to the effect of geometric frustration in the triangular lattice. Frustration also plays a significant role in the case of face-centred cubic (FCC) lattice because the network of tetrahedra in an FCC lattice causes the lattice to be 'fully frustrated'. Since the Hamiltonian of a binary alloy can be exactly mapped to an Ising antiferromagnet [1], an FCC antiferromagnet is doubly important from materials science point 1 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. of view. There are numerous binary alloys that crystallise in FCC lattices. Therefore, methods developed to study frustrated antiferromagnets can be exploited to study ordering in binary alloys too.
The first phase diagram for a CuAu alloy, i.e. an FCC system, was studied by Shockley [2] using the Bragg-Williams conventional mean-field (c-MFT) approximation. As we will discussed extensively in following sections, the c-MFT method fails to reproduce the correct topology of the phase diagram. Transition temperatures along phase boundaries of c-MFT are two times larger than experimental CuAu. Although Shockley attributed the disagreement with experiment to neglecting of next (or higher)-nearest neighbours interactions, it was later found that c-MFT disagrees with numericallycalculated Monte-Carlo results with only nearest neighbours. The main reason for the failure of conventional mean-field theory (c-MFT) in frustrated systems is 'the unnatural' averaged effective field felt by a spin [3]. The subsequent study of the CuAu system was the quasi-chemical calculation of Li [4]. The quasi-chemical phase diagram was improved in comparison to c-MFT. However, there are missing features in the phase diagram which makes it unsatisfactory. The first significant advance in calculation of phase diagram was the earliest calculations of a binary FCC alloy [5][6][7] using the cluster variation method (CVM) of Kikuchi [8]. In these calculations, a tetrahedron was used as the basic cluster of the CVM. They argued that because of the frustration effects, the point, the bond (Bethe), or the triangle approximation are not satisfactory and one must consider at least a tetrahedron as the basic cluster in order to obtain sensible results [9]. Denoting ↑ spins as A atoms and ↓ spins as B atoms, phase diagram of an FCC antiferromagnet depending on the chemical potential is made of A(B), A 3 B(B 3 A), and AB. In the material science literature, AB phase is called L1 0 , and A 3 B is called L1 2 , please see figures 1(a)-(c) for schematic illustrations. Other orderings on the lattice are possible, e.g. DO 22 which has the same ground-state energy as L1 2 in a near-neighbour model.
Binder's paper [10] was the first calculation of the phase diagram from Monte-Carlo (MC) which sparked a controversy. According to Binder's calculation [10], the phase boundaries of L1 0 and L1 2 meet at T = 0 where they neighbour the disordered-phase at a 'triple point'. Real-space renormalization group analysis also predicted the triple point to be at the field H = 4, the temperature T = 0 [11]. Contrary to MC results of Binder, the tetrahedron CVM calculations [5][6][7] predict a triple point at a finite temperature, T = 1.6. Intrigued by Binder's results, Sanchez et al [12] extended the basic cluster of CVM to a larger cluster of a tetrahedron with an octahedron (TO). The temperature of triple point from their calculation was predicted to be about 1.2. A question could be raised: if one increases the cluster size, will the triple point will be calculated as zero as Binder predicted?
Ground states of an FCC binary alloy (antiferromagnet) with nearest neighbour interaction are infinitely degenerate. The reason is that antiphase boundaries (APBs) can be created in L1 0 and L1 2 phases with no energy penalty. Presence of APBs in a system can make the interpretation of MC results rather difficult. This is because a straightforward averaging over the entire system will lead to a vanishing long-range order parameter. However, when a small amount of positive (ferromagnetic) next-nearest neighbour (NNN) interaction is added, APBs will disappear and ground states become unique up to symmetries.
Gahn [13,14] noted this fact and observed that the triple point could be accurately determined by adding a weak next nearest neighbour (NNN) interaction. Lebowitz and coworkers [15] calculated the position of the triple point as a function of NNN interaction. In the limit of vanishing NNN interaction, a linear extrapolation of their results gave a triple point at T ≈ 1. However, it is not obvious if the triple point would behave linearly as a function of NNN interaction. In a subsequent MC calculation, Diep et al [16] employed the Edwards-Anderson order parameter [17] which measures time correlations to circumvent the problem of determination of long-range order parameter. The triple point from their results was predicted to be at T = 1.00 ± 0.1 confirming the results of Lebowitz et al [15]. The discrepancy between Binder's phase diagram and the CVM results was attributed to the appearance and the fluctuation of APBs around the triple point [18,19].
In this paper, the phase diagram of an FCC binary alloy (antiferromagnet) is calculated using the effective field theory of Honmura and Kaneyoshi [20] (HK). As discussed above, the effect of NNN interactions on the phase diagram is substantial and will be considered here. A brief account of the HK theory and its derivation for an FCC system will be given in section 2. In section 3, we will contrast our results with available MC and CVM calculations and discuss accuracy and effectiveness of the HK theory. The phase diagram for the nearneighbour case was also determined from MC calculations and an attempt was made to understand structure emerging at the triple point, which will be discussed with more details in section 3.

Effective field theory
In the language of a magnetic system, the Hamiltonian, H, for the Ising model with the nearest-neighbour interaction, J 1 , and the next-nearest-neighbour interaction, J 2 , is where, H is the magnetic field, stands for summation over the nearest neighbours (NNs), and for NNNs each ij pair being counted only once. The Hamiltonian in the above equation is a function of, Assuming a cluster of spins, N, being separated from the rest of the system, the effective field theory of HK is based on the exact identity of Callen [21]-Suzuki [22] for the cluster from which the thermal average of an observable, O N , can be calculated from, where Tr N is a partial trace over the cluster of spins and β = 1/k B T. The observable O N here is a function of spins in the cluster and . . . stands for the canonical ensemble average. A cluster of two spins, spin 1 and 2 in figure 1(a), is considered in an FCC lattice to study the effect of NNN interactions. Rationally, increasing the size of the cluster should lead to a more accurate approximation. In addition to a cluster of two spins, the transition temperature of L1 0 phase at zero field, H = 0, in the case of NN interaction, α = −J 2 /J 1 = 0, will be analysed using clusters of three and four spins (e.g. spins 1 to 4 in figure 1(a)).

Two-spin cluster
The Hamiltonian for a two-spin cluster is with, where, N i is the set of NNs and M i is the set of NNNs of site i. For a cluster of two-spins, the number of NNs is 11 and the number of NNNs is 6, and a total cluster of 20 sites must be considered. The average magnetization for the spin i (see, for example, reference [23]) can be calculated from which, after taking the trace, can be written as where,J 1 = βJ 1 ,H = βH, andÃ i = βA i . How the field, H, is taken into account in the current formalism is similar to reference [24][25][26]. Making use of the identity exp(aD x + bD y )g(x, y) = g(x + a, y + b), where D x = ∂ ∂x is the differential operator, the above equations can be expressed as with Equations (9) and (10) can be further expanded using the exact van der Waerden identity for a two-state spin (S i = ±1), exp(aS i ) = cosh(a) + S i sinh(a). This will make the calculations intractable. However, one can progress by employing the decoupling approximation [27]. The decoupling (or Zernike) approximation ignores higher-order correlations and assumes the ensemble average of a correlation function can be decomposed to the product of averages of it constituents such as Although the above approximation neglects correlations between spins, it satisfies the 'hard-spin' condition, (S i ) 2 = 1.
In c-MFT, the effective field created from unequal magnetizations of neighbouring spins (of a given spin) inappropriately removes the frustration. Nevertheless, in a frustrated system, a spin feels an effective field from all the neighbouring spins S i = ±1. And this is the reason that a frustrated system can remain disordered even at low temperatures; the net effective field from its neighbours is zero [3,28]. Therefore, it seems to be imperative to take into account the 'hard-spin' condition for frustrated systems.
Once the approximation of equation (13)  . We obtain equations of state for L1 0 given by where ilarly, the equations of state for L1 2 can be given by And the equation of state for the disordered phase is The definition of coefficients δ i n,m , (i = A, B) and γ j (H, T, α) are similar to λ in the L1 0 case. Solutions of equations (14)-(18) will be discussed in section 3.

Three-spin cluster
The Hamiltonian contribution from the equilateral triangle three-spin cluster is where and N i is the set of 10 NNs of site i, excluding the two in the cluster itself. The site magnetization is calculated from the partial trace over the cluster, Following a procedure similar to section 2.1, the equation of state can be obtained. The transition at zero field for L1 0 is of interest here. Following the procedure in previous, one can equations of state are Omitted details are presented in appendix A.

Four-spin cluster
The Hamiltonian contribution from the tetrahedron four-spin cluster is where The number of external NNs is 9 for a cluster of four-spins. Similar to section 2.2, equations of state for L1 0 phase at the zero field, H = 0, can be obtained which are Further details are given in appendix B The coefficients, λ α i, j , δ α i, j , γ i , φ α i, j and ψ α i, j , in the above equations are rather too lengthy to be reproduced here.

Calculation of free energy
Determination of the phase diagram requires the knowledge of Helmholtz free energy. While the HK theory leads us to equations of state, it does not provide any means to calculate the free energy [30,31]. Once solutions of the equations of state are known, the free energy, F, is calculated by integration of the Gibbs-Helmholtz relation, along a constant field parallel to the T-axis, where U is the internal energy and M is total magnetization [32]. The integral is taken from a high temperature limit: the absolute value of free energy is not known, but nor is it needed, because only free energy differences are important for the phase diagram. At high temperature all solutions collapse to the disorderedphase, so the free energy differences between them are zero. The equation has a singularity at zero temperature, so in practice, the free energy is calculated up to a small cut-off (T cut ≈ 0.1) and then extrapolated to zero.

Results and discussion
To obtain the free energy and subsequently the phase diagram, the equations of state presented in the previous section are solved numerically on an uniform grid of temperature and field. Then, around special points, e.g. the triple point, a much denser grid is used. It should be added that for a given set of equations of state at a field and a temperature, there are multiple solutions, some of which are unphysical (e.g. sublattice magnetizations exceeding one). Table 1 shows the calculated transition temperatures for L1 0 at zero field. We observe that as the cluster size increases, the calculated value gets closer to the MC results. The agreement of 4-spin cluster with other theoretical results is rather remarkable; the difference is between 0.85%-2.37%. It should, however, be noted that a large finite size effect in MC calculations was reported [36] for this transition temperature where, depending on the system size, T N can fall in the range 1.69-1.76. It is known that CVM overestimates transition temperatures [5,12,18,34], so the difference there is as expected.  [33]. b MC result of reference [16]. c CVM calculation of reference [12]. d CVM calculation of reference [34]. e SE analysis of reference [35].

Nearest neighbour interactions only (α = 0)
The c-MFT equation of state for the site i is given by where, N i is the NNs of spin i and M i is the NNNs of spin i. For a cluster of 2-spins, the above equations were numerically solved for the case of α = 0. Free energies were calculated following the procedure described in section 2.4. The phase diagram for c-MF can be seen in figure 2(a) where it is compared to the modified MF calculation of Beath and Ryan [29] (BR). Transition temperatures of BR are lower than c-MF, though the overall topology of phase diagram seems to be almost identical. A feature of the phase diagram which is unexplained in the literature is the strong re-entrant behaviour of L1 2 phase at H 12. The modified MF treatment of BR captured the same behaviour. To clarify this, sublattice magnetizations for various fields around H = 12 were examined (typical examples in figure 2(b)). For a given temperature (below the transition temperature), each sublattice magnetization has two solutions. The solution with higher free energy is 'unstable' and is denoted by dotted lines in all of figures related to site magnetization. At zero temperature, the structures must be fully ordered; so the free energy of L1 2 is −H/2, for L1 0 it is −2J 1 while the homogeneous 'disordered' phase is ferromagnetic and has free energy of −H + 6J 1 . Therefore the phase boundaries should intercept the axis at H = 4 and H = 12 where two energies are equal.
The L1 2 sublattice magnetization of the minority (A) site goes to zero at low temperature. But at H = 12, a very small increase of the field stabilises the homogeneous phase, which can be described in the 'L1 2 ' setting with m A = m B = −1. Increasing temperature causes the magnetization to drop, due to the entropy, but surprisingly it reaches a minimum after which the minority magnetization increases again with the increase of temperature! The re-entrance of L1 2 around H = 12 is because of the behaviour of the minority magnetization: in the ground state it is opposite to the majority spin, but with increased temperature it samples states aligned with the majority spins. It is also worth mentioning that sublattice magnetizations having the same sign for H > 12 raises the question whether or not that region should be called an antiferromagnet.
The phase diagram calculated from the EFT is plotted in figure 3(a) and compared with the MC calculations of references [10,16] and the TO-CVM result of reference [39]. Comparing figures 2(a) and 3(a), one can see how strongly c-MFT disagrees with other theories in both topology and magnitude points.
The EFT triple point comes out at zero temperature similar to Binder's [10]. Away from the triple point, the phase diagram is quite close to MC and CVM calculations. To check whether or not increasing cluster size will affect the EFT diagram, we also calculated a number of points on the phase diagram using 4-spin cluster of section 2.3. We observed that the phase diagram is the same as 2-spin cluster with the distance similar to the value reported in table 1.
Another point of view regarding the role of APBs was proposed by Falicov and Chrzan [40]. They introduced an exactly soluble simplified model for an FCC nearest neighbour Ising model and studied APBs. They conclude that APBs appear in MC calculations because of finite size effect. According to their results, the 'excitation' gap for an APB depends on field and system size and such APBs should not appear in an infinite system.
The points H = 4 and H = 12 have been called 'superdegenerate' points [15,41], because there is a finite ('residual') entropy at zero temperature, i.e. violating the third law of thermodynamics. While this is true, we note that for any H, α = 0 the ground state is degenerate. More precisely, the stable L1 0 phase for |H| < 4 and L1 2 for 4 < |H| < 12 are degenerate with respect to the same phases containing any APBs. For a finite system with linear size of L, there can be 2 L degenerate states by creation of APBs. Moreover, interfaces between L1 0 and L1 2 have zero energy, so any two-phase microstructure is also degenerate with the ground state-an example is shown in figure 3.
No ordered phase was found from the EFT at the superdegenerate points; the solutions are disordered all way from the cutoff to finite temperature. However, a small deviation from these points, e.g. H = 3.95 or H = 11.95, shows two types of solutions, as expected. Considering the fact that zero temperature L1 0 and L1 2 intersects at H = 4, the phase boundaries are extrapolated from each side to meet at the super-degenerate points.

Ground state degeneracy and antiphase boundaries at α = 0
A number of the discrepancies between methods at α = 0 can be related to the ground state degeneracies (table 2). Both L1 0 and L1 2 structures can have APBs with zero energy, which will contribute a configurational entropy even at T = 0. For non-zero α, these have an extensive energy, so will only appear in MC simulations due to finite-size effects.
The 'superdegenerate' points have even more degeneracy. At H = 4 there is degeneracy between L1 0 and L1 2 and either structure containing APBs, and any mixture of L1 0 and L1 2 has zero (100) interfacial energy. The This means that at H = 12 all 'L1 2 ' structures between m A = 1 (normal L1 2 ) and m A = m B = −1 (homogeneous solution) are degenerate. At H = 12 the c-MFT averages over these degenerate states giving m B = 0, while for H = 12 − δH the ground state is non-degenerate so m B = +1, and for H = 12 + δH we get m B = −1. At finite temperatures, H = 12, the degenerate m B states are sampled with no energy penalty, giving disorder on the minority L1 2 sublattice only. This gives L1 2 a higher entropy than the homogeneous phase whereas any entropy gained from disorder m > −1 also incurs an energy penalty in c-MFT. This extra entropy leads to the stabilization of L1 2 above H = 12 at finite T, and is evidenced by the leap in m B from −1 towards zero in figure 2(b).
These degeneracies raise difficulties for Monte Carlo simulation to define long-range order parameters for the structures. Whereas we use the long range order to identify a crystal structure, the APBs make this impossible in stochastic simulations. In practice Binder uses either a short ranged, or a two dimensional order parameter, whereas Diep looks at time correlation functions. We ran our of Monte Carlo simula-tions using a Metropolis algorithm on an 8000 site fcc lattice, on a grid of (H, T) points with 0.1 units resolution and 8 × 10 7 stochastic updates, with 1.2 × 10 5 equilibration steps, and located the phase boundaries by peaks in the specific heat (vs temperature) and susceptibility (vs field). The results are in agreement with previous work. A typewriter update scheme gives a similar phase diagram, but tends to sweep antiphase boundaries out of the system. Extrapolation from a NNN model is another option, and it is unclear whether these approaches are equivalent.
The inset of figure 3(a) shows a cross section of our MC calculations along 001 direction annealed at H = 4. One can see a mixture of L1 0 and L1 2 phases in the snapshot where it is tiled with many APBs such that the structure might seem disordered. The long-and short-range order parameter for such a structure is very similar to the disordered phase. This is the reason that the Binder's phase diagram is different from reference [16]; the long-range order parameters corresponding to the conventional L1 2 vanish around H = 4 because of the emergence of such structures but the structure still has a form of order. Therefore, one can consider the Binder's phase diagram as the boundary of the conventional L1 2 phase.

Effect of 2nd neighbour (α = 0)
Introducing the NNN interaction leads to a finite energy for APBs and removes the infinite degeneracy. Now, the L1 0 phase is sixfold degenerate for |H| < 4 and the L1 2 phase is fourfold degenerate for 4 < |H| < 12. At H = 4 the degeneracy is tenfold, sum of either side, and at H = 12 it is fivefold (the substantial reduction of degeneracy makes an accurate interpretation of MC calculations much more feasible). The phase diagram for several values of α is plotted in figure 3(b). While the agreement with different MC results are rather good, the agreement with the CVM calculation of reference [39] is excellent.
Adding a very small amount of α changes the shape of phase diagram; the triple point moves from (H t = 0, T t = 0) | α=0 to (H t = 3.44, T t = 0.97) α=0.025 as can be seen in figure 3(c) where the variation of the triple temperature T t is shown as function of α and compared to the MC values of reference [15]. It is apparent from the figure 3(b) that increasing α pushes the triple point towards zero field, H t → 0. While the triple point calculated from the EFT agrees rather well with the MC results, T t calculated from EFT has a non-linear behaviour as α → 0; contrary to the MC trend, T t goes to zero as α approaches zero. As discussed in sections 3.2 and 3.3, this is related to superdegeneracy and emergence of structures with a non-conventional form of order.
We can see in the figure 3(b) that for α = 0.3, the triple point merges with the transition temperature T c of the L1 0 phase at H = 0 and creates a multicritical point which will be discussed in the following.
For the near-neighbour case, MC calculations and CVM results predict that all phase boundaries are first order. The L1 0 phase in a magnetic field can be mapped to the four-state Potts model in three-dimension [42]. Therefore, the phase boundary between L1 0 and the disordered phase is first-order. Following a similar argument, we found that the L1 2 can be mapped to the four -state Potts model too, and hence it has a first-order transition line with the disordered phase. In addition, the transition from L1 0 to L1 2 can be mapped to the three-state Potts model. It is known that the Potts model with three or more components, q 3, in three dimensions exhibits a firstorder phase transition [43]. Our EFT calculations have first order transitions for both phases and across the whole span of the α = 0 phase diagram. An example is shown in inset of figure 3(c) for the L1 0 phase at H = 0, where the discontinuity in the derivative of the free energy shows that the transition for the nearest neighbour case is first-order.
Adding second-neighbour interactions can change the order of transition. It was discussed above that for α 0.3 the triple point moves toward H = 0 and finally merges with the L1 0 transition point. This changes the nature of the phase transition from first-order to second-order at H = 0 as can be seen in inset of figure 3(c) in full agreement with the MC calculation [15]. Reference [15] reports the point of merger for α 0.35 obtained from linear extrapolations of T t (α) and T c (α) and finding where they cross. It is, however, noted in the paper that merging of T t (α) and T c (α) should have a totally different behaviour than of linear exploration. Therefore, α 0.35 should be considered a very rough estimate.
The phase diagram for α = 1 is plotted in figure 3(d), and compared with the MC calculation of reference [38] and CVM treatment of reference [39]. We see that transition temperatures are higher than MC predictions and closer to the CVM results. For example, the transition temperature of multicritical point at H = 0 is T = 8.21 which should be compared to 7.23 and 7.60 calculated from MC and CVM techniques, respectively. It is usually argued that for large value of α, the FCC lattice can be considered as two decouple interpenetrating SC lattices and c-MFT should suffice to describe the phase diagram [1,39]. However, we see that the c-MFT is still far from other calculations hinting that c-MFT is not satisfactory even for a case as extreme as α = 1. Again, the c-MFT and the EFT predict a re-entrant behaviour for the L1 2 at H 12. Contrary to α = 0 case, the magnetization for both cases here behave as expected. The absence of re-entrance from CVM and MC calculations suggests that the re-entrance of L1 2 phase is possibly an artefact caused by over-stabilization of the site magnetizations. Figures 4(a) and (b) show sublattice magnetizations at a constant field for the L1 0 and L1 2 , respectively. Similar to the c-MFT α = 0 case discussed in section 3.2, the sublattice magnetization of L1 0 for α = 0 has two sets of solutions: 'stable' and 'unstable'. An interesting observation is the disappearance of the unstable solution for α 0.3. We should emphasise that at H = 0 and for α 3, L1 0 always transforms to L1 2 which is a first-order transition, please see figures 3(b) and (d).

Sublattice magnetization
On the other hand, the unstable solution always exists for L1 2 phase as can be seen in figure 4(b). The sublattice magnetizations, in this case, jump discontinuously to the disordered phase at T N , indication a first-order transition. The sublattice magnetizations of the L1 0 and the L1 2 for T N /2 are plotted in figures 4(c) and (d), respectively. For α = 0, the solutions are confined to two regions for either phase. Beyond these regions, there is only the disordered solution. In contrast, for α = 0.3 and α = 1.0, the solutions extend to a much larger field; for α = 0 the L1 0 solution stretches away from H = 4. In the case of L1 2 a nonzero α causes the solution to exist well beyond H = 12. This is the cause for the reentrant behaviour for α = 1.0 discussed in section 3.4, though a much milder re-entrance also exists for lower values of α, see figure 3(b).

Conclusion
The phase diagram of an FCC antiferromagnet, a fully frustrated system, is calculated using the effective field theory (EFT) of Honmura and Kaneyoshi (HK); the effect of first and as well as second neighbours are considered in this study. We observed that the predicted phase diagrams are very different from c-MFT and comparable to MC and CVM results away from the triple point for the nearest-neighbour (α = 0) case,.
We ran our MC calculations and explain why the EFT phase diagram for the case of α = 0 is similar to the Binder's calculations and differs from reference [16]; around H = 4, due to superdegeneracy, a set of fluctuating structures emerge which are mixture of L1 0 and L1 2 with many APBs and have different order parameters than the conventional L1 2 , meaning that unconventional order parameters are required to identify that region. The long-range order parameter of these structure is similar to the disordered phase. Therefore, the Binder's phase diagram should be considered as the boundary of L1 2 phase with conventional long-range order parameter.
In EFT calculations, because of the definition of L1 2 , emerging structures around H = 4 are excluded. Consequently, the EFT is capable of determining the phase boundaries for the L1 2 phase with conventional long-range order parameter. The very good agreement between the EFT and the Binder's phase diagram show that the EFT is successful in capturing the frustration effects in FCC lattice.
Curious reentrant behaviour of the L1 2 phase is seen beyond H = 12. This occurs because there are multiple degenerate states involving spin-flips on the minority sublattice. These states can be accessed at zero energy cost but finite entropy gain at finite temperature-a curious situation where the spontaneous symmetry breaking is driven by entropy. M i j is the set of NNs shared by site i and j. The f i functions are f 1 =e 6J 1 2 cosh(Y 34 )(sinh(2H + X 12 ) + e 2J 1 sinh(Y 12 )) + 2 sinh(Y 12 ) cosh(2H + X 34 ) where Z 0 = x + y + z + w, Z 34 = x + y − z − w, X 12 = x + y, X 34 = z + w, Y 12 = x − y, and Y 34 = z − w. Applying the decoupling approximation of equation (13) and the boundary conditions for L1 0 (m A = = , m B = = ) to equation (B.1), one can get the equations of state given by equations (24) and (25). Please see reference [44] for a more pedagogical introduction for a multi-spin formalism of EFT.