Frustrated honeycomb-lattice bilayer quantum antiferromagnet in a magnetic field: Unconventional phase transitions in a two-dimensional isotropic Heisenberg model
Taras Krokhmalskii Vasyl Baliha Oleg Derzhko Jörg Schulenburg Johannes Richter





Institute for Condensed Matter Physics,
          National Academy of Sciences of Ukraine,
          Svientsitskii Street 1, 79011 L'viv, Ukraine
Department for Theoretical Physics,
          Ivan Franko National University of L'viv,
          Drahomanov Street 12, 79005 L'viv, Ukraine


Institute for Condensed Matter Physics,
          National Academy of Sciences of Ukraine,
          Svientsitskii Street 1, 79011 L'viv, Ukraine


Institute for Condensed Matter Physics,
          National Academy of Sciences of Ukraine,
          Svientsitskii Street 1, 79011 L'viv, Ukraine
Institut für theoretische Physik,
          Otto-von-Guericke-Universität Magdeburg,
          P.O. Box 4120, 39016 Magdeburg, Germany
Department for Theoretical Physics,
          Ivan Franko National University of L'viv,
          Drahomanov Street 12, 79005 L'viv, Ukraine
Abdus Salam International Centre for Theoretical Physics,
          Strada Costiera 11, 34151 Trieste, Italy


Universitätsrechenzentrum,
          Otto-von-Guericke-Universität Magdeburg,
          P.O. Box 4120, 39016 Magdeburg, Germany


Institut für theoretische Physik,
          Otto-von-Guericke-Universität Magdeburg,
          P.O. Box 4120, 39016 Magdeburg, Germany




75.10.-b, 
75.10.Jm  





We consider the spin-1/2 antiferromagnetic Heisenberg model on a bilayer honeycomb lattice including interlayer frustration
in the presence of an external magnetic field.
In the vicinity of the saturation field,
we map the low-energy states of this quantum system onto the spatial configurations of hard hexagons on a honeycomb lattice.
As a result,
we can construct effective classical models (lattice-gas as well as Ising models) on the honeycomb lattice
to calculate the properties of the frustrated quantum Heisenberg spin system in the low-temperature regime.
We perform classical Monte Carlo simulations for a hard-hexagon model and adopt known results for an Ising model
to discuss the finite-temperature order-disorder phase transition
that is driven by a magnetic field at low temperatures.
We also discuss an effective-model description around the ideal frustration case 
and find indications for a spin-flop like transition in the considered isotropic spin model.


Frustrated honeycomb-lattice bilayer quantum antiferromagnet in a magnetic field: 
       Unconventional phase transitions in a two-dimensional isotropic Heisenberg model
    Johannes Richter
    December 30, 2023
==========================================================================================================================================================================





§ INTRODUCTION




An important class of quantum Heisenberg antiferromagnets consists of the so-called two-dimensional dimerized quantum antiferromagnets.
They can be obtained by placing strongly antiferromagnetically interacting pairs of spins 1/2 (dimers) on a regular two-dimensional lattice 
and assuming weak antiferromagnetic interactions between dimers.
Among such models one may mention the J-J^' model
with the staggered arrangement of the strong J^' bonds 
(defining dimers and favoring singlet formation on dimers)
on a square lattice<cit.>
(see also Ref. valenti-eggert for related dimerized square-lattice models).
Other examples are the bilayer models:
They consist of two antiferromagnets in each layer with a dominant nearest-neighbor interlayer coupling which defines dimers.<cit.>
By considering additional frustrating interlayer couplings 
the bilayer model can be pushed in the parameter space to a point 
which admits a rather comprehensive analysis of the energy spectrum.<cit.>
For this special set of coupling parameters, the frustrated bilayer is a system with local conservation laws
(the square of the total spin of each dimer is a good quantum number)
that explains why it is much easier to examine this specific case.
On the other hand,
the frustrated bilayer belongs to the class of so-called localized-magnon spin systems,<cit.>
which exhibit some prominent features around the saturation field, 
such as 
a ground-state magnetization jump at the saturation field,
a finite residual entropy at the saturation field,
and 
an unconventional low-temperature thermodynamics, 
for a review see Refs. review1,review2,review3.
The singlet state of the dimer is the localized-magnon state which belongs to a completely dispersionless (flat) one-magnon band.
Over the last decade  a large  variety of flat-band systems with unconventional physical properties was found, 
see Refs. other1,other2,other3 and references therein.
For the flat-band systems at hand, 
the local nature of the one-magnon states allows to construct also localized many-magnon states 
and to calculate their degeneracy by mapping the problem onto a classical hard-core-object lattice gas;
the case of the frustrated bilayer was discussed in Refs. prb2006,prb2010.
In the strong-field low-temperature regime
the independent localized-magnon states are the lowest-energy ones 
and therefore they dominate the thermodynamics.
The thermodynamic properties in this regime can be efficiently calculated using classical Monte Carlo simulations for a lattice-gas problem.
Even in case of small deviations from the ideal flat-band geometry 
a description which is based on the strong-coupling approach<cit.> can be elaborated.<cit.> 
Again the effective theory is much simpler than that for the initial problem. 

From the theoretical side, frustrated bilayer systems have been studied by several authors.
Thus,
the frustrated square-lattice bilayer quantum Heisenberg antiferromagnet was studied in Refs. lin,prb2006,prb2010,chen,albuquerque,murakami,alet,
whereas the honeycomb-lattice bilayer with frustration was studied in 
Refs. oitmaa,zhang,bishop (intralayer frustration) 
and 
Refs. classical,brenig (interlayer frustration).
For the system to be examined in our paper,
i.e., 
the spin-1/2 antiferromagnetic Heisenberg model on a bilayer honeycomb lattice including interlayer frustration,
H. Zhang et al.<cit.> have determined the quantum phase diagram at zero magnetic field 
for a rather general case of an arbitrary relation between the nearest-neighbor intralayer coupling and the frustrating interlayer coupling.
Another recent study reported in Ref. classical 
concerns the antiferromagnetic classical Heisenberg model on a bilayer honeycomb lattice in a highly frustrated regime 
in the presence of a magnetic field.
Its main result is the phase diagram of the model in the plane “magnetic field – temperature”.
However, this analysis cannot contain any hallmarks caused by the localized magnons, 
since localized-magnon features represent a pure quantum effect which disappears in the classical limit.

