Phase transitions, order by disorder and finite entropy in the Ising antiferromagnetic bilayer honeycomb lattice

We present an analytical and numerical study of the Ising model on a bilayer honeycomb lattice including interlayer frustration and coupling with an external magnetic field. First, we discuss the exact $T=0$ phase diagram, where we find finite entropy phases for different magnetisations. Then, we study the magnetic properties of the system at finite temperature using complementary analytical techniques (Bethe lattice), and two types of Monte-Carlo algorithms (Metropolis and Wang-Landau). We characterize the phase transitions and discuss the phase diagrams. The system presents a rich phenomenology: there are first and second order transitions, low-temperature phases with extensive degeneracy, and order-by-disorder state selection.


I. INTRODUCTION
The continuous exploration of frustrated spin systems in the last years has been driven by the role of frustration to induce unconventional magnetic orders or macroscopic degeneracy in the ground state with no long-range ordering [1]. However, this macroscopic degeneracy will depend critically of the coordination number and/or the spin representation. For instance, on one hand in the antiferromagnetic (AF) triangular lattice the classical Heisenberg model has a unique ordered ground state [2] while the classical Ising model shows a large degeneracy in the ground state [3]. On the other hand, in the AF kagome lattice, both models (Heisenberg and Ising) present a disordered ground state with macroscopic degeneracy [4]. The interaction with a magnetic field lowers the symmetries of these systems, and may lead to a total or partial reduction of the ground state degeneracy.
In the case of the honeycomb lattice, since it is bipartite, the AF model with nearest-neighbor interactions is not frustrated. Additional interaction terms, for example next-nearest neighbors, are needed to introduce magnetic frustration [5][6]. In the last few years, several works [7][8][9][10][11][12][13] have been published motivated by non-trivial phases found in the material Bi 3 Mn 4 O 12 (NO 3 ) [14]. Experimental evidence shows that this material can be modeled as a weakly coupled bilayer honeycomb lattice where magnetic frustration, suggested by the large negative value of the Curie-Weiss temperature Θ CW = −257K, could play an important role in low-temperature properties. For that reason, in a recent paper we studied the antiferromagnetic bilayer honeycomb lattice in the highly frustrated case for classical spins [15]. Frustration in this model is given by a competition between intralayer nearest-neighbors, and two interlayer couplings (all antiferromagnetic). In that work, we found that due to the the high level of frustration, an external magnetic field induced the selection of non trivial low temperature phases.
In this work, we study the Ising model in the honeycomb bilayer lattice and explore the magnetic properties for the full range of the antiferromagnetic couplings. In order to do this, we resort to a combination of analytical and numerical techniques (Bethe lattice approximation, Metropolis and Wang-Landau Monte-Carlo simulations). Surprisingly, we find a very good agreement between the numerical results and the mean field technique. As we will show in the following sections, the interplay between these methods has proved essential for a thorough study of the model. The paper is structured as follows: first, we introduce the model and present the T = 0 phase diagrams in Sec. II. In Sec. III we present the methods and the order parameters used to study the low temperature behavior of the system, discussing the characteristics, advantages and limitations of each technique. In Sec. IV we present the low temperature phase diagrams for different regimes of the model. We find different types of phase transitions, highly degenerate phases and selection of states by thermal fluctuations. Concluding remarks are presented in Sec. V.

