Thermodynamic properties of the macroscopically degenerate tetramer-dimer phase of the spin-1/2 Heisenberg model on the diamond-decorated square lattice
Katarina Karlova Andreas Honecker Nils Caci Stefan Wessel Jozef Strecka Taras Verkholyak








katarina.karlova@cyu.fr
Laboratoire de Physique Théorique et
Modélisation, CNRS UMR 8089, CY Cergy Paris Université, Cergy-Pontoise, France


Laboratoire de Physique Théorique et
Modélisation, CNRS UMR 8089, CY Cergy Paris Université, Cergy-Pontoise, France


Laboratoire Kastler Brossel, Collège de France, CNRS, École Normale Supérieure - Université PSL, Sorbonne Université, 75005 Paris, France


Institute for Theoretical Solid State Physics, RWTH Aachen University,
Otto-Blumenthal-Str. 26, 52074  Aachen, Germany


Institute of Physics, Faculty of Science, P. J. Šafárik University, Park Angelinum 9, 04001 Košice, Slovakia


Institute for Condensed Matter Physics, NASU, Svientsitskii Street 1, 790 11, L'viv, Ukraine




§ ABSTRACT 
The spin-1/2 Heisenberg antiferromagnet on the diamond-decorated square lattice in the presence of a magnetic field displays  various quantum phases including the Lieb-Mattis ferrimagnetic, dimer-tetramer, monomer-dimer, and spin-canted phases, in addition to the trivial fully saturated state. Thermodynamic properties of this model are investigated using several complementary analytical and numerical methods such as exact diagonalization up to the systems of 40 spins, an effective monomer-dimer description, sign-problem-free quantum Monte Carlo simulations for up to 180 spins, and a decoupling approximation. Our particular attention is focused on the parameter region favoring the dimer-tetramer phase. This ground state can be represented by a classical hard-dimer model on the square lattice and retains a macroscopic degeneracy even under a magnetic field. However, the description of the low-temperature thermodynamics close to the boundary between the macroscopically degenerate dimer-tetramer and the non-degenerate monomer-dimer phases requires an extended classical monomer-dimer lattice-gas model. Anomalous thermodynamic properties emerging in the vicinity of the dimer-tetramer phase are studied in detail. Under the adiabatic demagnetization we detect an enhanced magnetocaloric effect promoting an efficient cooling to absolute zero temperature, provided that the system reaches the dimer-tetramer ground state at zero field. 




Thermodynamic properties  of the macroscopically degenerate tetramer-dimer phase of the spin-1/2 Heisenberg model on the  diamond-decorated square lattice
    Taras Verkholyak
    September 2, 2024
==========================================================================================================================================================





§ INTRODUCTION

Over the past decades significant attention has been dedicated to the investigation of low-dimensional quantum magnets situated on geometrically frustrated lattices <cit.>. The hallmark of classical highly frustrated magnets are
macroscopically degenerate  ground states <cit.>. 
Quantum fluctuations instead 
possess the remarkable capability to manipulate and potentially eliminate such degeneracies through a phenomenon known as the quantum order-by-disorder mechanism <cit.>.

Quantum spin systems thus often exhibit  ordered ground states such as long-range antiferromagnetism <cit.> or valence bond solids <cit.>. However, the interplay of  frustration and quantum fluctuations can occasionally also give rise to the emergence of exotic quantum disordered states, most prominently in quantum spin liquids <cit.>. Characterized by  long-range entanglement, quantum spin liquids indeed exhibit a plethora of intriguing phenomena, including the fractionalization of excitations and the emergence of topological order <cit.>. Nevertheless, it is also feasible for frustrated quantum magnets to stabilize quantum disordered phases that retain a peculiar macroscopic degeneracy in the absence of quantum spin-liquid characteristics. 


Frustrated quantum spin systems constitute a challenge for theory, for example because the infamous sign problem <cit.>
precludes straightforward and efficient quantum Monte Carlo simulations.
For this reason, tailor-made models such as the Rokhsar-Kivelson <cit.> and Kitaev <cit.> models have been proposed in order to allow for obtaining an exact solution.
However, geometric frustration can sometimes also be beneficial rather than detrimental: if kinetic energy is significantly or even completely suppressed by destructive quantum interference, the quantum problem may be mapped to an effective classical dimer model <cit.> or at least well approximated by a quantum dimer model <cit.>, thus rendering a controlled approximation if not an exact treatment of the low-energy physics possible.
Examples of such phenomena include but are not limited to
models with exact dimer ground states in one <cit.> and two dimensions <cit.>,
the one-third plateau in the kagome quantum antiferromagnet <cit.>,
and localized magnons appearing just below the saturation field of highly frustrated antiferromagnets
<cit.>.



In this respect, the Heisenberg model on the diamond-decorated square lattice, sketched in Fig. <ref>, serves as an illustrative example of the interplay between frustration, quantum fluctuations, and emergent phenomena.
Previous investigations into the ground-state characteristics of this model
in the absence of a magnetic field
demonstrate the presence of a Lieb-Mattis ferrimagnetic phase, a monomer-dimer phase, and a dimer-tetramer phase <cit.>. The latter two phases exhibit a macroscopic ground-state degeneracy in the absence of a magnetic field and thus deserve particular attention.
The dimer-tetramer ground state retains its macroscopic degeneracy even in the presence of a finite magnetic field (in contrast to the monomer-dimer state), highlighting its exceptional character <cit.>. This robust degeneracy can furthermore be related to the classical hard-dimer model on the square lattice, emphasizing its connection to well-known models in condensed-matter physics <cit.>.
Hence, the dimer-tetramer phase, characterized by its short-range ordering and robust macroscopic degeneracy over an extended parameter regime, offers an intriguing perspective on frustrated quantum magnetism. 
In this paper, we therefore examine the thermodynamic and magnetocaloric properties of the dimer-tetramer phase of the spin-1/2 Heisenberg model on the diamond-decorated square lattice, including a finite magnetic field. 
Indeed, we will show that this phase presents an ideal candidate for applications such as magnetic cooling during the process of adiabatic demagnetization <cit.>.

The further organization of this paper is as follows: the Hamiltonian and ground-state phase diagram of the spin-1/2 Heisenberg model on the diamond-decorated square lattice are reviewed in Sec. <ref>. Our analytical and numerical approaches will be introduced in Sec. <ref>. Results are then presented in Sec. <ref> and a final summary is provided in Sec. <ref>.






																		



§ HEISENBERG MODEL ON THE DIAMOND-DECORATED SQUARE LATTICE



In the following, we consider the spin-1/2 Heisenberg model on the diamond-decorated square lattice depicted in Fig. <ref>, and defined through the  Hamiltionian 

    H   =    J_1 ∑_i=1^N [ S_i,1·(S_i,2+S_i,3+S_i,4+S_i,5
              +
    	S_i-x̂,2+S_i-x̂,3+S_i-ŷ,4+S_i-ŷ,5)
    	]
          + J_2 ∑_i=1^N (S_i,2·S_i,3 +
    	S_i,4·S_i,5)
          - h∑_i=1^N∑_μ=1^5 S_i,μ^z .