From the experimental side,
one may mention several layered materials, which can be viewed as frustrated bilayer quantum Heisenberg antiferromagnets.
Thus, 
the compound Ba_2CoSi_2O_6Cl_2 could be described as a two-dimensionally antiferromagnetically coupled spin-1/2 XY-like spin dimer system
in which Co^2+ sites form the frustrated square-lattice bilayer.<cit.>
The interest in the frustrated honeycomb-lattice bilayers stems from experiments on Bi_3Mn_4O_12(NO_3).<cit.>
In this compound,
the ions Mn^4+ form a frustrated spin-3/2 bilayer honeycomb lattice.<cit.>
Finally, let us mention that a bilayer honeycomb lattice can be realized using ultracold atoms.<cit.>

The present study has several goals.
Motivated by the recent paper of H. Zhang et al.,<cit.> we wish to extend it
to the case of nonzero magnetic field.
On the other hand, with our study we complement the analysis of the classical case<cit.> to the pure quantum case of s=1/2.
Finally, 
the present study can be viewed as an extension of our previous calculations<cit.> to the honeycomb-lattice geometry.
Although we do not intend to provide a theoretical description of Bi_3Mn_4O_12(NO_3),
our results may be relevant for the discussion of the localized-magnon effects in this and in similar materials.

The outline of the paper is as follows.
Section <ref> contains the spectroscopic study of the frustrated honeycomb-lattice bilayer spin-1/2 Heisenberg antiferromagnet:
By exact diagonalization for finite quantum systems and direct calculations for finite hard-core lattice-gas systems
we show the correspondence between the ground states in the large-S^z subspaces and the spatial configurations of hard hexagons on an auxiliary honeycomb lattice.
Based on the established correspondence,
in Section <ref>
we report results of classical Monte Carlo simulations for hard hexagons on the honeycomb lattice
and use them to predict the properties of the frustrated honeycomb-lattice bilayer spin-1/2 Heisenberg antiferromagnet 
in the strong-field low-temperature regime.
The most intriguing outcome is an order-disorder phase transition which is expected at low temperatures just below the saturation
field.
This transition  is related to the ordering of the localized magnons on the two-sublattice honeycomb lattice as the density of the localized magnons increases.
Section <ref> deals with some generalization of the independent localized-magnon picture:
We show how to take into account the contribution of a low-lying set of other localized states
as well as discuss the effect of deviations from the ideal frustration
case.
We end with a summarizing discussion in Section <ref>.
Several technical details are put to the appendixes.



§ INDEPENDENT LOCALIZED-MAGNON STATES








In the present paper,
we consider the spin-1/2 Heisenberg antiferromagnet with the Hamiltonian

    H=∑_⟨ i j⟩J_ijs_i·s_j -hS^z,
       
    J_ij>0,
       
    S^z=∑_i s_i^z

defined on the honeycomb-lattice bilayer shown in Fig. <ref>.
The first sum in Eq. (<ref>) runs over all bonds of the lattice 
and hence J_ij acquires three values: 
J_2 (dimer bonds), 
J_1 (nearest-neighbor intralayer bonds), 
and 
J_ (frustrating interlayer bonds), 
see Fig. <ref>.
In what follows 
we consider  the case J_=J_1 and call it the “ideal frustration case”
(or  “ideal flat-band case”).
Only in Sec. <ref> we discuss deviations from the ideal frustration case, 
i.e., J_ J_1.
Since the z component of the total spin S^z commutes with the Hamiltonian
we can consider the subspaces with different values of S^z separately.

In the strong-field regime the subspaces with large S^z are relevant.
The only state with S^z=N/2 is the fully polarized state
|…↑…⟩
with the energy
E_FM=N(J_2/8+3J_1/4).
In the subspace with S^z=N/2-1 (one-magnon subspace)
N eigenstates of H (<ref>) belong to four one-magnon bands,
E_FM+Λ_k^(α), α=1,2,3,4,
with the dispersion relations:

    Λ_k^(1)=Λ_k^(2)
    =-J_2-3J_1,
       Λ_k^(3,4)
    =-3J_1∓ J_1|γ_k|,
    
    |γ(k)|
    =
    √(3+2[cos k_a +cos k_b +cos(k_a+k_b)]).

Here  
k=(k_x,k_y),
k_a=√(3)a_0 k_x, 
k_b=3a_0k_y/2-√(3)a_0k_x/2,
where a_0 is the hexagon side length,
and k acquires N/2 values from the first Brillouin zone, 
see Appendix A.
The N states from the two flat bands α=1 and α=2 can be chosen as a set of localized states 
where the spin flip is located on one of the N vertical dimers,
see Fig. <ref>.
The remaining N states 
(i.e., from the two dispersive bands α=3 and α=4) 
are extended over the whole lattice.
As can be seen from Eq. (<ref>),
the two-fold degenerate dispersionless (flat) one-magnon band becomes the lowest-energy one,
if J_2> 3J_1,
i.e., if the  strength of the vertical bond J_2 is sufficiently large.
In what follows, 
we assume that this inequality holds.
From the one-magnon spectra (<ref>)
we can also get  the value of the saturation field:
h_sat=J_2+3J_1.

We pass to the many-magnon ground states.
Because the localized one-magnon states have the lowest energy in the one-magnon subspace,
the ground states in the subspaces with S^z=N/2-n, n=2,…, n_max,
n_max=N/2=N/4
can be obtained by populating the dimers.
However, for the ground-state manifold a hard-core constraint is valid,
i.e.,
the neighboring vertical dimers cannot be populated simultaneously, 
since the occupation of nearest neighbors leads to an increase of the energy.
Thus, we arrive at the mapping onto a classical lattice-gas model of hard hexagons on an auxiliary honeycomb lattice:
Each ground state of the quantum spin model can be visualized as a spatial configuration of the hard hexagons on the honeycomb lattice 
excluding the population of neighboring sites (hard-core rule),
see Fig. <ref>.





The occupation of neighboring sites, 
excluded for the ground-state manifold at S^z=N/2-2,…, N/4,
provides another class of localized states 
which can be visualized as overlapping hexagons on the honeycomb lattice:
These states were completely characterized in Refs. andreas and prb2010.
Each overlapping pair of hexagons 
(i.e., occupation of neighboring dimers by localized magnons) 
increases the energy by J_1.
If J_2/J_1 is sufficiently large,
the overlapping hexagon states are the lowest excited states in the subspaces with S^z=N/2-2,…,N/4, 
but they are the ground states in the subspaces with lower S^z.
From exact-diagonalization data for N=24, 32, 36, 48 
we determined the required values of J_2/J_1 as 3.687, 3.781, 3.813, 3.874, respectively:
For these values the first excited state in the subspace with  S^z=N/2-2 are the overlapping hexagon states.