II. MODEL AND T = 0 PHASE DIAGRAM
Let us define the Ising model on the antiferromagnetic bilayer honeycomb lattice as: where r runs over unit cells, r, r ′ denotes interactions within the cell and between nearest-neighbor cells and i is the spin cell index i = A, B, C, D (with sites A, B, C, D shown in Fig. 1). The lattice structure and exchange antiferromagnetic couplings J 1 , J x and J p are shown in Fig. 1. J p is the coupling joining the honeycomb layers (above and below), J 1 , J x are the nearest-neighbors inplane and intraplane couplings respectively. Note that the model in Eq. (1) can be mapped onto an identical one replacing J 1 ↔ J x by exchange of opposite sites in each plaquette: A ↔ C or B ↔ D. This Z 2 symmetry will play an important role in the characterization of the low-temperature phases.
In the particular case of h = 0 and J x = 0, the ground state has long-range Néel order composed of two opposite Néel states in each layer (for instance σ A = σ D = −σ B = −σ C = +1). A non zero value of J x > 0 introduces frustration, which leads to interesting phenomena even at zero temperature. The magnetic properties of this model are controlled then by two factors :the level of frustration and the magnetic field. We first study the zero temperature T = 0 phase diagram. In order to do this, we need to describe the groundstate configurations for an individual plaquette. There are 16 possible states in each 4-spin plaquette. We consider only those with zero or positive magnetisation, and we are left with five types of plaquette arrangements. These are listed in Table I [16] with their energy (E 0 ), magnetisation (m 0 = 1 4 (σ A + σ B + σ C + σ D )) and degree of degeneracy ( D). We also introduce the notation for each configuration that we will use throughout this work: AF stands for antiferromagnetic ordering in each layer, U for uniform order in a layer. Having constructed the different plaquettes configurations, listed in Table I, we now discuss the T = 0 phase diagram for different cases: We show the J x /J p vs J 1 /J p T = 0 phase diagram of the model for the case of zero magnetic field in the upper-left panel of Fig. 2. Interesting features appear on the specific line J x = J 1 < J p /3 where the ground state is highly degenerate because each plaquette can either be in a U 2 or AF 2 configuration. The energy does not depend on J p , and thus each J p pair can be flipped without energy cost. As a consequence the system has macroscopic degeneracy, and therefore a non zero entropy, at T = 0. The determination of the value of this entropy can be easily computed: it is the same as that of a random spin configuration in one of the layers (the spins in the other layer are simply opposite). Therefore, the entropy per spin s = S/N in units of ln 2 (the value in the paramagnetic case) is s = 1/2. In the highly frustrated point J 1 = J x = J p /3, in addition to the degenerate configurations we just described there are two more possible states: all plaquettes with an AF 1 configuration.

h > 0 case
The effect of an external field is summarized in Fig. 2 where we show the h/J p vs T /J p phase diagrams for three different regions of interest for the frustrating relation   We found that there are three possible magnetisation plateaux: at m = 0 (with structure U 2 , AF 1 or AF 2 ), at m = 1/2 (UAF) and the saturation plateau m = 1 (U) for h > h sat . Explicit expressions for the critical fields are given in Appendix A . As we just stated, the m = 1/2 is given by UAF plaquettes. This configuration is highly degenerate: one J p pair of the spins of a plaquette has positive (+) spin, and the other one has two opposite spins. This pair has no particular orientation, so in each unit cell there is one degree of freedom. Therefore, at this plateau there is finite entropy s = 1/4. For J 1 = J x , the m = 0 plateau is simply doubly degenerate, and therefore has s = 0 in the thermodynamic limit. The magnetic field then takes the system from a s = 0 phase to a magnetised highly degenerate one, where the entropy is finite at T = 0. Wang-Landau), which we briefly describe in the following subsections.