The spin-1/2 operator S_i,μ= (S_i,μ^x, S_i,μ^y, S_i,μ^z) is assigned to the μth spin of the ith unit cell. The indices i-x̂  and i-ŷ  refer to the unit cells immediately to the left and below the ith unit cell, respectively. We consider a finite lattice consisting of N unit cells, containing N_s = 5N spins, and we apply periodic boundary conditions. Typically, we utilize square lattices with lattice sizes L_x, L_y, and the total number of unit cells N=L_xL_y. The two distinct exchange interactions J_1  and J_2 are depicted in Fig. <ref> by thin and thick lines, respectively. The magnetic-field  (h) term  represents the standard Zeeman coupling of the spin degrees of freedom to an external magnetic field.





Our investigations center on the thermodynamic properties of the dimer-tetramer phase. However,  we start by providing a brief overview of the ground-state phase diagram and the characteristics of various phases occurring in the low-field regime. The ground-state phase diagram of the spin-1/2 Heisenberg model on the diamond-decorated square lattice is depicted in Fig. <ref> and was extensively described previously <cit.>. It  exhibits five phases in the parameter plane  spanned by  the interaction ratio J_2/J_1 and the reduced magnetic field h/J_1. In addition to a saturated paramagnetic (PM) state and the spin-canted (SC) phase, the phase diagram features three distinct phases: the Lieb-Mattis (LM) ferrimagnetic, the dimer-tetramer (DT), and the monomer-dimer (MD) phases. In the LM phase the Heisenberg dimers form triplets and the overall state is reminiscent of classical ferrimagnetic order.
The MD and DT phases exhibit fragmentation. In both cases, the overall wave function can be expressed as a product of individual fragments <cit.>. In the MD phase, singlet dimers (depicted in Fig. <ref> as yellow ovals) coexist with free monomeric spins in zero field (h=0) and fully polarized monomeric spins in nonzero field.
More specifically,  the ground state can be presented as follows: 

    |MD⟩={[ ∏_i=1^N|σ⟩_i,1⊗∏_d=1^2N|s⟩_d,  σ∈{↑,↓},  h=0;        ∏_i=1^N|↑⟩_i,1⊗∏_d=1^2N|s⟩_d,   h> 0, ].

where |s⟩_d=(|↑⟩|↓⟩-|↓⟩|↑⟩)/√(2) represents the dimer-singlet state.
The DT phase is particularly intriguing because it comprises a macroscopic number of states that contain both singlet dimers and singlet tetramers.
The singlet tetramers contain in turn triplet dimers, as depicted in Fig. <ref> by small and large yellow ovals, respectively. Furthermore, each singlet tetramer is surrounded by singlet dimers within the DT ground-state manifold. 
As a result, overall these DT states are singlets and are of the following form:

    |DT⟩   =   '∏_{d}^∏_d'≠ d |t⟩_d |s⟩_d',
    	
    
    	|t⟩_d    =   1/√(3)(|↑⟩_i,1|↓↑⟩_d|↓⟩_i',1+|↓⟩_i,1|↑↓⟩_d|↑⟩_i',1) 
          -1/2(|↑⟩_i,1|↑↓⟩_d|↓⟩_i',1 + |↑⟩_i,1|↓↓⟩_d|↑⟩_i',1. 
            +.|↓⟩_i,1|↑↑⟩_d|↓⟩_i',1+ |↓⟩_i,1|↓↑⟩_d|↑⟩_i',1)  ,

where '∏_{d}^ denotes the product over a subset of dimers that are not nearest neighbors. 
This phase  exhibits a macroscopically high degeneracy that is related to the hard-dimer model on the square lattice <cit.>.
We will be especially interested in the regime  near the quantum phase transition line between the MD and DT phases. This transition line is given by the condition

    h = 2J_1 - J_2,   if 1.487 J_1≲ J_2 ≤ 2J_1,

and here,  the model is expected to exhibit an enhanced magnetocaloric effect, as we explain further below.  

We note that the finite-temperature properties in the vicinity of the MD-LM transition are well captured by a
Ising-Heisenberg variant of the model <cit.>. By contrast, the DT phase (which is the focus of the present investigation) owes its existence to quantum fluctuations on the monomeric spins and is thus absent in the Ising-Heisenberg variant.




§ METHODS



The overall 
thermodynamic properties and thermal phase transitions for various phases of this model were examined previously <cit.>. Here, we focus on the macroscopically degenerate DT phase, which is investigated using exact diagonalization, sign-problem-free quantum Monte Carlo simulations, an effective monomer-dimer (EMD) description, as well as a spin-star decoupling approximation. 
We also show some density-matrix-renormalization group (DMRG) ground-state
results, but no new data has been produced for the present purposes such
that we refer to Ref. <cit.> for technical details on these DMRG
calculations.



 §.§ Effective monomer-dimer (EMD) model


In order to formulate the effective model at the boundary between the MD and DT phases, we 
 recall some notation <cit.>. 
The energy of singlet dimers of the MD state (<ref>) can be calculated as follows:

    E^(0)_MD = 2Nε_sd = -3N/2J_2,
 
where ε_sd=-3/4J_2 is the energy of one singlet dimer.
To obtain the energy of the MD state at nonzero fields, we need to add the contribution of the Zeeman term for the monomer spins to Eq. (<ref>):

    E_MD = E^(0)_MD - Nh/2.
 
The DT phase of the states (<ref>) consists of N/2 singlet tetramers and 3N/2 singlet dimers.
Thus, its energy is

    E_ DT = N(3/2ε_sd + 1/2ε_st)
    	= - N (J_2 + J_1),

where ε_st=-2J_1+1/4J_2 is the energy of the singlet tetramer.
The condition E_MD = E_ DT yields the phase-transition line given by Eq. (<ref>). 
The ground state is macroscopically degenerate on the phase-transition line. It contains all possible configurations of two states of decorated diamonds: (i)  singlet dimers, (ii)  singlet tetramers containing a dimer triplet. 
For all states the condition applies that the singlet tetramers cannot have common edges. 

The singlet-dimer-to-singlet-tetramer excitation is given by
the gap

    Δ_H,V = ε_st - ε_sd = J_2^H,V - 2J_1.

Here we have introduced a distinction between horizontal (H) and vertical dimers (V)
with two different interaction constants J_2^H and J_2^V
in anticipation of a transfer-matrix approach that will treat the two spatial directions differently.

It should be noted that the two monomer spins in a diamond do not contribute to the Zeeman term in case they are included in the singlet tetramer. 
As a result, the effective low-temperature model can be considered as a monomer-dimer problem on the square lattice <cit.>, where dimers can be located on the bonds of a square lattice connecting nearest-neighbor vertices (which corresponds to the singlet tetramer state of the  diamond-decorated square lattice) and each vertex may host no more than one dimer (i.e., different singlet tetramers cannot have common spins). This correspondence is illustrated in Figs. <ref>(a), (b). 
Defining the energy with respect to the dimer singlets, the partition function of such a model can be written as follows:

    Z = e^-β E^(0)_ MD∑_C x^H y^V z^M,

where ∑_C denotes the sum over all possible configurations of dimers on the square lattice; 
M, H, V are the numbers of monomers and horizontal and vertical dimers (which obey the restriction 2H+2V+M = N); 
x=e^-βΔ_H, y=e^-βΔ_V, z=2cosh(β h/2) are the respective activities,
β = 1/T is the inverse temperature (where we set k_B=1). 
The monomer activity z corresponds to the partition function of a single monomer spin in an external field.
Recall that we have introduced the anisotropic excitations Δ_H and Δ_V for the horizontal and vertical dimers in order to identify the corresponding contributions in Eq. (<ref>) explicitly, but fix Δ_H=Δ_V (J_2^H=J_2^V=J_2) in all final calculations.

After introducing α=x/y, γ=z/y^1/2, we represent the partition function as 

    Z = e^-β E^(0)_ MDy^N/2∑_Cα^H γ^M,

where α=e^-β(Δ_H-Δ_V),  γ=2cosh(β h/2)e^βΔ_V/2.

One can prove without solving the model that the free energy of the monomer-dimer system is analytic such that
this model cannot exhibit any finite-temperature phase transition <cit.>.



In order to calculate the partition function, we resort to the transfer-matrix formalism <cit.> (see also Refs. <cit.>). The transfer matrix proceeds on two neighboring rows of vertical bonds. Vertical dimers are identified with spins up (↑), and the absence of dimers on vertical bonds corresponds to spin down (↓) (see Fig. <ref>(c)). The transfer matrix in operator form contains three contributions, V=V_3 V_2 V_1, with

    V_1     =    ∏_i=1^L_xσ^x_i,
    
    
    V_2     =    exp(γ∑_i=1^L_xσ^-_i),
    
    
    V_3     =    exp(α∑_i=1^L_xσ^-_iσ^-_i+1),

where σ^x_i and σ^±_i = σ^x_i ± iσ^y_i are the corresponding Pauli matrices. 
Here, V_1 counts all possible configurations of closely packed vertical dimers, where σ^x_i fixes the condition that the configurations with two spins up in the vertical direction, i.e., two neighboring dimers, are forbidden. 
The second term V_2 describes the creation of a monomer instead of a vertical dimer upon the action of the σ^-_i operator. 
In a similar manner, V_3 corresponds to the creation of horizontal dimers by the annihilation of pairs of neighboring vertical dimers. 
It is easy to check that V= V^† is hermitian, since V_1 V^†_i = V_i V_1 for i=2,3. 
For the algebraic calculations it is more convenient to consider V^2,

    V^2 = V_3 V_2 V_1 V_3 V_2 V_1 = V_3 V_2 V^†_3 V^†_2,

where V^†_2 =exp(γ∑_i=1^L_xσ^+_i), 
V^†_3 =exp(α∑_i=1^L_xσ^+_iσ^+_i+1).
The partition function can then be written as 

    Z = e^-β E^(0)_ MD y^N/2 Tr V^L_y.







An exact solution of this problem is available only for the pure dimer model (γ=0) <cit.>,
where one can establish an equivalence to a generalized XY model <cit.>. 
Then, the transfer matrix is  of the form V=V_3 V^†_3, and  can be diagonalized within the Jordan-Wigner fermionization scheme <cit.>. In the isotropic case, α=1, one obtains the entropy per unit cell of the closely-packed dimer model,  
s_d≈ 0.29156 <cit.>. 

Here, we consider the EMD model where all the terms (<ref>) in the transfer matrix are taken into account. Therefore, we proceed with a numerical calculation of the maximal eigenvalue Λ_max of the transfer matrix.
Using the Lanczos procedure <cit.>, we are able to find Λ_max for systems with a number of unit cells in the horizontal direction L_x up to 16. 
In the limit L_y→∞, only the maximal eigenvalue Λ_max remains essential compared to the  others, thus giving us the 
free energy per unit cell as follows:

    f    =    -1/Nβln Z = f_0 -1/β L_xlnΛ_max,
    	
    
    	f_0 
    	   =    -3/4 J_2^H - 1/4 J_2^V  - J_1.

All other thermodynamic quantities are finally obtained by numerical differentiation. 






 §.§ Spin-star decoupling approximation





The spin-star decoupling approximation is developed from the idea of decomposing the total Hamiltonian of the spin-1/2 Heisenberg model on the diamond-decorated square lattice (<ref>) into a sum of non-commuting cluster Hamiltonians H^0=∑_j=1^N  H_j^0 ascribed to individual spin-star clusters. As a starting point of our calculation we will therefore consider the Hamiltonian H_j^0 of a single spin-star cluster, which is composed of one monomeric spin and its four neighboring dimers. The considered spin-star cluster is schematically illustrated on the left-hand-side of Fig. <ref> including simplified notation and is mathematically defined as follows

    H_j^0   =    J_1S_0·(S_1+S_2+S_3+S_4+S_5+S_6+S_7+S_8)
          + J_2/2(S_1·S_2+S_3·S_4+S_5·S_6+S_7·S_8).

The inclusion of the factor of 1/2 in the second term ensures that the intra-dimer interaction J_2 will not be double-counted if one performs a summation over the spin-star Hamiltonians H_j^0 at the end of our calculation. Note furthermore that we have omitted the Zeeman term because this term can be trivially added at the end of the calculation as it only shifts the respective eigenenergies thanks to the total Hamiltonian and the Zeeman term commuting with each other. Within the spin-star cluster decoupling, the total energy of the spin-1/2 Heisenberg diamond-decorated square lattice will be obtained as the sum of eigenenergies of the spin stars given by the cluster Hamiltonian (<ref>) upon simply disregarding the non-commutative nature of the cluster Hamiltonians [ H_i^0, H_j^0]≠ 0 for neighbors i and j.

 Owing to the local conservation of the total spin on each dimer, it is convenient to express the cluster Hamiltonian (<ref>) in terms of composite spin operators of four dimers denoted as S_D_i = S_2i-1 + S_2i (i=1, …,4), as well as the total spin of the dimers S_D = ∑_i=1^4S_D_i
 and the total spin of the spin star denoted as S_t=S_0+S_D. Consequently, the eigenenergies of the spin-star cluster Hamiltonian (<ref>) can be obtained from the eigenenergies of a simpler composite spin star illustrated on the right-hand-side of Fig. <ref>, which depend only on the total quantum spin number S_t, the composite quantum spin numbers of the individual dimers S_D_i, and the total quantum spin number of all four dimers S_D, 

    E_S_t, S_D, S_D_i^0   =   J_1/2[S_t(S_t+1)-3/4- S_D
    (S_D+1)]
          - 3/2J_2+J_2/4∑_i=1^4 S_D_i (S_D_i + 1).

We note
that each composite spin dimer  can be either in a singlet state characterized  by S_D_i=0 or a triplet state with S_D_i=1, which means that the spin star can have from zero up to four triplets.  Assuming  antiferromagnetic couplings J_1>0 between the central monomeric spin S_0 and the four neighboring spin dimers S_D_i, the energetically most favorable states of the spin star in each sector correspond to the highest
multiplicity of the dimer spins S_D^ max = ∑_i S_D_i, given by the eigenenergies E_S_t, S_D^ max^0 with

    E_1/2,0^0   =    -3/2J_2,
    
    E_1/2,1^0   =    -J_1-J_2,
    
    E_3/2,2^0   =    -3/2J_1-J_2/2,
    
    E_5/2,3^0   =    -2J_1,
    
    E_7/2,4^0   =    -5/2J_1+J_2/2.

The total number of triplets within this set of eigenstates of the spin star thus coincides with the sum of the composite quantum spin numbers S_D. 

The total energy of the spin-1/2 Heisenberg model on the diamond-decorated square lattice can be then determined by summing the energies of individual spin stars, whereby the states of the dimers can be effectively described within a quasi-particle formalism when assigning the occupation number n_i=1 to a dimer-triplet state and n_i=0 to a dimer-singlet state, respectively. By considering only the five lowest-energy eigenstates of the spin star, ranging from zero up to four triplet-dimer states as given by Eq. (<ref>), the eigenenergies can be expressed in terms of the occupation numbers of the dimer-triplet states,

    E_n_1,…, n_4^0 =-3/8J_1-3/2J_2+J_2/2∑_i=1^4n_i+J_1/2∏_i=1^4(1-n_i)
    
          +J_1/2[(∑_i=1^4  n_i-1/2)(∑_i=1^4  n_i+1/2)
    -∑_i=1^4  n_i (∑_i=1^4  n_i+1)],

which can be further simplified into the following compact form

    E_n_1, …, n_4^0=   -   1/2(J_1+3J_2)+1/2(J_2-J_1)∑_i=1^4  n_i 
       +   J_1/2∏_i=1^4(1-n_i).

The last term represents the correction to the energy of free monomeric spins surrounded by four adjacent dimer singlets. 
The eigenenergy E_n_1, …, n_4^0 corresponds to the spin multiplet with the total quantum spin number

    S_t= ∑_i=1^4 n_i+ ∏_i=1^4  (1-n_i)-1/2.

Summing the eigenenergies of the spin stars and adding the respective Zeeman term yields the following formula for the overall energy of the spin-1/2 Heisenberg model on the diamond-decorated square lattice
E_T({n_i}) = E_T^0({n_i})-h S_T^z, where E_T^0({n_i}) represents the respective zero-field energy depending on the set of all occupation numbers { n_i}:


    E_T^0({n_i})   =   -N/2(J_1+3J_2)+(J_2-J_1)∑_i=1^2Nn_i
          +(J_1/2)∑_j=1^N ∏_i∈ (1-n_j,i).

The last term in Eq. (<ref>) provides a correction to the energy when all four dimers surrounding the jth monomeric spin are in a singlet state. This situation is depicted in Fig. <ref>, where the central monomeric spin is surrounded by four occupation numbers assigned to  composite spins of adjacent dimers. If all four composite spins of the dimers are in the singlet state, the correction term accounts for this specific configuration, ensuring that the energy is accurately represented under this specific condition. In Eq. (<ref>), the notation i∈ indicates that the index i runs over the four sites of a square plaquette surrounding the central monomeric spin, as depicted in Fig. <ref>. Finally, S_T^z is the z-component of the total spin S_T of the whole system, which can be expressed in terms of the occupation numbers of the dimer-triplet states

    S_T= -N/2 +∑_i=1^2Nn_i+∑_j=1^N ∏_i∈ (1-n_j,i).

The first term accounts for the contribution of the monomeric spins, while the second term accounts for the contribution of the dimeric spins.  Hence, the second term involves a summation over the occupation numbers assigned to the dimer-triplet states and thus corresponds to the total number of triplets N_ trip.. The summation in the third term is a correction stemming from all free
monomeric spins. In fact, the third term involves a projection operator for the monomeric spin surrounded by four dimer-singlet
states and thus corresponds to the total number of free monomeric spins N_ free. 


The free energy of the spin-1/2 Heisenberg model on the diamond-decorated square lattice within the spin-star decoupling approximation can then be calculated according to the formula:

    f_ star = -1/Nβln∑_{n_i}∑_S_T^z=-S_T^S_Texp[-β E_T^0({n_i})+β h S_T^z].

We emphasize that the formula (<ref>) for the total quantum spin number S_T tacitly assumes that all free monomeric spins are fully polarized by the magnetic field. Therefore, in the zero-field limit  the formula (<ref>) for the free energy does not take into account  the degeneracy from the free monomeric spins (i.e., monomeric spins surrounded by four singlet dimers), which have a significant role in determining the zero-point entropy and low-temperature thermodynamics in  the absence of a magnetic field.

To improve the zero-field estimate of the free energy within the spin-star decoupling approximation, we therefore modified the formula (<ref>) for the free energy by accounting for the respective degeneracy factor of the free monomeric spins

    f_ star^0 = -1/Nβln∑_{n_i}∑_S_T^z=-S_T^0^S_T^02^N_ freeexp[-β E_T^0({n_i})],

where N_ free=∑_j=1^N ∏_i∈ (1-n_j,i) denotes the total number of the free monomeric spins and the overall spin multiplicity S_T^0 in zero field can be expressed as   

    S_T^0=N_ free-N/2 + ∑_i=1^2Nn_i.


The spin-star decoupling approximation recovers the transition fields around the DT  phase in the ground-state phase diagram (see red broken lines in Fig. <ref>(b)). By setting all n_i=0  in Eq. (<ref>), we actually get the exact value for the ground-state energy of the MD  phase (<ref>). Requiring that there is a single triplet dimer in each spin-star configuration, we instead obtain the exact energy of the DT phase (<ref>).
Calculating the ground-state energy of the LM phase (all n_i=1) within the spin-star approach, E_ LM=-N/2(5J_1 - J_2 + 3h), yields the transition field between the LM and DT phases as h=J_2 - J_1, which  is very close to the DMRG results shown in Fig. <ref>(b) (cf. the red broken and blue solid lines in Fig. <ref>(b)).  
Recall that on the phase boundary the energy of any configuration that contains at least one triplet in the spin-star cluster has the same energy. This  leads to an additional macroscopic degeneracy (higher than the degeneracy of the DT phase). This artifact of the decoupling approximation influences the low-temperature behavior entropy discussed in Sec. <ref>.


The decoupling approximation based on the spin-star clusters was developed under two fundamental assumptions, (i) ignoring the non-commutative nature of the cluster Hamiltonians (<ref>) and (ii) subsequently considering only the five energetically most favorable states of the spin star given by Eq. (<ref>). Despite these simplifications, the spin-star cluster decoupling  offers a reasonable approximation at sufficiently low temperatures whenever the system is driven towards the fragmented MD  or DT phases, whereas the collective nature of the LM ferrimagnetic phase is  
not rigorously accounted for within this treatment. Although we have attempted to improve precision of the spin-star approximation in zero magnetic field by considering degeneracies pertinent to free monomeric spins, this calculation procedure should be regarded as much more reliable and precise in
the presence of a magnetic field due to Zeeman splitting of the relevant spin multiplets. It indeed turns out that the quantitative agreement between the results derived from the spin-star approximation and a numerical treatment of the full model progressively extends to higher temperatures as the magnetic field increases,
see Sec. <ref>.




 §.§ Exact diagonalization






The full N_s=30 spectra were obtained previously <cit.> and can be re-evaluated for the present purposes.
However, here we push exact diagonalization (ED) further to the N_s=40 system sketched in
Fig. <ref>. We have followed the same strategy as in Ref. <cit.>, exploited conservation of the total spin
on each dimer and thus passed to an effective Lieb lattice, represented by the yellow spheres in the background of Fig. <ref>.
For the present N_s=40 system, computer enumeration then yields 191 topologically nonequivalent arrangements of spin triplets and singlets
on the dimers with degeneracies ranging from 1 (for all dimers in the triplet state) to 1152. 
We were then able to compute the full spectra in the sectors with up
to N_ trip.=8 triplets. This includes in particular the sector with N_ trip.=4 triplets that contains the dimer-tetramer phase for the N_s=40 geometry.
For the sectors with 9 triplets or more, we further diagonalized high-S_z sectors completely and obtained low-energy levels
using the Lanczos procedure <cit.> combined with the strategy of Ref. <cit.>. Specifically, in the Lieb-Mattis regime that corresponds to 16 triplets in the N_s=40 system,
we were only able to obtain full spectra for S^z ≥ 13 and had to contend ourselves with low-lying levels for S^z ≤ 12,
i.e., the ferrimagnetic state with magnetization 3/5 and below.
Furthermore, we reconstructed the spin multiplets using S^z classification and spin inversion symmetry in the S^z=0 sector;
the latter permits us to avoid diagonalizations in the S^z=1 sector, which would usually have the highest dimension.
All in all, we were able to recover an entropy of 36.082 ln(2), i.e., 90.2% of the total entropy of the N_s=40 system.

We note that going to the N_s=40 system is particularly important in the vicinity of the Lieb-Mattis phase.
Specifically, the lowest state in the sector with N_ trip.>N_s/10 triplets turns out to be in the sector with
spin S=N_ trip.-N_s/10 and approaches the Lieb-Mattis limit S=3N_s/10 for N_ trip.=2N_s/5.
At J_2=0, we find a ground-state energy of -2.46115 J_1 per unit cell in the S=12, N_ trip.=16 sector for N_s =40.
The finite-size error is about an order of magnitude smaller than the N_s=30 estimate of -2.46809 J_1 as compared to the estimates
-2.46083 J_1 <cit.> and -2.46000 J_1 per unit cell <cit.>, respectively. This indicates that the
N_s=40 system can be considered significantly closer to the thermodynamic limit than the N_s=30 one, even if we increased
the number of spins only by a third. This large change can be attributed to the N_s=30 system still having additional symmetries that arise from the periodic
boundary conditions imposed on this small system whereas such particular symmetries are absent in the N_s=40 geometry,
rendering the latter geometry significantly more generic.
As another indicator, we mention that the ground-state entropy per site is s≈0.0880 and 0.0795 for N_s=30 and 40, respectively,
compared to s≈0.0583 in the thermodynamic limit <cit.> (the corresponding ground-state degeneracy is 14 and 24, respectively).




Figure <ref> shows the gaps to the lowest excitation above the DT ground-state for the N_s=40 system (N_ trip.=4).
For J_2 ≳ 1.5 J_1,
the gap to the N_ trip.=3 excitation coincides with the DT-MD phase transition in the phase diagram Fig. <ref>, as in fact
in this case states with any number N_ trip.≤ N_s/10 collapse.
On the other hand, the field-induced DT-LM transition in Fig. <ref> is a true first-order transition. Thus, the closing of the N_ trip.=5 gap in Fig. <ref> is preempted by a direct transition to a N_ trip.=16 configuration, reproduced from Fig. <ref> by the dashed line in
Fig. <ref>. As a consequence, the maximum of the gap is found for J_2/J_1 ≈ 1.4
whereas the maximal extension of the DT phase in a magnetic field occurs for J_2/J_1 ≈ 1.5.
Furthermore, one observes an N_ trip.=6  bound state coming down in energy when approaching the DT-LM transition at h=0. The nature of the DT-LM transition at h=0 is less clear to decide, at least based on the N_s=40 ED data, i.e., one does observe higher-N_ trip. / higher-S states coming down in energy, but it is not clear if they actually collapse in the thermodynamic limit, or if the DT-LM transition remains a conventional first-order transition at h=0.






 §.§ Quantum Monte Carlo

In addition to ED, we can also use QMC in order to access thermodynamic properties of the spin-1/2 Heisenberg antiferromagnet on the diamond-decorated square lattice. In particular, this approach allows us to obtain unbiased results for system sizes that  extend beyond those accessible to ED. The QMC method that we used for our investigation has already been employed in our previous finite-field study of the spin-1/2 Heisenberg antiferromagnet on the diamond-decorated square lattice <cit.> and therefore we summarize it only briefly.
 
More specifically, our QMC approach is based on the stochastic series expansion (SSE)  
method with directed loop updates <cit.>. In order to avoid 
the  sign problem, i.e., an exponential drop of the statistical accuracy of the QMC simulations at low temperatures and large system sizes <cit.> due to the presence of geometric frustration, we avoid working in the conventional local spin-S^z basis for the model at hand. 
 Instead, we eliminate the sign problem upon using  appropriate basis states after decomposing the Hamiltonian into separate terms based on dimers (or trimers) <cit.>. 
For certain so-called fully frustrated models
the sign problem can indeed be completely eliminated, cf. Refs. <cit.> for the spin-dimer and Ref. <cit.> for the spin-trimer basis, respectively. 
 
In case of the diamond-decorated square lattice a finite value of the  coupling J_2 leads to geometric frustration. 
In order to avoid the associated sign problem, we follow Ref. <cit.> and  treat all J_2-dimer spins in the spin-dimer basis, but keep using the local S^z-basis for the monomer spins. Using this combined five-site cluster-basis, we can simulate our model without a sign-problem within the SSE framework, based on the abstract operator loop update introduced in Ref. <cit.>. We refer to Ref. <cit.> for further details on this QMC approach. 
Here, we performed QMC simulations for systems with 
up to 6×6 unit cells.








§ RESULTS






 §.§ Dimer triplet density, entropy, and specific heat above the DT phase





Let us start our discussion by examining the dimer triplet density, entropy, and specific heat of the spin-1/2 Heisenberg model on the diamond-decorated square lattice in zero magnetic field. We will compare results from the different methods
introduced in Sec. <ref>. For the spin-star decoupling approximation we will use the zero-field formula Eq. (<ref>) for the free energy.

Figure <ref> presents results  for the interaction ratio J_2/J_1=1.7. 
At zero temperature, the dimer triplet density is 0.25, reflecting the nature of the DT state, where 1/4 of the dimers are in the dimer-triplet state within the singlet tetramers and 3/4 of the dimers remain in the dimer-singlet state. Due to the relatively lower energy cost of creating a dimer-singlet state compared to an additional dimer-triplet state for the interaction ratio J_2/J_1=1.7, a decrease in the dimer-triplet density can be observed at 
low temperatures terminated at a round local minimum. With further increase of temperature, the dimer-triplet density tends towards the asymptotic value 0.75 reached in the limit of infinite temperature due to the three-fold degeneracy of the triplet state in zero magnetic field. 

A comparison of the dimer-triplet density obtained for N_s=30 and N_s=40 spins from exact diagonalization (ED) reveals 
small
finite-size effects. ED for N_s=40 spins provides an interpolation between the spin-star decoupling approximation and the EMD model at low temperatures and QMC at high temperatures.
The excitation from a dimer-singlet state to a dimer-triplet state is neglected within the EMD model and the dimer-triplet density should accordingly always decrease with increasing temperature. On the other hand, the local minimum observed in the dimer-triplet density arises from the competition of the lower excitation energy for the dimer-singlet state and the higher degeneracy of the triplet state.
In fact, the lowest states not taken into account by the EMD model are at energies slightly above
0.6J_1 for J_2=1.7J_1 and these are responsible for the triplet density increasing again with increasing temperature after going through a minimum. However, the spin-star decoupling approximation takes only some of these into account (as reflected also by the lower entropy s compared to ED and QMC, see Fig. <ref>(b)) and approximates their energies, thus explaining the deviation from N_s=40 ED and QMC, which are in good agreement for T/J_1 ≳ 0.2.


Figure <ref>(b) shows
the entropy per spin s
as a function of temperature at zero magnetic field for the interaction ratio J_2/J_1=1.7. The zero-temperature entropy of the spin-1/2 Heisenberg model on the diamond-decorated square lattice in the parameter regime of the DT phase is given by s=1/N_SlnΩ_N_s, where Ω_N_s represents the degeneracy of the DT phase for the system size with N_s spins; for instance,
Ω_30=14, Ω_40=24, as already stated in the ED context in Sec. <ref>.
Figure <ref>(b) shows that the results obtained from the spin-star decoupling approximation for N_s=40 are in perfect agreement with ED data for the same system size at low enough temperature T/J_1 ≲ 0.1. Similar accuracy can be expected at low temperatures for the EMD model. Notably, the residual entropy of the DT phase decreases with increasing system size, but it still remains nonzero even in the thermodynamic limit N_s →∞. From the EMD model with system size 16×∞, one can for instance infer a residual entropy  s≈0.0587 that quite closely converges to the known value for the hard-dimer model on the square lattice <cit.>. 

Figure <ref>(c) displays
the specific heat per spin c as a function of temperature for the same value of the interaction ratio. A notable feature 
is the presence of a low-temperature peak respectively shoulder superimposed on the high-temperature maximum of the specific heat. As anticipated, the spin-star and EMD approaches fail to capture the behavior of the specific heat at higher temperatures as they neglect high-energy excitations above the DT phase.
On the other hand, these approaches accurately account for low-lying excitations which leads to a close agreement between the ED data for N_s=40 and the analytical results obtained from the spin-star decoupling approximation at low enough temperatures. A comparison of the peak height in the specific heat obtained from ED and QMC reveals relatively large finite-size effects at low temperatures. The EMD model provides further insight into the behavior of the specific heat at low temperatures, which effectively bridges the gap between the available QMC data (T/J_1>0.1) and extends it further towards lower temperatures. Figure <ref>(c) shows that thermal activation of the specific heat indeed
occurs at progressively lower temperatures as the system size increases such that the EMD model provides the best approximation to the thermodynamic limit at low temperatures.



We now turn to a discussion of the dimer-triplet density, entropy, and specific heat in zero magnetic field for the case  J_2/J_1=1.5  plotted in Fig. <ref>. 
Here the numerical data for the dimer-triplet density
exhibits only a shallow minimum, which is consistent with the gap for one triplet less being only slightly smaller than that for an additional triplet at  J_2/J_1=1.5, compare Fig. <ref>.
By contrast, the spin-star decoupling approximation fails to account for this minimum. The reason is that, within the spin-star decoupling approximation, the excitation energies from the DT phase to singlets and triplets are equal for this interaction ratio.
Consequently, the spin-star approach cannot explain the thermal population into singlet states for the interaction ratio J_2/J_1≤1.5. Comparing the results obtained from QMC simulations for a system size of 6 × 6 (total of 180 spins) with those from ED results for 40 spins reveals nearly perfect agreement, indicating negligible finite-size effects. 
However, the discrepancy between ED results for 30 spins and 40 spins at intermediate temperatures (T/J_1>0.1) may stem from artificial symmetries imposed for the system size of 30 spins by considering periodic boundary conditions, which are absent for the system size of 40 spins.

The entropy and specific heat
for the interaction ratio J_2/J_1=1.5 in zero magnetic field are shown in Fig. <ref>(b) and <ref>(c), respectively, and exhibit qualitatively similar behavior to that for the interaction ratio J_2/J_1=1.7. In the high-temperature regime, small finite-size effects are observed.
While the high-temperature peak of the specific heat shifts slightly towards higher values with increasing system size, the low-temperature peak is conversely suppressed towards lower values.
The results obtained from the spin-star decoupling approximation for 40 spins are slightly overestimated compared to ED data for N_s=40  at  low temperatures T/J_1 ≲ 0.1. The discrepancy at higher temperatures arises due to the presence of low-energy excitations to states that share the same energy as the Lieb-Mattis ferrimagnetic state at zero magnetic field (compare the discussion in Sec. <ref>), which is not accurately accounted for by the spin-star decoupling approximation simply ignoring the non-commutativity of the spin-star Hamiltonians. The QMC data are consistent with the intriguing double-peak temperature dependence of the zero-field specific heat with low- and high-temperature peaks located around T/J_1≈ 0.12 and 0.6, respectively.




With a yet lower value of the interaction ratio J_2/J_1, the analytical approaches provided by the spin-star and EMD models neglect more and more important low-energy excitations. For the interaction ratio J_2/J_1=1.3, both of these approaches only provide an accurate description
at very low temperatures as evidenced by the temperature dependencies of the dimer-triplet density, entropy, and specific heat depicted in Fig. <ref>. This limitation arises from the failure of these analytical approaches  to accurately capture the
degeneracy of the LM phase in zero magnetic field with possible values of the z-component of the total spin S_T^z = -S_T, -S_T+1, ..., S_T, and the low-energy excitations to these states are not properly accounted for either.
The dimer-triplet density
monotonically increases without any local minimum for the interaction ratio J_2/J_1=1.3, as evidenced by the ED data for N_s=40 spins. This behavior is due to the
lowest excitation in this parameter regime arising from an additional dimer-triplet,
compare Fig. <ref>.
The EMD approach neglects the lowest excitations and thus only gets the T=0 limit right, but fails to capture the leading low-temperature asymptotics. By contrast, the star approach does account for these excitations, even if only in an approximate fashion, and thus captures at least qualitatively the low-temperature behavior of the triplet density.

While the high-temperature peak of the specific heat shown in Fig. <ref>(c) exhibits a slight variation with increasing system size, the behavior of the low-temperature peak is less straightforward. For a system size of N_s=30, the peak is lower than for size N_s=40, but the specific heat obtained from QMC for a system size of 6×6 unit cells (N_s=180) shows the lowest peak. This indicates a more complex dependency on system size.





 §.§ Magnetization, entropy, and specific heat under magnetic field




We proceed to a discussion of the most intriguing findings
under finite magnetic fields.
Now we use the finite-field formula Eq. (<ref>) for the free energy
within the spin-star decoupling approximation.



The field dependencies of magnetization, entropy, and specific heat are illustrated for the interaction ratio J_2/J_1=1.7 in Fig. <ref> for two different temperatures T/J_1=0.01 and 0.1.
The magnetization curve at the interaction ratio J_2/J_1=1.7 exhibits magnetization plateaux at m=0 and 1/5 (full magnetization corresponding to 1/2) with the characteristics of the DT and MD phases, respectively. These can be almost perfectly described by the analytical spin-star or EMD approaches including the most important low-energy excitations above these ground states.
However, another transition from the 1/5-plateau to the 3/5-plateau corresponding to the Lieb-Mattis ferrimagnetic phase cannot be described at all by the EMD model and only qualitatively by the spin-star decoupling approximation. This finding relates to the fact that the ground-state energy of the Lieb-Mattis phase is determined within the spin-star decoupling approximation with a relative error of approximately 4.3% (the ground-state energy of the LM phase within DMRG and spin-star decoupling methods in zero field are E_ LM^ DMRG=-2.46J_1+2J_2  <cit.>
and E_ LM^ STAR =-2.50J_1+2J_2, respectively, compare also the related discussion in
Sec. <ref>).

It can be observed from Fig. <ref>(b) that the entropy starts from its nonzero residual value, which is progressively converging with increasing system size to the residual entropy of the dimer model on the square lattice per unit cell s_d≈ 0.29156 (per spin ≈ 0.058312) <cit.>. Contrary to this, one detects a sizable entropy gain reaching the specific value s=0.13256
at the transition field between the DT and MD phases, which is almost independent of the system size as corroborated by ED, the EMD model, as well as the spin-star decoupling approximation. Similar qualitative agreement between results of all three aforementioned methods is observed in the peculiar asymmetric double-peak structure of the specific heat, which can be detected in the vicinity of the magnetic-field-driven phase transition between DT and MD phases (see the inset of Fig. <ref>(c)). The lower peak height observed below the relevant transition field can be attributed to the macroscopic degeneracy of the DT ground state, which is not lifted by the magnetic field in contrast to the MD phase  that becomes non-degenerate at any nonzero magnetic field.


The results obtained from the spin-star decoupling approximation are in very good agreement with the ED data even for magnetic fields higher than  h/J_1≳ 1.0, up to the end of the Lieb-Mattis ferrimagnetic phase, and even at moderately high temperatures T/J_1=0.1 still constitute a good approximation of the field dependence of the entropy and specific heat, as depicted in Fig. <ref>(b) and (c).





Figure <ref>(a) shows  the magnetization curves
for the interaction ratio J_2/J_1=1.5.
There is a distinct magnetization jump from the DT phase to the MD phase at zero temperature, as shown in the inset of Fig. <ref>(a). However, the 1/5-plateau is sufficiently small to be smeared out at finite temperatures. Consequently, there is a sharp rise in magnetization from the 0-plateau to the 3/5-plateau even at temperatures as low as T/J_1=0.01, which becomes smoother with increasing temperature. Due to the presence of the MD phase within a narrow range of magnetic-field strengths, the EMD model is in excellent agreement with the ED data up to a magnetic field of approximately h/J_1≈ 0.54 at a very low temperature T/J_1=0.01 though it fails to account for further increases in magnetization with increasing magnetic-field strength. On the other hand, the spin-star decoupling approximation is still capable of qualitatively reproducing the ED data even at higher magnetic fields. The comparison of ED data for N_s=40 spins with QMC data for 36 unit cells (i.e., N_s=180 spins) indicates negligible finite-size effect in this parameter region.

The results for the low-temperature entropy exhibit a relevant discrepancy between the ED and spin-star decoupling approximation only in the close vicinity of the transition field h/J_1=0.5 (see Fig. <ref>(b)). The excessive macroscopic degeneracy predicted by the spin-star decoupling approximation relates to a triple phase coexistence point of the LM, DT, and MD phases, which is located exactly at the magnetic field h/J_1=0.5 for J_2/J_1=1.5.
Recall that the triple point obtained within the spin-star approach is slightly shifted towards a lower value of the magnetic field with respect to its true coordinates (see Fig. <ref>(b) and the discussion in Sec. <ref>). At a higher temperature T/J_1=0.1 the discrepancy between ED and the spin-star decoupling approximation is significantly reduced, because all relevant low-lying excited states are thermally activated. Contrary to this, we see the opposite effect for the EMD model with much lower entropy, because high-energy excitations are missing within this description. A similar effect can also be observed in the field dependence of the specific heat (
Fig. <ref>(c)). While the EMD model reproduces the ED data at low temperatures, the spin-star decoupling approximation correctly describes the double peak of the specific heat also at a higher temperature of T/J_1=0.1. At this higher
temperature, the double-peak structure of the specific heat is also confirmed by the QMC simulations
for the bigger system of 6×6 unit cells.





Last but not least, magnetization, entropy, and specific heat as a function of temperature are depicted in Fig. <ref> for the interaction ratio J_2/J_1=1.3. As discussed in Sec. <ref>, for the interaction ratio J_2/J_1≲ 1.4, the lowest excitation above the DT phase corresponds to the creation of an additional triplet (see Fig. <ref>). For this reason, the EMD model is unable to describe the finite-temperature properties
for the interaction ratio J_2/J_1=1.3.
On the other hand, our second analytical approach derived from the spin-star decoupling approximation provides a quite reliable description of all quantities above the transition field h/J_1=0.3.
The most important artifacts of the spin-star decoupling approximation are just below the transition field
and at T/J_1=0.01 
where it yields a spike in the entropy s (Fig. <ref>(b)) and a second maximum
in the specific heat c (Fig. <ref>(c)) while such pronounced maxima are absent in the ED results.
The reason for this difference is that the field-induced DT-LM transition is a true first-order transition in the full model
such that excitations are scattered over a certain energy range below the transition
while the spin-star decoupling approximation treats these as a degenerate manifold that comes down at the transition field.
Conversely, a temperature T/J_1=0.1 is sufficiently high to thermally excite these low-lying excitations
in the full model such that at this higher temperature the N_s=40 ED results also exhibit a pronounced maximum
below the transition field, both in the entropy (Fig. <ref>(b)) and the specific heat (Fig. <ref>(c)).
The latter behavior is confirmed by our QMC simulations on the bigger system of 6× 6 unit cells.




 §.§ Enhanced magnetocaloric effect due to the macroscopically degenerate DT phase






Finally, we discuss the magnetocaloric effect with a focus on the transition region between the macroscopically
degenerate DT phase and the MD state (that is non-degenerate in a finite field).
Figure <ref> shows 
density plots of the entropy in the magnetic field versus temperature plane
for the interaction ratio  J_2/J_1=1.7.
The superimposed constant-entropy curves correspond to the variation of the temperature T during an adiabatic
(de)magnetization process.

First, we discuss the ED results for a system size of N_s=40 spins shown in Fig. <ref>(a).
An enhanced magnetocaloric effect can clearly be seen in the vicinity of the magnetic field h/J_1=0.3, where the temperature sharply decreases during an adiabatic demagnetization process, whereas an inverse magnetocaloric effect associated with an increase in temperature with further reduction of the magnetic field can be mostly observed below this value of the magnetic field. 
The specific value of the magnetic field h/J_1=0.3 corresponds to the magnetic-field-driven phase transition from the MD phase to the DT phase for the interaction ratio J_2/J_1=1.7. Recall also that the residual entropy of the DT phase (per spin) for a system size of N_s=40 is equal to s=ln 24/40≈ 0.079 and, hence, the system cools down to absolute zero temperature at h/J_1=0.3, whenever the entropy at the beginning of the adiabatic demagnetization is tuned below this value. Under this condition, the temperature stays at zero below the transition field and there is no consecutive increase in temperature upon further demagnetization,
as Fig. <ref>(a) demonstrates by the absence of isentropes below the magnetic field value h/J_1≲ 0.3 (i.e., below the cyan curve). This asymmetry 
of the entropy around the transition field is caused by the special character of the DT phase, which preserves its macroscopic degeneracy 
at nonzero magnetic fields. From this point of view, the spin-1/2 Heisenberg model on the diamond-decorated square lattice offers 
promising
perspectives
for achieving ultra-low temperatures during adiabatic demagnetization.


For comparison, Fig. <ref>(b) presents the density plot of the entropy, for the same set of parameters and the same system size N_s=40, as obtained from Eq. (<ref>) within the spin-star decoupling approximation.
Remarkably, the spin-star approach is capable of predicting the aforementioned asymmetry of the entropy, and the relevant density plot reveals qualitatively the same features of the magnetocaloric effect as obtained from ED.
Moreover, the enhanced magnetocaloric effect and cooling to absolute zero temperature should persist even in the thermodynamic limit due to the macroscopic degeneracy of the DT phase.
The EMD model shows indeed
that the asymmetric behavior of the isentropes as a function of magnetic field and temperature is preserved 
for a system size comprising of 16 elementary units in one spatial direction and infinite in the other direction, see
Fig. <ref>(c). For this system size, the residual entropy of the DT phase is  s ≈ 0.0587, which is consistent with the macroscopic degeneracy of the hard-dimer model on the square lattice
<cit.>. 
 Hence, the enhanced magnetocaloric effect accompanied with cooling to  zero temperature is detected for h/J_1<0.3 whenever the entropy is tuned below its zero-point value s ≲ 0.0587  at the beginning of an adiabatic demagnetization process (i.e., all curves lying below the cyan curve). 





§ CONCLUSION


In this study, we have investigated the thermodynamic properties of the spin-1/2 Heisenberg model on the diamond-decorated square lattice, particularly focusing on the macroscopically degenerate DT phase. Using a combination of exact diagonalization, sign-problem-free quantum Monte Carlo simulations, an effective monomer-dimer model, and a spin-star decoupling approximation, we have provided a comprehensive analysis of this model's behavior under varying interaction ratios and external magnetic fields.
Our results confirm that the DT phase exhibits unique thermodynamic properties due to its macroscopic degeneracy, which is linked to the classical hard-dimer model on the square lattice. The boundary between the MD and DT phases is well described by the effective monomer-dimer model, highlighting the significance of localized excitations in determining the low-temperature behavior of this system.


A particularly interesting
result of our study is the enhanced magnetocaloric effect observed upon approaching the DT phase. This effect is particularly pronounced near the transition 
field between the MD and the DT phase, where the temperature decreases sharply during adiabatic demagnetization. This phenomenon arises due to the macroscopic degeneracy of the DT phase, which allows the system to cool to absolute zero temperature as long as the entropy content is set below the residual entropy value s ≈ 0.0583 derived from the effective hard-dimer model on the square lattice in the thermodynamic limit. This renders the DT phase a promising candidate for achieving ultra-low temperatures
by magnetic refrigeration.
The third law of thermodynamics <cit.> requires this macroscopic degeneracy to be lifted
in an experimental realization, for example by a lattice distortion that breaks the conservation of the total
spin of a dimer.
Our conclusions will nevertheless remain relevant
as long as the resultant splitting of the ground-state manifold remains at an energy scale that is below the temperature scale of interest.
The challenge is more to find a suitable compound to begin with. In this respect, it may be interesting
to consider the diamond-decorated honeycomb lattice since in this case there are known compounds that exhibit
the corresponding crystal structure, see for example Refs. <cit.>.

Our findings 
demonstrate the usefulness of effective classical statistical physics models for the low-energy and low-temperature
behavior of a highly frustrated quantum spin system and
shed light on the rich phase diagram and thermodynamic phenomena of the spin-1/2 Heisenberg model on the diamond-decorated square lattice.
Furthermore, they suggest potential applications of such highly frustrated quantum spin systems in magnetic cooling technologies.

 

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 945380. We acknowledge support under the Štefánik+ program for Slovak-France bilateral under contract SK-FR-22-0011/49880PG. 
Computing time for the QMC simulations was provided by the IT Center at RWTH Aachen University,
part of the ED computations were carried out on the “osaka” cluster at the Centre de Calcul (CDC) of CY Cergy Paris Université.
T.V. is supported through the EURIZON project (#3025 “Frustrated quantum spin models to explain the properties of magnets over wide temperature range”), which is funded by the European Union under grant agreement No. 871072. J.S. acknowledges financial support provided by Slovak Research and Development Agency under the contract No. APVV-20-0150. N.C. acknowledges support by the ANR through grant
  LODIS (ANR-21-CE30-0033).