We check our statements on the character of the ground states and the excited states by comparison with exact-diagonalization data.
Clearly,
exact diagonalizations are restricted to finite lattices,
which are shown in Fig. <ref>.
We use the spinpack package  <cit.> 
and exploit the local symmetries to perform numerical exact calculations for large sizes of the Hamiltonian matrix. 
The ground-state degeneracy coincides with the number of spatial configurations of hard hexagons on the honeycomb lattice for all considered cases,
see Table <ref>.
In Table <ref>
we also report the energy gap Δ  to the first excited state and the degeneracy of the first excited state.
While in the one-magnon subspace we have Δ=J_2-3J_1, see Eq. (<ref>),
the energy gap in the subspaces S^z=N/2-n, n=2,…,N/2-1 agrees with the conjecture 
that for large enough J_2/J_1 ≳ 4 the first excited states are other localized-magnon states
for which two of the localized magnons are neighbors
(two hard hexagons overlap), see above.
Further evidence for this picture is provided by the value Δ=2J_1 for S^z=N/4:
The first excited state with respect to the localized-magnon-crystal state corresponds to three overlapping hard hexagons 
resulting in an increase of energy by 2J_1
[see also Eq. (<ref>) in Appendix B].









The zero-temperature magnetization curve is shown by the thick solid red curve in Fig. <ref>.
The magnetization curve probes the ground-state manifold and it is in a perfect agreement with
the above described picture.
There are two characteristic fields,
h_2=J_2
and
h_sat=J_2+3J_1,
at which the ground-state magnetization curve has a jump.
To demonstrate  the robustness of the main features of the magnetization curve against deviations from the ideal frustration case, 
we also show the curve when J_ slightly differs from J_1.
A more detailed discussion of this issue is then provided in Sec. <ref>.





In the next section we use the established correspondence between the spin model and the hard-hexagon model 
to calculate the thermodynamic properties of the frustrated honeycomb-lattice bilayer quantum Heisenberg antiferromagnet
in the strong-field low-temperature regime.



§ HARD HEXAGONS ON THE HONEYCOMB LATTICE




The lowest eigenstates in the subspaces with large S^z become ground states for strong magnetic fields.
Thus, 
the energy of these lowest eigenstates in the subspaces with S^z=N/2-n, n=0,1,…,n_max in the presence of the field h
is

    E^lm_n(h)=E_FM-hN/2-(ϵ_1-h)n,
     ϵ_1=J_2+3J_1 .

At the saturation field, 
i.e., at h=h_sat=ϵ_1, all these energies 
become independent of n,
E^lm_n(h_sat)=E_FM-ϵ_1 N/2.
Therefore, the system exhibits  a huge ground-state degeneracy at
h_sat which grows exponentially with the system size N:<cit.>
W=∑_n=0^N/2g_N(n)≈exp(0.218 N), see Eq. (<ref>) below.
Here g_N(n) denotes the degeneracy of the ground state for the 2N-site frustrated honeycomb-lattice bilayer in the subspace with S^z=N/2-n.

Furthermore,
following Refs. mike and review2,
the contribution of the independent localized-magnon states to the partition function is given by the following formula:

    Z_lm(T,h,N)=∑_n=0^N/2g_N(n)exp[-E^lm_n(h)/T].

Since g_N(n)=Z_hc(n,N) is the canonical partition function of n hard hexagons on the N-site honeycomb lattice,
Eq. (<ref>) can be rewritten as

    Z_lm(T,h,N)=exp(-E_FM-hN/2/T) Ξ_hc(T,μ,N),
    
    Ξ_hc(T,μ,N)=∑_n=0^N/2Z_hc(n,N)exp(μ n/T),
    μ=ϵ_1-h.

As a result, 
we get the following relations:

    F_lm(T,h,N)/N
    =
    E_FM/N-h/2+1/2Ω_hc (T,μ,N)/N,
    
    Ω_hc (T,μ,N)
    =
    -TlnΞ_hc(T,μ,N)

for the free energy per site f(T,h),

    M_lm(T,h,N)/N
    =
    1/2+1/2∂/∂μΩ_hc(T,μ,N)/N

for the magnetization per site m(T,h),

    S_lm(T,h,N)/N
    =
    1/2S_hc(T,μ,N)/N

for the entropy per site s(T,h),

    C_lm(T,h,N)/N
    =
    1/2C_hc(T,μ,N)/N

for the specific heat per site c(T,h).
Note that h and μ are related by μ=h_sat-h.
The hard-hexagon quantities in the r.h.s. of these equations depend on the temperature and
the chemical potential only through the activity z=exp(μ/T).
That means that for the frustrated quantum spin system at hand all thermodynamic
quantities depend on temperature and magnetic field only via x=(h_sat-h)/T=ln z,
i.e., a universal behavior emerges in this regime.

To check the formulas for thermodynamic quantities  given in Eqs. (<ref>) – (<ref>)
we compare the exact-diagonalization data with the predictions based on the hard-hexagon picture.
We set J_1=1, J_2=5 and perform exact-diagonalization calculations 
for thermodynamics for the frustrated quantum spin system of N=24 sites,<cit.> see Fig. <ref>, 
where the total size of the Hamiltonian matrix is already 
16 777 216 × 16 777 216.
We also perform the simpler calculations for the corresponding hard-hexagon systems, 
see Appendix B.

Our results for temperature dependences of the specific heat around the saturation are collected in Figs. <ref> and <ref>.
As can be seen from these plots,
the hard-hexagon picture perfectly reproduces the low-temperature features of the frustrated quantum spin model around the saturation field.
Deviations from the  hard-core-model predictions in the upper panel of Fig. <ref> become visible only at T=0.2.
From the middle panel of Fig. <ref> one can conclude 
that the temperature profiles for specific heat at h=7.95 and h=8.05 are well described by the hard-core model again up to about T=0.1.









Using the correspondence between the frustrated quantum spin model and the classical hard-core-object lattice-gas model,
we can give a number of predictions for the former model based on the analysis of the latter one.
For example,
we can calculate the ground-state entropy at the saturation field:

    S(T→ 0,h=h_sat,N)/2N
    =lnΞ_hc(z=1,N)/2N≈ 0.218.

This number follows by direct calculations for finite lattices up to N=64 sites.
On the other hand,
for the problem of hard hexagons on a honeycomb lattice
κ=exp[lnΞ_hc(z=1,N)/N]= 1.546…
plays 
the same role as 
the hard-square entropy constant
κ=1.503048082…
for hard squares on the square lattice
or 
the hard-hexagon entropy constant 
κ=1.395485972…
for hard hexagons on the triangular lattice.<cit.>
Such constants determine the asymptotic growth and are also of interest to combinatorialists.
A more precise value of this constant for hard hexagons on a honeycomb lattice can be found in Ref. baxter2.