A. Bethe lattice
The Bethe lattice (BL) is a mean-field approach that, for first neighbor interactions, is equivalent to the Bethe approximation [17]. From the point of view of correlations, in a simple mean-field calculation no correlations are taken into account, while in Bethe lattice solutions short-range correlations are considered. In particular, for the model studied here, the correlations between spins in a plaquette are taken into account in exact form.
The Bethe lattice consists in the exact solution of a statistical model in the core of a Cayley tree. In order to approximate the bilayer honeycomb lattice, we define a bilayer Cayley tree, similar as the one used previously for different Ising-like models by Hu et al. [18] and Albayrak and co-workers [19], and for self-avoiding walks by Serra and Stilck [20]. The simplest tree-like approximation for our model is a bilayer Cayley tree with the same coordination q = 3 of the honeycomb lattice, and interlayer first and second interactions. It is important to note that the model in this lattice reproduces the correct T = 0 phase diagram (ground state) of the two-dimensional model described in Sec. II. In order to take statistical averages in the central zone of the Cayley tree (the Bethe lattice approach), we define it as a central plaquette, as shown in Fig. 3.
The honeycomb lattice is bipartite, then, in the notation of Fig. 3, the points A and C belong to a sublattice a, and the points B and D to the other sublattice b. In order to describe the different phases we need to define eight partial partition functions (PPF) Z a σAσC , and Z b σB σD each one for the four possible values of the spins corresponding to a sublattice bilayer, (1, 1), (1, −1), (−1, 1), and (−1, −1). The tree-like structure of the lattice allows us to write down recursion re-lations (RR) for the PPF. Once the RR are obtained, the thermodynamic phases are given by their stable fixed points. Where the stability line of two different fixed points are coincident, it represents a second order line. For first order lines we have to calculate the Bethe lattice free energy [21 and 22], and look for the line in the coexistence zone where the free energies of both phases are equal. The PPF, free energy and other details of the calculations are developed in Appendix B.
In general, the stability lines of the several thermodynamic fixed points are given as solutions of sets of nonlinear coupled equations. However, for the important case h = 0, J 1 = J x , the paramagnetic stability line takes the simple form The other curves in the phase diagram were calculated solving numerically the set of exact algebraic coupled equations (see Appendix B). The order parameters and phase diagrams obtained within the Bethe lattice approximation are displayed in the following subsections, which describe the MC techniques.

B. Monte-Carlo simulations
We simulated lattices with periodic boundary conditions in the two directions with Metropolis update [23] (MC-M) and Wang-Landau [24] (MC-WL) methods. In both cases we used a single spin flip algorithm. We exploited the strength of each technique to fully understand the system, as commented below.

Metropolis
We performed MC-M simulations on lattices of N = 4 × L 2 sites (L = 24 − 60). To avoid the problem of low temperature "freezing" of the simulations, we used the the annealing technique, lowering the temperature as T n+1 = 0.9 × T n , from T i /J p = 5 to T f /J p ≈ 0.1. We also averaged results over 100 copies of the simulations, generated from different random seeds. Data was taken in each copy averaging 4 × 10 5 Monte-Carlo Steps (mcs), after discarding 2 × 10 5 mcs for thermalization. We measured the energy per spin, the magnetisation per spin (Eq. (2)), the specific heat per spin and three different order parameters to detect the three possible zero magnetisation arrangements of the plaque-ttes (see Table I) defined as follows, where r runs over unit cells. With the previous definitions, the local order parameters take values −1 or 1 when the plaquette is in the specific configuration (AF 1 , AF 2 or U 2 ), and zero if they are in any of the other two.
In each Monte-Carlo step, the order parameter is calculated for every unit cell. It should be noted that when averaging these parameters for different copies we will take the absolute value of the MC measurement, since otherwise it will average to zero.

Wang-Landau
The Wang-Landau (MC-WL) algorithm [24] has emerged as an efficient Monte-Carlo technique in statistical physics. In the last years, this technique has been applied to a variety of studies of classical statistical models such as the Ising [25] and Potts [26] spin model, Heisenberg ferromagnetic systems [27] and antiferromagnetic frustrated models [28]. In this work, we performed simulations on lattices of N = 4 × L 2 sites (L = 2 − 10). To optimize the convergence of the algorithm we use the modification proposed in Ref. [29]. This algorithm allows the estimation of the energy density of states (DoS) g(E) performing a random walk in energy space. Then, from the DoS we can construct the partition function and obtain thermodynamic quantities like entropy and free energy, which are not easily accessible through conventional Monte-Carlo methods based on Metropolis algorithm. To better characterize the system, we often need to calculate a joint density of states (JDoS) g(E, OP), where OP is an order parameter. This allows us to explore the phases of the system and also to calculate thermodynamic quantities like the Landau free energy [30].
Once obtained g(E, OP), the partition function can be computed as with β = 1/k B T and µ some Lagrange multiplier (for example, µ may be the magnetic field and thus the OP would correspond to the total magnetisation) [31]. From the partition function, we can obtain thermodynamic quantities in the canonical ensemble for all values of β (temperature) and µ. For instance, the mean value of the energy E and the order parameter OP may be calculated as: In addition to the standard averages, it is straightforward to determine some important quantities like Helmholtz's free energy and entropy defined as Finally, to obtain more information on the global behavior of the order parameter around the phase transition, in this work we have calculated the Landau free energy as The study of the Landau free energy will be further discussed in the following section.