The most interesting consequence of the correspondence between the frustrated quantum bilayer and the hard-core lattice gas is the existence of an order-disorder phase transition.
It is generally known
that for the lattice-gas model on the honeycomb lattice with first neighbor exclusion
the  
hard hexagons spontaneously occupy one of two sublattices of the honeycomb lattice
as the activity z exceeds the critical value z_c=7.92…, see Ref. baxter2.
In the spin language,
this corresponds to the ordering of the localized magnons as their density increases. 
This occurs at low temperatures just below the saturation field.
For the fixed (small) deviation from the saturation field, 
h_sat-h,
the formula for the critical temperature reads:

    T_c=h_sat-h/ln z_c≈ 0.48(h_sat-h).

Furthermore,
the critical behavior falls into the universality class of the two-dimensional-Ising-model.<cit.>
That means, the specific heat at T_c (<ref>) shows a logarithmic singularity.
Of course, the calculated T_c (<ref>) must be small,
otherwise the elaborated effective low-energy theory fails,
see Figs. <ref> and <ref>.



§ BEYOND INDEPENDENT LOCALIZED-MAGNON STATES






 §.§ Other localized-magnon states



Following Ref. prb2010,
in addition to the independent localized-magnon states 
(which obey the hard-hexagon rule) 
we may also take into account another class of localized-magnon states which correspond to overlapping hexagon states 
(i.e., they violate the hard-hexagon rule), 
see also our discussion in Sec. <ref>.
The corresponding lattice-gas Hamiltonian has the form:

    H({n_m})
    =
    -μ∑_m=1^Nn_m
    +V∑_⟨ mn⟩ n_mn_n.

Here 
n_m=0,1 is the occupation number attached to each site m=1,…,N of the auxiliary honeycomb lattice,
the first (second) sum runs over all sites (nearest-neighbor bonds) of this auxiliary lattice,
and μ=h_sat-h, V=J_1.
The interaction describes the energy increase if two neighboring sites are occupied by hexagons.
In the limit V→∞ the hard-core rule is restored.

The partition function is given by

    Z_LM(T,h,N)
    =
    exp(-E_FM-hN/2/T) Ξ_lg(T,μ,N),
    
    Ξ_lg(T,μ,N)
    =∑_n_1=0,1…∑_n_N=0,1exp[-H({n_m})/T].

Since  Z_LM contains not only the contribution from independent localized-magnon states, 
but also from  overlapping localized-magnon states, 
it is valid in a significantly wider region of magnetic fields and temperatures.

Evidently,
new Ising variables σ_m=2n_m-1 may be introduced in Eqs. (<ref>) and (<ref>)
and as a result we face the antiferromagnetic honeycomb-lattice Ising model in a uniform magnetic field:

    H
    =
    N(-μ/2+3V/8)
    -Γ∑_m=1^Nσ_m
    +J∑_⟨ mn⟩σ_mσ_n,
    
    Γ=μ/2-3V/4,
       J=V/4>0.

The Ising variable σ_m acquires two values ± 1,
the nearest-neighbor interaction J=J_1/4>0 is antiferromagnetic,
and the effective magnetic field Γ=(h_sat-h)/2 -3J_1/4=(J_2+3J_1/2-h)/2 is zero when h=J_2+3J_1/2.
The zero-field case (i.e., Γ=0) is exactly solvable, 
see Ref. jozef-michal and references therein.
For example, 
the critical temperature is known to be
T_c/J_1=1/[2ln (2+√(3))]≈ 0.380.
The ground-state antiferromagnetic order in the model (<ref>) survives at T=0
at small fields 
|Γ|<3J, i.e.,
for h_2<h<h_sat, h_2=J_2, h_sat=J_2+3J_1.
The antiferromagnetic honeycomb-lattice Ising model in a uniform magnetic field
was a subject of several studies in the past.<cit.>
In particular,
several closed-form expressions for the critical line in the plane “magnetic field – temperature”
which are in good agreement with numerical results were obtained,
see Refs. honeycomb-ising1,honeycomb-ising2 and also Refs. honeycomb-ising3,honeycomb-ising4.
On the basis of these studies we can construct the phase diagram, see Fig. <ref>.
Here we have used the two closed-form expressions for the critical line of the antiferromagnetic Ising model in a magnetic field  
suggested in Refs. honeycomb-ising1 and honeycomb-ising2, 
where both are indistinguishable in the scale used in Fig. <ref>.
Although the two-dimensional Ising model in a field has not been solved
analytically,
the results of Refs. honeycomb-ising1,honeycomb-ising2 are known to be very accurate.<cit.>





In Fig. <ref> 
the temperature profiles for the specific heat in a wide range of magnetic fields
are shown.
The comparison with the exact-diagonalization data demonstrates a clear improvement of the hard-hexagon description
after using the lattice-gas model (<ref>).
Furthermore,
in the lower panel of Fig. <ref> we report classical Monte Carlo data for h=5.05 
[lattice-gas model (<ref>)]
which shows how the temperature profile C(T,h,N)/N modifies and develops a singularity
as the system size increases
[see C(T,h,N)/N for N=288^2 in the lower panel of Fig. <ref>].



 §.§ Deviation from the ideal flat-band geometry



Following Ref. around,
we can consider an effective low-energy description
when the flat-band conditions are slightly violated and the former flat band acquires a small dispersion.
To this end,
we assume that the intralayer nearest-neighbor interaction J_1 and the interlayer frustrating interaction J_ are different,
but the difference is small | J_1-J_|/J_2≪ 1.
Then in the strong-field low-temperature regime there are two relevant states at each dimer:
| u⟩=|↑↑⟩ 
and 
| d⟩=(|↑↓⟩ - |↓↑⟩)/√(2).
Their energies, 
ϵ_u=J_2/4-h
and
ϵ_d=-3J_2/4,
coincide at h=h_0=J_2.
Now the 2^N-fold degenerate ground-state manifold is splitted by the perturbation,
which consists of the Zeeman term -(h-h_0)∑_i s_i^z
and the interdimer interactions with the coupling constants J_1 and J_.
The effective Hamiltonian acting in the ground-state manifold can be found perturbatively:<cit.>

    H_eff=PHP+…,

where P=|φ_0⟩⟨φ_0| is the projector onto the ground-state manifold,
|φ_0⟩=∏_m=1^N|v⟩,
where |v⟩ is either the state | u⟩ or the state | d⟩.
After some straightforward calculations and introducing the (pseudo)spin-1/2 operators
T^z=(| u⟩⟨ u| - | d⟩⟨ d|)/2,
T^+= | u⟩⟨ d|,
T^-= | d⟩⟨ u|
at each vertical bond
we arrive at the following result:

    H_eff
    =
    N(-h/2-J_2/4+3J/8)
    -h∑_m=1^NT_m^z
    
    
    +∑_⟨ mn⟩[J^z T_m^zT_n^z+J(T_m^xT_n^x+T_m^yT_n^y)],
    
    h=h-J_2-3J/2,
    J=J_1+J_/2,
    J^z=J,
    J=J_1-J_.

The second sum in Eq. (<ref>) runs over all 3N/2 nearest-neighbor bonds of the auxiliary honeycomb lattice.
Note that the sign of the coupling constant J is not important, since the
auxiliary-lattice model (<ref>) is bipartite.
Again the effective Hamiltonian (<ref>) 
which corresponds to the spin-1/2 XXZ Heisenberg model in a z-aligned field on the honeycomb lattice
is much simpler than the initial model and it can be studied further by, e.g., the
quantum Monte Carlo method.<cit.>

For the ideal flat-band geometry (ideal frustration case) the effective Hamiltonian (<ref>) transforms into the above discussed lattice-gas or Ising models.
To make this evident we have to take into account that 
J=J_1=V, 
h=h-h_sat+3J/2=-μ+3V/2, 
J^z=J_1=V, 
J=0, 
and replace T^z by -σ/2:

    H_eff
    =
    N(-h/2-J_2/4+3V/8)
    
    
    -(μ/2-3V/4)∑_m=1^Nσ_m
    +V/4∑_⟨ mn⟩σ_mσ_n
    
    
    =E_FM-hN
    +N(-μ/2 +3V/8)
    
    
    -(μ/2-3V/4)∑_m=1^Nσ_m
    +V/4∑_⟨ mn⟩σ_mσ_n,

cf. Eqs. (<ref>) and (<ref>).





To illustrate the quality of the effective description,
we compare the results for the ground-state magnetization curve obtained by exact diagonalization
for the initial model of N=32 sites
(thin solid curves in Fig. <ref>)
and
for the effective model of N=16 sites
(thin dashed curves in Fig. <ref>).

It is worth noting the symmetry present in the Hamiltonian (<ref>):
If one replaces 
h=J_2+3J+δ h to h=J_2- δ h
and
all T_m^z to -T_m^z
the Hamiltonian (<ref>) (up to the constant) remains the same.
This symmetry of the effective model is also present in the exact-diagonalization 
data for the initial model,
if deviations from the flat-band geometry are small,
see the thin solid black curve and the thin solid brown curve in Figs. <ref> and <ref>.
Moreover, it is also obvious in the lattice-gas Hamiltonian (<ref>):
After the replacement 
μ=δμ to μ=3J_1-δμ 
and
all n_m to 1-n_m
the Hamiltonian (<ref>) (up to the constant) remains the same.

As can be seen in Fig. <ref>,
the magnetization jumps survive even for moderate deviations from the ideal frustration
case.
The nature of the jump is evident from the effective model (<ref>):
It is a spin-flop transition, 
which is present in a two-dimensional Ising-like XXZ Heisenberg antiferromagnet in an external field along the easy axis,
see, e.g., Refs. spin-flop-yunoki,spin-flop-selke,spin-flop.
Note that according to Eq. (<ref>)  
the effective easy-axis XXZ model becomes isotropic for J_1+J_=2(J_1-J_), 
i.e., the spin-flop transition disappears as increasing the deviation from the ideal frustration case J_1=J_.
Although, we are not aware of previous studies of the spin-flop transition for 
the honeycomb-lattice spin-1/2 XXZ model
(and such a study is beyond the scope of the present paper),
we may mention here that the square-lattice case was examined in Ref. spin-flop-yunoki.
In particular, 
one may find there the dependences of the height of the magnetization jump and of the
transition field on the anisotropy.
Furthermore, for temperature effects, see Ref. spin-flop-selke.
Supposing, that for the honeycomb-lattice case the same scenario as for the square-lattice case is valid, 
we may expect that the spin-flop like transition in our model only disappears at the isotropic point J_1+J_=2(J_1-J_). 
Our quantum Monte Carlo data shown in the inset of Fig. <ref> support this conclusion.

Let us complete this section with a general remark on effective models around the ideal flat-band geometry (the ideal frustration case).
Recalling the findings of Ref. around, 
where several localized-magnon systems including the square-kagome model were examined,
we conclude that the effective model around the ideal flat-band geometry essentially depends on the universality class of the localized-magnon system.
(For a comprehensive discussion of the various universality classes of localized-magnon
systems, see Ref. univ_class.)
While for the square-kagome model falling into the monomer universality class we
obtained the (pseudo)spin-1/2 XXZ models with easy-plane anisotropy,<cit.>
for the considered frustrated honeycomb-lattice bilayer model, 
which belongs to a hard-hexagon universality class,
we get the (pseudo)spin-1/2 XXZ models with easy-axis anisotropy.
Clearly, the magnitude of the Ising terms in the effective Hamiltonian are related to the specific hard-core rules.



§ CONCLUSIONS




In this paper
we examine 
the low-temperature properties of the frustrated honeycomb-lattice bilayer spin-1/2 Heisenberg antiferromagnet in a magnetic field.
For the considered model,
when the system has local conservation laws,
it is possible to construct a subset of 2^N eigenstates
(N=N/2)  
of the Hamiltonian and to calculate their contribution to thermodynamics.
For sufficiently strong interlayer coupling, 
these states are low-energy ones for strong and intermediate fields 
and therefore they dominate the thermodynamic properties.

The most interesting features of the studied frustrated quantum spin model are:
The magnetization jumps as well as wide plateaus,
the residual ground-state entropy,
the extra low-temperature peak in the temperature dependence of the specific heat around the saturation,
and the finite-temperature order-disorder phase transition of the two-dimensional Ising-model universality class.
The phase transition occurs just below the saturation field h_sat.
However, for large enough J_2/J_1,  
there is a line of phase transitions which occur below T_c/J_1=1/[2ln (2+√(3))]≈ 0.380
for h in the region between h_2=J_2 and h_sat=J_2+3J_1.
Finally,
for deviations from the ideal frustration case we observe for the
isotropic Heisenberg model at hand magnetization jumps which can be understood
as spin-flop like transitions.