IV. RESULTS AND PHASE DIAGRAMS
In this section we explore the low temperature behavior of the model introduced in Eq. (1) combining the three approaches described in the previous section. We first focus on the model in the absence of a magnetic field, where a rich phenomenology is found. Then, we discuss the effect of an external field.
In the absence of the magnetic field, the system presents different phase transitions and selection mechanisms, which we will discuss below.

1.
J1 = Jx: second order phase transitions for broken Z2 symmetry ground states As a typical example, we focus on the J 1 = J p case, where for J x < J p /3 the ground state is AF 2 and for J x > J p /3, AF 1 . We show the transition lines from a paramagnetic to an ordered phase in a J x /J p vs T /J p phase diagram in Fig.4 [32]. These lines were obtained from the Bethe lattice analysis and from MC-M simulations. According to the Bethe lattice technique, these transitions are of second order. We check this with MC-M following the standard procedure: locating the crossing point of the corresponding susceptibility χ OP and Binder cumulant U OP (measured for different system sizes) for the relevant order parameters, defined as: These phase transitions are associated with the breaking of Z 2 symmetry. Therefore, the critical exponent for the susceptibility near the critical temperature is known, η = 1 4 [33]. In that region, χ OP = L 2−η (f (|1− T Tc |L 1/ν )). Thus the scaled susceptibility is size independent at the critical temperature T c , where different system sizes should show a crossing point. We illustrate this method for the AF 2 in Fig. 4 (b and c) where we show the Binder cumulant (b) and the scaled susceptibility (c) for different system sizes as a function of T /J p (J x /J p = 0.2,J 1 /J p = 1). We can observe that indeed the curves for different L plotted as functions of T /J p for both the Binder cumulant and the normalized susceptibility exhibit a crossing point at the critical temperature, confirming that the transition is of second order. Finally, we remark that the agreement between the Bethe lattice results and the MC-M simulations is both qualitative and quantitative, as can be seen in Fig. 4: the difference in the values of the critical temperatures is 10%. This difference is probably due to both the finite size effect in the MC-M simulations and the fact that the BL analysis is a mean field calculation. We now focus on the J 1 = J x line, where the system exhibits a range of interesting and different phenomena at low temperatures. Our main results are summarized in the J 1 /J p = J x /J p vs T /J p phase diagram in Fig.  5. The transition lines and the tricritical point were obtained with the Bethe lattice technique. The points correspond to the maximum of the specific heat in the MC-M simulations for L = 48. We can identify three types of behavior: (i) a cooperative paramagnet phase for J 1 = J x < J p /3, (ii) a first order phase transition from the paramagnetic phase to the AF 1 phase for J p /3 < J 1 = J x < J * ≈ 0.45J p , and (iii) a second order phase transition from a paramagnetic to a broken Z 2 symmetry phase (AF 1 ) for J 1 = J x > J * . The highly frustrated point J 1 = J x = J p /3 will be discussed in the next subsection. To further characterize these three types of behaviour, we study several variables. Besides the OP AF1 parameter, the specific heat and the entropy, we introduce a variable that we call "J p correlator" (CJ p ) defined in terms of the scalar product of spins connected by J p as, This variable is defined so that CJ p = −1 for both AF 2 and U 2 , and CJ p = +1 in AF 1 . This will be relevant to identify the cooperative paramagnet phase, as we discuss below. We now comment specifically on the physics of each regime of the couplings. To illustrate each case, we show for a specific point in each range of the couplings (J 1 = J x = 0.2J p , 0.4J p , 0.6J p ) the OP AF 1 parameter, the specific heat, the entropy and CJ p as a function of temperature in Fig.6.
(i) For J 1 = J x < J p /3, at low temperatures, the system is in a cooperative paramagnet phase. It is degenerate: the plaquettes are in a mixture of AF 2 and U 2 states. A first indicator of this phase is the behavior of the order parameters and the CJ p correlator at low temperatures. OP AF 1 , OP AF 2 , OP U 2 average up to 0, but the CJ p correlator tends to −1. This indicates that the pairs of spins joined by J p are either +− or −+, and thus that each plaquette is either in an AF 2 or a U 2 state. Another important indicator of the cooperative paramagnet phase is the shape of the specific heat. It can be seen that in this case it shows a broad maximum, there is no sharp feature. This kind of behavior is seen for example in spin ice systems [1,[34][35][36][37][38][39]: it is an indication that no long range order is developed through a thermodynamic phase transition. This cooperative paramagnet phase has extensive entropy. We show the entropy as a function of temperature, obtained from Bethe lattice and MC-WL calculations (Fig.6, third panel from the left). All the curves have the same T → ∞ limit (1 in units of N log 2) as is expected. However, at low temperatures, the red curve shows that for J 1 = J x < J p /3 the system remains disordered as a consequence of plaquette degeneration, and that at low temperatures s = 1 2 . (ii) For J p /3 < J 1 = J x < J * , there is a sharp first order transition to the ordered phase (AF 1 ) at T /J p ∼ 0.55. This is clearly seen as a characteristic jump in the parameters shown in Fig. 6 (blue lines). The specific heat obtained from both the MC-M and MC-WL simulations shows a clear discontinuity (the peak in the simulations is off-scale, and therefore not shown in the figure).
(iii) For J 1 = J x > J * , the system is also ordered at low temperatures: all the plaquettes are in the AF 1 case. The transition from the paramagnetic to the ordered phase is of second order. We have confirmed this by computation of the scaled susceptibility and the Binder cumulant for different system sizes, as done in the previous subsection.
We now exploit the power of the MC-WL technique, studying the Landau free-energy (Eq. (13)) as a function of the relevant order parameter. We will discuss In each region a different behavior is seen. For the larger values of Jx/Jp the system orders in the AF1 phase through a second order transition from the paramagnetic phase. For a smaller range of J1 = Jx > Jp/3, the system orders in the AF1 at low temperatures through a first order transition, clearly seen as a jump in the order parameter and an off-scale peak in the specific heat. For J1 = Jx < Jp/3, the system does not order at low temperature remaining in a macroscopically degenerate state. In the latter case, the specific heat shows a broad peak, as seen for example in spin ice systems. The entropy in this region does not vanish as the temperature is lowered towards T = 0, but tends to the residual value s = 1/2.
how using this variable it is possible to provide further evidence of the different types of phase transitions for J 1 = J x > J p /3. The Landau free-energy as a function of OP AF 1 , for two values of the couplings characteristic of each region, is shown in Fig. 7. On one hand, in Fig. 7 (a) we see that for J x /J p = 0.4 the behavior near the critical temperature T c the position of the global minimum changes abruptly from OP AF 1 min = 0, for T > T c , to OP AF1 min = 1 (normalized) for T < T c . At this point it should be clarified that this shape of the Landau free-energy is characteristic of a finite system size (see for example Ref. [40]), as can be seen in Fig. 7 (c), where the Landau free energy for different system sizes is shown. The curves are flatter with increasing system size L.
In the thermodynamic limit, the curve at T = T c is expected to tend to the dotted curve which has a flat portion between the two minima. On the other hand, typical second-order transition behavior is observed for J x /J p = 0.6 ( Fig. 7 (c)) with a gradual increase of OP AF 1 min from zero, for T > T c to OP AF 1 min ≈ 0 for T T c . This analysis supports previous results obtained with BL and MC-M, and provides a different way to study the nature of the phase transition.
Let us focus on the highly frustrated point, J 1 = J x = J p /3. At this point, the three m 0 = 0 configurations (AF 1 , AF 2 and U 2 ) in each plaquette have the same energy. However, there is a substantial difference between the possible configurations. In the AF 1 case, all plaquettes are in the same arrangement, the degeneracy is only two-fold. In the AF 2 and U 2 case, as was discussed for J 1 = J x < J p /3, the plaquettes can be in any of the two arrangements, thus leading to extensively degenerate ground states. In order to clarify if the system chooses one of these ground states with temperature (and thus thermal order by disorder [41 and 42] is at play), we study the three order parameters (Eqs. (5), (6), (7)), the CJ p correlator (Eq. (15)), the specific heat and the entropy as a function of temperature at h = 0. The results are analogous to those in Fig. 6 (dotted red line) for the case J x /J p < 1/3. This indicates a phenomenon that is not present in the previously discussed case: (partial) thermal order-by-disorder. At low temperatures the system chooses the cooperative paramagnet AF 2 -U 2 phase over the AF 1 one, since the AF 2 -U 2 phase has a lower contribution to the free energy. Now, we center our study in the line J 1(x) = J p /3, J x(1) > J p /3 (horizontal and vertical lines in Fig.  2), where the ground state of the system is either in a AF 1 state or a U 2 (AF 2 ) state. Since the two cases are equivalent, we direct our attention to the J 1 = J p /3, J x > J p /3 case. Contrary to the highly frustrated point, in this case both phases have entropy s = 0 in the thermodynamic limit. However, MC-M simulations show that for lower values of J x the system choses the U 2 phase, whereas for higher values it can be in either state. This is evidence of order-by-disorder state selection at low J x . The origin of this selection can be understood in the following way. In the order-by-disorder phenomenon, in the T → 0 limit, the system chooses states which have a lower contribution to the free energy, even though they are degenerate at T = 0. In this case, the U 2 state has lower energy fluctuations, and thus it has a lower free energy at T → 0. For both U 2 and AF 1 states, the first excitation is simply to flip one spin, and in both cases the energy change is 6J x . To explore the next excitations to the ground states, one can consider the energy difference of flipping spins joined by the different couplings. For the AF 1 state, the next lowest energy excitation is flipping along a J x bond with an energy cost of ∆E (Jx) = 8J x . For the U 2 case for low values of J x flipping along a J p bond has a lower energy cost ∆E (Jp) = 12(J x − J p /3). Therefore, when the temperature is low enough to make these excitations available, the systems chooses the U 2 state. The Bethe lattice technique provides an interesting way of checking this, since the free energy of each state can be calculated. Fig. 7 (d) shows the difference in the free energy between the AF 1 state and the U 2 state as a function of J x for different temperatures. It can be clearly seen that the U 2 state has a lower free energy, but that this difference is smaller at lower temperatures with increasing J x .
In previous sections we have studied the effect of the thermal fluctuations in the stability and transitions to the low temperature phases. The coupling with an external magnetic field can tune the system into high energy plaquette phases inaccessible at h = 0. Here we will compare a typical case with a well defined m = 0 state, for example J 1 = 0.2J p , J x = 0.5J p with the highly frustrated case J 1 = J x = J p /3. The magnetisation and the entropy as a function of temperature and magnetic field by MC-WL simulations are shown in Fig. 8. The m = 0 phases follow the physics previously discussed for h = 0. The transition from the paramagnetic to the m = 1/2 phase is predicted to be of second-order according to the Bethe lattice calculations, which we confirmed with MC-M simulations using the susceptibility and Binder-cumulant analysis as before. A remarkable feature is that the m = 1/2 plateau has finite extensive entropy for both the highly frustrated case and the more typical case presented here, as commented in Sec.II. Finally, high magnetic fields stabilize the saturation plateau m = 1.

V. CONCLUSIONS
In this work, we present a complete study of the Ising model on a bilayer honeycomb lattice including interlayer frustration and coupling with an external magnetic field. We first present and discuss the exact T = 0 phase diagram and we highlight specific points and lines where we expect interesting physics. Then, we study the effect of temperature using a combination of analytical meanfield-like considerations (Bethe lattice) and Monte-Carlo (Metropolis and Wang-Landau) simulations. The interplay between these techniques has been essential to obtain magnetic and thermodynamic properties of the system. We have found a very rich phase diagram with nontrivial regions characterized by broken symmetries and non zero entropy values.
In the case of zero magnetic field, along the highly frustrated line where the intralayer (J 1 ) and frustrating interlayer (J x ) couplings are equal, for J 1 = J x < J p /3 there is a crossover from a paramagnetic phase to a cooperative paramagnet, where there is finite entropy. At the highly frustrated point J 1 = J x = J p /3, order-by-disorder is at play, and this cooperative paramagnet phase is the one selected at low temperatures. For higher values of the couplings, the system is ordered at low temperatures. The Bethe lattice technique shows that this transition is of first order for a certain range of parameters, and then it is second order. We checked this studying thermodynamic variables with both types of simulations. We also used the Wang-Landau technique to study the Landaufree-energy with the corresponding order parameter, to illustrate the difference between a first and a second order phase transitions. Through Monte-Carlo Metropolis simulations and the free energy obtained from the Bethe lattice approximation we showed that order by disorder is also at play in the coexistence lines between two ordered phases (J 1(x) = J p /3, J x(1) > J p /3).
In the presence of a magnetic field, there are three plateaux: at zero and 1/2 magnetisation, and the saturation plateau. The 1/2 plateau is highly degenerate for all antiferromagnetic couplings. This implies that there is a non zero entropy induced by the field through a second order phase transition, even for sets of parameters where in the absence of the external field the system is ordered.
In summary, we have shown that the analytical (Bethe lattice) and numerical (Metropolis and Wang-Landau simulations) techniques are complementary and provide a solid way of exploring the different non trivial phases of this system. We expect to extend this sort of study to other highly degenerate systems, like the highly frustrated honeycomb and kagome lattices in the extended Heisenberg model, exploiting the richness of these techniques to obtain the complete phase diagrams The RR Eqs. (B1) are divergent, and, as usual, we pro-ceed to define new recursion relations that converge in the thermodynamic limit dividing by Z a −− or Z b −− , Eqs. (B1) and (B2) give the following RR where 1.

Fixed points and the thermodynamic phases
The thermodynamic phases are given by the stable fixed points R * of Eqs.(B3). The continuous, or second order lines are defined as coincident stability lines of different fixed points (phases).
At zero magnetic field, the fixed points of the different phases are given by the conditions • Paramagnetic phase • AF 2 phase.
For H = 0 the symmetry is broken, and these relations between the R * i are not valid.

The partition function, thermodynamic averages and the free energy
In order to classify the different thermodynamics phases, we need the partition function and thermody-namics averages, as the magnetisation per site. When the stability lines of two (or more) fixed points are not coincident, the overlap region is a coexistence zone, and we also need to calculate the first-order line as the line where the free energies of the corresponding phases take the same value.
As usual, the thermodynamics averages and the free energy must be calculated on the central region. There is not a unique manner to define the central zone, we define it as a central plaquette where four subtrees are attached as it shown as shown in (Fig. 3).
Putting a plaquette as central zone, always a bond belong to the sublattice a (b), and we attach to them two subtrees with the root belonging to the sublattice b (a), then, at any generations, including the surface, a half of points belong to a sublattice and the other half to the other sublattice. Then, we obtain for the partition function for a M −generations tree, where all the PPF corresponds to M-generations subtrees.