There might be some relevance of our study for the magnetic compound
Bi_3Mn_4O_12(NO_3).
The most intriguing question is:
Can the phase diagram from Fig. <ref> be observed experimentally?
First, the exchange couplings for Bi_3Mn_4O_12(NO_3) are still under debate<cit.> but the relation J_2/J_1≈ 2 looks plausible.
In this case the flat band is not the lowest-energy one, see Eq. (<ref>). 
Second, the spin value is s=3/2 for this compound 
(each Mn^4+ ion carries a spin s=3/2)
and the localized-magnon effects are less pronounced in comparison with the s=1/2
case.
For example, the magnitude ground-state magnetization jump at the saturation is still
N/2, but this magnitude is only 1/6 of the saturation value
(in contrast to 1/2 of the saturation value for the s=1/2 case). 
Thus, further studies on this compound are needed 
to clarify the relation to the localized-magnon scenario presented in our paper.



§ ACKNOWLEDGMENTS


The present study was supported by the Deutsche Forschungsgemeinschaft (project RI615/21-2).
O. D. acknowledges the kind hospitality of the University of Magdeburg in October-December of 2016.
The work of T. K. and O. D. was partially supported by Project FF-30F (No. 0116U001539) from the Ministry of Education and Science of Ukraine.
O. D. would like to thank the Abdus Salam International Centre for Theoretical Physics (Trieste, Italy) 
for partial support of these studies through the Senior Associate award.

00

jjprime
R. R. P. Singh, M. P. Gelfand, and D. A. Huse,
Phys. Rev. Lett. 61, 2484 (1988);
S. E. Krüger, J. Richter, J. Schulenburg, D. J. J. Farnell, and R. F. Bishop,
Phys. Rev. B 61, 14607 (2000);
S. Wenzel, L. Bogacz, and W. Janke,
Phys. Rev. Lett. 101, 127202 (2008);
L. Fritz, R. L. Doretto, S. Wessel, S. Wenzel, S. Burdin, and M. Vojta,
Phys. Rev. B 83, 174416 (2011);
J. Richter and O. Derzhko,
to appear in Eur. J. Phys.
[arXiv:1611.09655].

valenti-eggert
U. Tutsch, B. Wolf, S. Wessel, L. Postulka, Y. Tsui, H. O. Jeschke, I. Opahle, T. Saha-Dasgupta, R. Valenti, A. Brühl,
K. Remović-Langer, T. Kretz, H.-W. Lerner, M. Wagner, and M. Lang,
Nature Communications 5, 5169 (2014);
D. Straßel, P. Kopietz, and S. Eggert,
Phys. Rev. B 91, 134406 (2015).

bilayer
K. Hida,
J. Phys. Soc. Jpn. 59, 2230 (1990);
L. Wang, K. S. D. Beach, and A. W. Sandvik,
Phys. Rev. B 73, 014431 (2006);
R. Ganesh, D. N. Sheng, Y.-J. Kim, and A. Paramekanti,
Phys. Rev. B 83, 144414 (2011);
R. Ganesh, S. V. Isakov, and A. Paramekanti,
Phys. Rev. B 84, 214412 (2011).

andreas
A. Honecker, F. Mila, and M. Troyer,
Eur. Phys. J. B 15, 227 (2000). 

prl2002
J. Schulenburg, A. Honecker, J. Schnack, J. Richter, and H.-J. Schmidt,
Phys. Rev. Lett. 88, 167207 (2002).

review1
J. Richter,
Fizika Nizkikh Temperatur (Kharkiv) 31, 918 (2005)
[Low Temperature Physics 31, 695 (2005)].

review2
O. Derzhko, J. Richter, A. Honecker, and H.-J. Schmidt,
Fizika Nizkikh Temperatur (Kharkiv) 33, 982 (2007) 
[Low Temperature Physics 33, 745 (2007)].

review3
O. Derzhko, J. Richter, and M. Maksymenko,
Int. J. Mod. Phys. B 29, 1530007 (2015).

other1
A. Mielke, J. Phys. A 24, L73 (1991); 
A. Mielke, J. Phys. A 24, 3311 (1991); 
A. Mielke, J. Phys. A 25, 4335 (1992); 
A. Mielke, Phys. Lett. A 174, 443 (1993);
H. Tasaki, Phys. Rev. Lett. 69, 1608 (1992); 
A. Mielke and H. Tasaki, Commun. Math. Phys. 158, 341 (1993);
H. Tasaki, J. Phys.: Condens. Matter 10, 4353 (1998); 
H. Tasaki, Prog. Theor. Phys. 99, 489 (1998).

other2
D. Leykam, S. Flach, O. Bahat-Treidel, and A. S. Desyatnikov,
Phys. Rev. B 88, 224203 (2013);
S. Flach, D. Leykam, J. D. Bodyfelt, P. Matthies, and A. S. Desyatnikov,
EPL 105, 10001 (2014);
W. Maimaiti, A. Andreanov, H. C. Park, O. Gendelman, and S. Flach,
arXiv:1610.02970.

other3
L. Morales-Inostroza and R. A. Vicencio,
Phys. Rev. A 94, 043831 (2016).

prb2006
J. Richter, O. Derzhko, and T. Krokhmalskii,
Phys. Rev. B 74, 144430 (2006). 

prb2010
O. Derzhko, T. Krokhmalskii, and J. Richter,
Phys. Rev. B 82, 214412 (2010);
O. Derzhko, T. Krokhmalskii, and J. Richter,
Teoret. Mat. Fiz. 168, 441 (2011)
[Theoretical and Mathematical Physics 168, 1236 (2011)].

aa
A. Honecker and A. Läuchli,
Phys. Rev. B 63, 174407 (2001).

around
O. Derzhko, J. Richter, O. Krupnitska, and T. Krokhmalskii,
Phys. Rev. B 88, 094426 (2013);
J. Richter, O. Krupnitska, T. Krokhmalskii, and O. Derzhko,
J. Magn. Magn. Mater. 379, 39 (2015).

lin 
H. Q. Lin and J. L. Shen, 
J. Phys. Soc. Jpn. 69, 878 (2000).

chen
P. Chen, C. Y. Lai, and M. F. Yang, 
Phys. Rev. B 81, 020409(R) (2010).

albuquerque
A. F. Albuquerque, N. Laflorencie, J. D. Picon, and F. Mila, 
Phys. Rev. B 83, 174421 (2011).

murakami
Y. Murakami, T. Oka, and H. Aoki,
Phys. Rev. B 88, 224404 (2013).

alet
F. Alet, K. Damle, and S. Pujari,
Phys. Rev. Lett. 117, 197203 (2016).

oitmaa
J. Oitmaa and R. R. P. Singh,
Phys. Rev. B 85, 014428 (2012).

zhang
H. Zhang, M. Arlego, and C. A. Lamas,
Phys. Rev. B 89, 024403 (2014).

bishop
R. F. Bishop and P. H. Y. Li,
arXiv:1611.03287.
 
classical
F. A. Gómez Albarracin and H. D. Rosales,
Phys. Rev. B 93, 144413 (2016).

brenig
H. Zhang, C. A. Lamas, M. Arlego, and W. Brenig,
Phys. Rev. B 93, 235150 (2016).

tanaka
H. Tanaka, N. Kurita, M. Okada, E. Kunihiro, Y. Shirata, K. Fujii, H. Uekusa, A. Matsuo, K. Kindo, and H. Nojiri,
J. Phys. Soc. Jpn. 83, 103701 (2014).

smirnova
O. Smirnova, M. Azuma, N. Kumada, Y. Kusano, M. Matsuda, Y. Shimakawa, T. Takei, Y. Yonesaki, and N. Kinomura, 
J. Am. Chem. Soc. 131, 8313 (2009);
S. Okubo, F. Elmasry, W. Zhang, M. Fujisawa, T. Sakurai, H. Ohta, M. Azuma, O. A. Sumirnova, and N. Kumada, 
J. Phys.: Conf. Ser. 200, 022042 (2010);
M. Matsuda, M. Azuma, M. Tokunaga, Y. Shimakawa, and N. Kumada,
Phys. Rev. Lett. 105, 187201 (2010).

kandpal
H. C. Kandpal and J. van den Brink,
Phys. Rev. B 83, 140412(R) (2011).

cold_atoms
S. Dey and R. Sensarma,
Phys. Rev. B 94, 235107 (2016).

spin
J. Richter and J. Schulenburg,
Eur. Phys. J. B 73, 117 (2010);
https://www-e.uni-magdeburg.de/jschulen/spin/

prb2004
O. Derzhko and J. Richter,
Phys. Rev. B 70, 104415 (2004).

mike
M. E. Zhitomirsky and H. Tsunetsugu,
Phys. Rev. B 70, 100403(R) (2004);
M. E. Zhitomirsky and H. Tsunetsugu,
Prog. Theor. Phys. Suppl. 160, 361 (2005).

baxter
R. J. Baxter, I. G. Enting, and S. K. Tsang,   
J. Stat. Phys. 22, 465 (1980);
R. J. Baxter and  S. K. Tsang,  
J. Phys. A 13, 1023 (1980);
R. J. Baxter, 
Exactly Solved Models in Statistical Mechanics
(Academic Press, New York, 1982). 

baxter2
R. J. Baxter,
Ann. Comb. 3, 191 (1999)
[arXiv:cond-mat/9811264].

debierre
J.-M. Debierre and L. Turban,
Phys. Lett. A 97, 235 (1983);
L. Turban and J.-M. Debierre,
Phys. Lett. A 103, 81 (1984).

jozef-michal
J. Strečka and M. Jaščur,
Acta Physica Slovaca 65, 235 (2015).

honeycomb-ising1
F. Y. Wu, X. N. Wu, and H. W. J. Blöte,
Phys. Rev. Lett. 62, 2773 (1989); 63, 696 (1989).

honeycomb-ising2
X.-Z. Wang and J. S. Kim,
Phys. Rev. Lett. 78, 413 (1997);
X.-Z. Wang and J. S. Kim,
Phys. Rev. E 56, 2793 (1997).

honeycomb-ising3
S.-Y. Kim,
Phys. Lett. A 358, 245 (2006).

honeycomb-ising4
S. L. A. de Queiroz,
Phys. Rev. E 87, 024102 (2013).

fulde
P. Fulde,
Electron Correlations in Molecules and Solids
(Springer-Verlag, Berlin, Heidelberg, 1993),
p. 77.

mila
F. Mila and K. P. Schmidt,
in 
Introduction to Frustrated Magnetism,
Springer Series in Solid-State Sciences 164,
edited by C. Lacroix, P. Mendels, and F. Mila
(Springer-Verlag, Berlin, Heidelberg, 2011),
pp. 537-559.

alps
A. F. Albuquerque et al. (ALPS collaboration), 
J. Magn. Magn. Mater. 310, 1187 (2007); 
B. Bauer et al. (ALPS collaboration), 
J. Stat. Mech. P05001 (2011).

spin-flop-yunoki
S. Yunoki,
Phys. Rev. B 65, 092402 (2002).

spin-flop-selke
M. Holtschneider, S. Wessel, and W. Selke,
Phys. Rev. B 75, 224417 (2007).

spin-flop
K. Balamurugan, S.-H. Lee, J.-S. Kim, J.-M. Ok, Y.-J. Jo, Y.-M. Song, S.-A. Kim, E. S. Choi, M. D. Le, and J.-G. Park,
Phys. Rev. B 90, 104412 (2014).

univ_class 
O. Derzhko and  J. Richter,
Eur. Phys. J. B 52, 23 (2006). 







§ APPENDIX A: ONE-MAGNON ENERGIES (<REF>)




In this appendix,
we present the calculation of the one-magnon energies (<ref>).
In the one-magnon subspace,
the Hamiltonian (<ref>) can be written in the following form (see Fig. <ref>):

    H=∑_m_a=0^L-1∑_m_b=0^L-1(
    J_2 h_m_a,m_b,1;m_a,m_b,3 +J_2 h_m_a,m_b,2;m_a,m_b,4.
    
    .
    +J_1 h_m_a,m_b,1;m_a,m_b,2
    +J_1 h_m_a,m_b,1;m_a,m_b,4
    +J_1 h_m_a,m_b,3;m_a,m_b,4
    +J_1 h_m_a,m_b,3;m_a,m_b,2.
    
    .
    +J_1 h_m_a,m_b,2;m_a,m_b+1,1
    +J_1 h_m_a,m_b,2;m_a,m_b+1,3
    +J_1 h_m_a,m_b,4;m_a,m_b+1,3
    +J_1 h_m_a,m_b,4;m_a,m_b+1,1.
    
    .
    +J_1 h_m_a,m_b,2;m_a+1,m_b+1,1
    +J_1 h_m_a,m_b,2;m_a+1,m_b+1,3
    +J_1 h_m_a,m_b,4;m_a+1,m_b+1,3
    +J_1 h_m_a,m_b,4;m_a+1,m_b+1,1),
    
    
    h_i;j
    =
    1/2(s_i^-s_j^+ + s_j^-s_i^+)
    -1/2(s_i^-s_i^+ + s_j^-s_j^+)
    +1/4.

Recall that N is the number of sites,
N=N/2 is the number of vertical bonds,
and N/2=L^2 is the number of sites of the triangular lattice which is used here:
The honeycomb-lattice bilayer is viewed as a triangular lattice with four sites in the unit cell.

We introduce the Fourier transformation,

    s^+_m_a,m_b,α
    =
    1/L∑_k_a∑_k_bexp[i(k_am_a+k_bm_b)] s^+_k,α,
    
    
    s^-_m_a,m_b,α
    =
    1/L∑_k_a∑_k_bexp[-i(k_am_a+k_bm_b)] s^-_k,α,
    
    k=k_a2/3a_0(√(3)/2i+1/2j)+k_b2/3a_0j,
    
    
    k_a=2π/Lz_a, z_a=0,1,…,L -1,
    k_b=2π/Lz_b, z_b=0,1,…,L -1,

a_0 is the hexagon side length.
After that, Hamiltonian (<ref>) can be cast into

    H
    =
    N/2(J_2/2+3J_1)
    +
    ∑_k(
    [ s^-_k,1 s^-_k,2 s^-_k,3 s^-_k,4 ])
    (
    [ H_11 H_12 H_13 H_14; H_21 H_22 H_23 H_24; H_31 H_32 H_33 H_34; H_41 H_42 H_43 H_44 ])
    (
    [ s^+_k,1; s^+_k,2; s^+_k,3; s^+_k,4 ])

with the following matrix H

    H
    =
    (
    [ -J_2/2-3J_1    J_1/2γ_k       J_2/2    J_1/2γ_k;  J_1/2γ^*_k -J_2/2-3J_1  J_1/2γ^*_k       J_2/2;       J_2/2    J_1/2γ_k -J_2/2-3J_1    J_1/2γ_k;  J_1/2γ^*_k       J_2/2  J_1/2γ^*_k -J_2/2-3J_1 ]),
    
    γ_k=1+exp(-ik_b)+exp[-i(k_a+k_b)].

The eigenvalues of the matrix H are as follows:

    { 
    -J_2-3J_1,
       
    -J_2-3J_1,
       
    -3J_1-J_1|γ_k|,
       
    -3J_1+J_1|γ_k|},
    
    |γ_k|
    =√(3+2[cos k_a +cos k_b +cos(k_a+k_b)]).

Therefore,
in the one-magnon subspace we have

    H
    =
    E_FM
    +
    ∑_k∑_α=1,2,3,4Λ_k^(α)s_k,α^-s_k,α^+,
    
    Λ_k^(1)=Λ_k^(2)
    =-J_2-3J_1,
       Λ_k^(3,4)
    =-3J_1∓ J_1|γ_k|.




§ APPENDIX B: FINITE LATTICES




In this appendix,
we collect some formulas for finite lattices.

Consider hard hexagons on the (periodic) N=12 site honeycomb lattice,
see Fig. <ref>.
Then

    Ξ_hc(T,μ, 12)
    =1 + 12z+ 48z^2+ 76z^3 +45z^4 +12z^5+ 2z^6,
    
    
    z=exp(μ/T),
       μ=h_sat-h,

see Table <ref>.
All thermodynamic quantities follow from Eq. (<ref>) according to standard prescriptions of statistical mechanics.

Next,
consider the lattice-gas model with finite nearest-neighbor repulsion on the (periodic) N=12 site honeycomb lattice,
see Fig. <ref>.
Then

    Ξ_lg(T,μ, 12)
    =
    ∑_n_1=0,1…∑_n_12=0,1exp[
    μ/T∑_m=1^12 n_m 
    .
    
    .
    -V/T(
    n_5n_1+n_1n_2+n_2n_3+n_4n_2+n_3n_5+n_9n_4+n_4n_6+n_5n_7+n_6n_7
    .
    .
    
    .
    .
    +n_8n_6+n_7n_9+n_12n_8+n_8n_10+n_9n_11+n_10n_11+n_1n_10+n_11n_12+n_12n_3
    )
    ]
    
    
    =
    ∑_n_1=0,1…∑_n_12=0,1
    z^∑_m=1^12 n_mexp[
    -V/T(
    n_5n_1+n_1n_2+n_2n_3+n_4n_2+n_3n_5+n_9n_4+n_4n_6+n_5n_7+n_6n_7
    .
    .
    
    .
    .
    +n_8n_6+n_7n_9+n_12n_8+n_8n_10+n_9n_11+n_10n_11+n_1n_10+n_11n_12+n_12n_3
    )
    ],
    
    μ=h_sat-h,
       
    V=J_1,
       
    z=exp(μ/T).

The partition function (<ref>) contains 4096 terms and can be easily calculated.
All thermodynamic quantities follow from Eq. (<ref>) according to standard prescriptions of statistical mechanics.

Clearly, Eq. (<ref>) implies a specific numbering of sites in Fig. <ref>.
However, it can be rewritten in the form that does not depend on the site numbering
[as Eq. (<ref>)].
If we introduce the function g(k_1,k_2) with 
the integer k_1=0,…, 12 to count the number of occupied sites
and 
the integer k_2=0,…,18 to count the number of bonds which connect the occupied sites,
Eq. (<ref>) can be cast into

    Ξ_lg(T,μ, 12)
    =
    ∑_k_1=0^12∑_k_2=0^18 g(k_1,k_2) z^k_1exp(-V/Tk_2).

The only non-zero values of the function g(k_1,k_2) are as follows:

    g(0,0)=1;
       
    g(1,0)=12;
       
    g(2,0)=48, g(2,1)=18;
       
    g(3,0)=76, g(3,1)=108, g(3,2)=36;
    
    
    g(4,0)=45, g(4,1)=168, g(4,2)=207, g(4,3)=72, g(4,4)=3;
    
    
    g(5,0)=12, g(5,1)=48, g(5,2)=276, g(5,3)=276, g(5,4)=168, g(5,5)=12;
    
    
    g(6,0)=2, g(6,1)=0, g(6,2)=42, g(6,3)=212, g(6,4)=342, g(6,5)=264, g(6,6)=62;
    
    
    g(7,3)=12, g(7,4)=18, g(7,5)=276, g(7,6)=276, g(7,7)=168, g(7,8)=12;
    
    
    g(8,6)=45, g(8,7)=168, g(8,8)=207, g(8,9)=72, g(8,10)=3;
    
    
    g(9,9)=76, g(9,10)=108, g(9,11)=36;
       
    g(10,12)=48, g(10,13)=189;
       
    g(11,15)=12;
       
    g(12,18)=1.


This representation allows to clarify the degeneracy of the first excited state reported in Table <ref>:
According to the elaborated picture it is given by the value of g(k_1,1), k_1=2,3,4,5.
Furthermore, 
g(6,1)=0, 
i.e., one cannot place 6 occupied sites on the 12-site lattice in Fig. <ref> to have only 1 bond connecting the occupied sites.
The smallest number of bonds connecting occupied sites is 2 and g(6,2)=42:
This explains the value of the energy gap Δ=2 and the 42-fold degeneracy of the first excited state.

