Mean-field phase diagram and spin glass phase of the dipolar Kagome Ising antiferromagnet
Leticia F. Cugliandolo Laura Foini Marco Tarzia





Sorbonne Université, Laboratoire de Physique Théorique et Hautes Energies, CNRS
UMR 7589, 4, Place Jussieu, Tour 13, 5ème étage, 75252 Paris Cedex 05, France.

Institut de Physique Théeorique, Université Paris Saclay, CNRS, CEA, F-91191
Gif-sur-Yvette, France

Sorbonne Université, Laboratoire de Physique Théorique de la Matière Condensée, CNRS UMR 7600, 4, Place Jussieu, Tour 13, 5ème étage, 75252 Paris Cedex 05, France


We derive the equilibrium phase diagram of the classical dipolar Ising antiferromagnet at the mean-field level on a geometry that mimics the two dimensional Kagome lattice. Our mean-field treatment is based on the combination of the cluster variational Bethe-Peierls formalism and the cavity method, developed in the context of the glass transition, and is
complementary to the Monte Carlo simulations realized in [Phys. Rev. B 98, 144439 (2018)]. 
Our results confirm the nature of the low temperature crystalline phase which is reached through a weakly first-order phase transition.
Moreover, they allow us to interpret the dynamical slowing down observed in the work of Hamp & al. as a remnant of a spin glass transition taking place at the mean-field level (and expected to be avoided in 2 dimensions).


xxx

Mean-field phase diagram and spin glass phase of 

the dipolar Kagome Ising antiferromagnet
    Marco Tarzia
    December 30, 2023
===========================================================================================





§ INTRODUCTION


Many interesting classes of classical and quantum magnetic systems are extremely constrained. Hard local constraints lead to frustration and to the impossibility of satisfying all competing interactions simultaneously <cit.>, giving rise to the existence of highly degenerate ground states <cit.>. Under certain conditions, these features produce a rich variety of collective behaviors <cit.>, unconventional phase transitions <cit.>, the emergence of a Coulomb phase with long-range correlations <cit.>, and other remarkably unusual and exotic phenomena.

On the other hand, frustration is also one of the key properties of glassy systems <cit.>, where it also arises from the fact that minimizing some local interactions leads to the impossibility of minimizing other ones <cit.>. This feature can generate rugged energy landscapes and slow dynamics even in the absence of disorder <cit.>.

It is therefore surprising at first sight that very little is known on glassy phases in geometrically frustrated magnetic systems.
One of the first tentative investigations on this subject has been performed in Ref. <cit.>,
where glassy behavior was observed in nonrandomly frustrated Ising models with competing interactions.
More recently, strong nonequilibrium effects, slow dynamics, and super-Arrhenius relaxation have also been reported in two-dimensional spin systems with competing long-range and short-range interactions <cit.>.


On a different front, a thermodynamic theory, called the “frustration-limited domain theory”  of the properties of supercooled liquids, and of the extraordinary increase of their characteristic structural relaxation times as the temperature is lowered, was formulated in terms of the postulated existence of a narrowly avoided thermodynamic phase transition due to geometric frustration <cit.>
(see Ref. <cit.> for a review). In this context frustration describes an incompatibility between extension of the locally preferred order in a liquid and tiling of the whole space. 
This picture is consistent with appropriate minimal statistical mechanical models, such as three-dimensional Ising Coulomb frustrated lattice models, which display a slowing down of the relaxation in Monte Carlo simulations <cit.> and an ideal glass transition within mean-field approximations <cit.>.
However, numerical simulations of these models in 3d are limited by the presence of a first-order transition to a modulated, defect-ordered phase <cit.>, and cannot be performed at sufficiently low temperatures.

Several frustrated spin (or Potts) lattice models without quenched disorder have also been introduced and studied over the past years to describe the key features of the glass transition. However, most of them are either mean-field in nature (and cannot be easily generalized to finite dimensions) <cit.>,  
or are characterized by (unphysical) multi-body interactions <cit.>.
The classical three-coloring model on the two-dimensional hexagonal lattice have been shown to undergo a dynamical freezing in metastable states very similar to the one observed in structural glasses <cit.>.
Slow dynamics also appears in electronic Coulomb liquids on the triangular lattice at quarter-filling <cit.>, as well as in spin-ice systems both in 2d <cit.> and in 3d <cit.>.
On the quantum side, it was shown in Ref. <cit.> that a valence bond glass phase emerges in the SU(N) Hubbard-Heisenberg model on a Bethe lattice in the large-N limit due to the interplay of strong magnetic frustration and quantum fluctuations. 

Yet, despite all these efforts over the past years, a clear and coherent picture of the glassy behavior that can arise due to the effect of geometric frustration in finite dimensional magnetic systems at low temperature is still missing.

Recently Hamp & al. <cit.> studied an Ising model on the Kagome lattice with short range antiferromagnetic interactions and dipolar interactions decaying as 1/r^3—the Dipolar Kagome Ising Antiferromagnet (DKIAFM) introduced in <cit.>. 
By means of extensive Monte Carlo simulations the authors first showed evidence for a first-order transition from the high temperature paramagnetic phase to a low temperature  crystal state that breaks time-reversal and sublattice symmetries, and coincides with the one previously proposed in Ref. <cit.> as the ground state.
Furthermore, upon cooling below the first-order transition, the system enters a supercooled liquid regime which exhibits all the characteristic features of fragile glasses: 
two-time autocorrelation functions decay as stretched exponentials and the relaxation time grows in a super-Arrhenius fashion as the temperature is decreased.
However, these conclusions were drawn out of numerical simulations of relatively small systems (about 300 spins) and might be affected by both strong finite-size effects and the difficulty of reaching thermal equilibrium in a reliable fashion due to strong metastability effects.
Moreover, a consistent picture of the physical origin of the dynamical slowing down at low temperatures has not been convincingly established yet.

In order to overcome, at least partially, these issues, in this paper we perform an analytical study 
of the equilibrium phase diagram of the DKIAFM in the thermodynamic limit at a mean-field level, focusing both on the ordered state and the glassy phase. Our results essentially confirm, support, and elucidate the observations reported in Ref. <cit.>.
Upon decreasing temperature, we first find a transition to a six-fold degenerate crystal state which breaks time reversal and rotation sublattice symmetry as the one observed in <cit.>.
The mean-field analysis indicates that the transition is indeed discontinuous. However, its first-order nature turns out to be extremely weak: the spinodal point of the crystal phase is very close to the transition point, resulting in a very large jump of the specific heat at the transition.  This feature provides a possible explanation of the fact that the finite-size scaling of the numerical data of the maximum of the specific heat with the system size performed in Ref. <cit.> did not find the usual behavior (C_ max∝ N) expected at a first-order transition due to very large finite-size effects.

When the system is supercooled below the first-order transition, we find that the paramagnetic state becomes unstable below a temperature at which the spin glass susceptibility diverges. Here, a continuous spin glass transition takes place at the mean-field level <cit.>. 


Note that the fact that the model displays a continuous spin glass transition in mean-field instead of a Random First-Order Transition <cit.> 

of the kind found in structural glasses (such as hard spheres in infinite dimensions <cit.> and lattice glass models on the Bethe lattice <cit.>) is perhaps not surprising. In fact, Ising spins with antiferromagnetic couplings on high-dimensional frustrated lattices and other related frustrated mean-field models with pairwise interactions are known to undergo a continuous transition to a spin glass phase 
when the temperature is lowered below the critical temperature of the antiferromagnetic phase <cit.>.

Beyond the fact that both spin glasses and structural glasses exhibit a pronounced slowing down of the dynamics upon cooling and aging in the low temperature phase, several important qualitative and quantitative differences characterize the dynamical behavior of these systems: In structural glasses two-time autocorrelation functions generically exhibit a two-step relaxation, characterized by a relatively fast decay to a plateau (i.e., the Edwards-Anderson order parameter) which appears discontinuously upon lowering the temperature, followed by a much slower decay, described by a stretched exponential. Moreover, the structural relaxation time is found to grow extremely fast, in a super-Arrhenius fashion, as the temperature is decreased <cit.>. Conversely, in spin glasses the Edwards-Anderson order parameter is continuous at the transition and vanishes in the paramagnetic phase. Hence two-time autocorrelation functions should display a simple exponential decay when the transition is approached from the high temperature phase, and an algebraic decay at the critical point. Furthemore, the relaxation time is expected to diverge (only) as a power-law at the critical point (and to stay infinite in the whole low temperature phase) <cit.>.
Nonetheless these differences are not clearly visible in numerical simulation of relatively small samples. A clear example of that is provided by the analysis of the dynamics of 3d Ising spin glasses performed in Ref. <cit.> using Monte Carlo simulations of systems with up to 64^3 spins. The two-time autocorrelation function was found to be very well fitted by stretched exponentials, with an exponent β that exhibits a temperature dependence extremely similar to the one reported in <cit.> for the DKIAFM. Moreover, although the divergence of the relaxation time as a power law, τ∼ (T - T_c)^-z ν, is consistent with the numerics, a Vogel-Fulcher law, τ∼ e^E_0/(T-T_c), was also found to account reasonably well for the data.

The lower-critical dimension of the spin glass transition is expected to be d_L ≈ 2.5 <cit.> (at least in the case of short-range interactions). Hence on general grounds we do not expect a genuine spin glass phase for the DKIAFM in 2d.
Yet, the manifestations of the vestige of the transition can be very strong also in two dimensional systems: The spin glass amorphous order can establish over very long (although not infinite) length scales, the spin glass susceptibility can become very large (although not infinite), and the relaxation time can grow very fast at low temperature.
Several experimental realizations of two-dimensional spin glasses using thin films do indeed show the same behavior as 3d spin glasses at sufficiently low temperature <cit.>.
In this sense, the existence of the spin glass phase in higher dimension, accompanied by the growth of long-range amorphous order and a rough free-energy landscape, provides a possible and natural explanation of the slow dynamics observed in the numerical simulations of Ref. <cit.> of the 2d model at low temperatures.

Despite the fact that our mean-field approach consists in studying the model on a random sparse graph of triangular Kagome plaquettes, and cutting-off the dipolar interactions beyond the second nearest-neighbour plaquettes (i.e., the 5th nearest-neighbor spins), it provides a remarkably good approximation for the equilibrium properties of the 2d DKIAFM. For instance, the mean-field approach yields a zero temperature entropy density of the nearest-neighbor Kagome spin ice model obtained for D=0 <cit.> equal to s_ GS≈ 0.75204, which turns out to be extremely close to the Pauling estimate s_ GS≈ 0.75225. Similarly, the ground state energy density of the crystal ground state within the mean-field approximation is e_ GS≈ -1.6116, which accounts reasonably well for the one found in the Monte Carlo simulations of systems with 300 spins, e_ GS≈ -1.515 <cit.>. As expected, the transition temperature to the crystalline phase is overestimated (by about a factor 3) by the mean-field treatment. Yet, the temperature dependence of the specific heat, the energy density, and the magnetization are remarkably similar, also at a quantitative level, to the ones found with Monte Carlo simulations (see Fig. <ref>).

The results presented here can serve at least two purposes:
(i) They help to support, understand, and clarify the numerical results of Ref. <cit.>.
(ii) They provide a first step to bridge the gap between the slow dynamics observed in geometrically frustrated magnetic systems and the theory of the glass transition formulated in terms of rough free-energy landscapes.
We believe that this analysis is of particular interest, especially in the light of the experimental relevance of the model, which could be potentially realized in several realistic set ups, including colloidal crystals <cit.>, artificial nanomagnetic arrays <cit.>, cold polar molecules <cit.>,  atomic gases with large magnetic dipole moments <cit.>, and layered bulk Kagome materials <cit.>.

The paper is organized as follows. In the next section we introduce the model. In Sec. <ref> we describe the mean-field approach, based on a cluster formulation of the problem on the Bethe lattice. In Sec. <ref> we show the results found within our analytical treatment, including the phase diagram 
and the equation of state. Finally, in Sec. <ref> we provide some concluding remarks and perspectives for future work.



§ THE MODEL


We consider the DKIAFM <cit.>
in which N classical spins S_i = ± 1 are
placed on the vertices of a two-dimensional 
Kagome lattice and point in a direction
perpendicular to the plane.
The Hamiltonian comprises an antiferromagnetic exchange term of
strength J between spins at nearest-neighbor lattice sites
⟨ ij ⟩ and long-range dipolar interactions of characteristic
strength D between all pairs of spins:

    H = J ∑_⟨ ij ⟩ S_i S_j + D/2∑_i ≠ jS_i S_j/r_ij^3  ,

where the distance r_ij = |  r_i -  r_j|/a between the spins i and j is measured in units of the lattice spacing a (that we set equal to 1 throughout).

In the following we will be interested in the case in which both interactions are antiferromagnetic, i.e., J>0 and D>0.
The case D = 0 is known to be fully frustrated and
does not order down to zero temperature <cit.>. 

The phase
diagram of the J = 0 model is less well understood
but the system is again strongly frustrated with any ordering (if present at all) suppressed down to temperatures
T ≪ D <cit.>.

The previous studies of the model <cit.> considered  the  coupling  parameters D= 1   K and J = 0.5   K (setting k_B= 1 and measuring all energies in Kelvin). A further advantage of developing an analytic (although approximate) treatment is that it is relatively simple to explore the parameter space. Without loss of generality we set J=0.5   K throughout (as in <cit.>) and study the phase diagram of the model and the constitutive equations in the different phases varying the dipolar coupling D and the temperature T.





§ THE MEAN-FIELD ANALYSIS
 

Our mean-field treatment is based on the combination of the cluster variational Bethe-Peierls
formalism (already successfully employed in the study of the equilibrium properties of geometrically frustrated magnetic systems <cit.>) and the cavity method <cit.>, developed in the context of glassy and disordered systems described by replica symmetry breaking (RSB). The latter concept is related to a complex free-energy landscape with special structure and the calculational meaning of it, in the context of the cavity method, will become clear below. 
In particular, we define the model on a Random Regular Graph <cit.> (RRG) of 
N_△ triangular plaquettes of
total coordination three (see Fig. <ref> for a sketch). In this case, the number of spins is equal to N=3 N_△/2 since since each spin belongs to two plaquettes and each plaquette contains three spins.
RRGs are a special class of sparse graphs, whose elements are chosen at random with uniform probability over the ensembles of all graphs of N_△ nodes, such that each node (i.e., a triangular plaquette) has exactly three neighbors.
RRGs have a local tree-like structure, which allows one to obtain exact self-consistent recursion relations for the probability
distributions of the spin configurations on each plaquette of the graph.
Yet, they have large loops, whose typical
length scales as log N_△ and diverges in the thermodynamic limit. Hence, a RRG 
is locally a tree, but it is frustrated, does not have a boundary, and is statistically translational invariant. For these reasons RRGs are
suitable lattices to study the thermodynamics of glassy and disordered systems at a mean-field level <cit.> (i.e., in the limit of infinite dimensions). 



 §.§ The “cavity” recursion relations
 

The standard way to obtain the recursion relations for the  marginal  probabilities of observing a given spin configuration on a given plaquette is provided by the cavity 
method <cit.>, which is equivalent to the Bethe-Peierls approximation at the replica-symmetric level. The cavity method is based on the assumption that, due to the tree-like structure of the lattice, in absence of a given plaquette (the cavity, e.g. the red triangle of Fig. <ref>), the neighboring plaquettes (the black triangles of Fig. <ref>) are uncorrelated and their marginal joint probabilities factorize.
Thanks to such factorization property one can write relatively simple recursion equations for the marginal probabilities of the cavity sites. Such equations have to be solved self-consistently, the fixed points of which yield the free-energy of the system along with all the thermodynamic observables (all the technical details of the method can be found in Refs. <cit.>). However, in order to be tractable, the cavity approach is formulated for systems with finite range interactions. Hence, before proceeding further we need to treat the dipolar interactions of Eq. (<ref>) in an approximate fashion. In practice, in the analytic calculations described below we choose to cut-off the dipolar couplings up to second nearest-neighboring plaquettes (i.e., the interactions between the spins belonging to the red plaquette α of Fig. <ref> and the spins belonging to green plaquettes are set to zero).

Consider now the cavity triangle α (red) in absence of one of its neighboring plaquettes β (dashed black). We define p_α→β ({S_α}|{S_β}) as the probability to observe the spin configuration {S_α}≡{ S_α,1,  S_α,2,  S_α,3} on the cavity triangle α of the (rooted) RRG, given that the spin configuration of the plaquette β is {S_β}≡{ S_β,1,  S_β,2,  S_β,3}.
We have adopted the convention that spin 1 is the root of the cavity plaquette and spins 2 and 3 are labeled anticlockwise. Using this convention one has that S_β,1≡ S_α,1.

The probabilities p_α→β ({S_α}|{S_β}) can be written in terms of the marginal probabilities defined on the cavity triangles γ and δ in absence of the triangle α, times the Gibbs' weight associated to each spin configuration:


    p_α→β({S_α}|{S_β})    = (  Z_α→β^( iter))^-1∑_γ, δ_α p_γ→α({S_γ}|{S_α})   p_δ→α({S_δ}|{S_α}) 
        
        
          e^-βH̃_α→β  ({S_α
        , 
        S_γ
        , 
        S_δ}
        | {
        S_β}) 
          ,

where Z_α→β^( iter) is a normalization factor ensuring that ∑_{S_α},{S_β} p_α→β({S_α}|{S_β}) = 1
and is associated to the “free-energy shift” involved in the iteration process:
- βΔ F_α→β^( iter)≡log Z_α→β^( iter).
Here we introduce the notation ∑_γ, δ_α that indicates the sum over all possible configurations {S_γ} and {S_δ} of the spin degrees of freedom of the plaquettes γ and δ, compatible with the constraints imposed by the spin configuration {S_α} on the plaquette α, i.e., S_γ,1 = S_α,2 and S_δ,1 = S_α,3 [see Eq. (<ref>)].
We will use this notation throughout this section.
The Hamiltonian 
H̃_α→β ({S_α
,
S_γ
,
S_δ} |
{
S_β}) appearing in the Gibbs factor of Eq. (<ref>) is a modified Hamiltonian, Eq. (<ref>), restricted to the cavity plaquette α→β:

    H̃_α→β ({S_α ,S_γ,S_δ} | { S_β})
       ≡ H_α^( AF) + H̃_γ→α^(2 D) + H̃_δ→α^(3 D) + H̃_(γ,δ)^(2 NND) + H̃_(γ,β)^(2 NND) + H̃_(δ,β)^(2 NND)  , 
     H_α^( AF)    = ( J + D  ) [ S_α,1 S_α,2 + S_α,1 S_α,3 + S_α,2 S_α,3]   , 
    H̃_γ→α^(2 D)    = D  [  S_α,1 S_γ,2 + S_α,3 S_γ,3/3 √(3) + S_α,1 S_γ,3 + S_α,3 S_γ,2/8]   , 
    H̃_δ→α^(3 D)    = D  [  S_α,1 S_δ,3 + S_α,2 S_δ,2/3 √(3) + S_α,1 S_δ,2 + S_α,2 S_δ,3/8]   , 
    H̃_(γ,δ)^(2 NND)    = D  [ S_γ,3 S_δ,2/8 +  S_γ,2 S_δ,2 + S_γ,3 S_δ,3/7 √(7) + S_γ,2 S_δ,3/27 ]   .

The meaning of this decomposition is the following. H̃_γ→α^(2 D) contains the 4 dipolar interaction terms between the spins belonging to the plaquette α and the spin belonging to its nearest-neighbor plaquette γ attached to the spin S_α,2, which are not already contained in H_α^( AF). H̃_(γ,δ)^(2 NND) contains the 4 dipolar interaction terms between the spins of the second nearest-neighbor plaquettes γ and δ which are not already contained in H̃_γ→α^(2 D) and H̃_δ→α^(3 D).

In order to obtain the marginal probabilities of the spin configurations on each plaquette of the (unrooted) RRG (where each triangular plaquette has exactly three neighbors), one needs to merge three cavity plaquettes (e.g., plaquettes β, γ, and δ of Fig. <ref>) onto their neighboring plaquette (e.g., plaquette α of Fig. <ref>). In this way one obtains:

    P_α({S_α})    = (  Z_α^(s))^-1∑_β, γ, δ_α p_β→α({S_β}|{S_α})   
        p_γ→α({S_γ}|{S_α})   p_δ→α({S_δ}|{S_α}) 
          e^-βH̃_α ({S_α
        ,
        S_β
        ,
        S_γ
        ,
        S_δ})
          ,

where Z_α^(s) is a normalization factor ensuring that ∑_{S_α} P_α({S_α}) = 1, and is associated to the “free-energy shift” involved in the process of joining three cavity plaquettes (β, γ, and δ) to a central plaquette (α): - βΔ F_α^(s)≡log Z_α^(s).
The plaquette Hamiltonian H̃_α ({S_α,S_β,S_γ,S_δ}) reads

    H̃_α ({S_α ,S_β, S_γ,S_δ})
       ≡ H_α^( AF) + H̃_β→α^(1 D) +
    H̃_γ→α^(2 D) + H̃_δ→α^(3 D) + H̃_(γ,δ)^(2 NND) + H̃_(γ,β)^(2 NND) + H̃_(δ,β)^(2 NND)  , 
    H̃_β→α^(1 D)    = D  [  S_α,2 S_β,3 + S_α,3 S_β,2/3 √(3) + S_α,2 S_β,2 + S_α,3 S_β,3/8]   ,

and the other terms are given in Eq. (<ref>).

The equilibrium averages of all local observables which involve the spin degrees of freedom of a given plaquette, including, e.g.,  the magnetization, can be expressed in terms of these marginal probabilities: 

    ⟨ O_α⟩ = 
    ∑_{S_α} O ({S_α}) P_α({S_α})   .

Similarly, the contribution to the average energy due to the plaquette α can be expressed as

    ⟨ e_α^(s)⟩    = ∑_α→ (β,γ,δ)   
        p_β→α({S_β}|{S_α})   
        p_γ→α({S_γ}|{S_α})   p_δ→α({S_δ}|{S_α}) 
         H̃_α ({S_α,S_β,S_γ,S_δ})  
        e^-βH̃_α ({S_α
        ,
        S_β
        ,
        S_γ
        ,
        S_δ})/∑_α→ ( β, γ, δ )  p_β→α({S_β}|{S_α})   
        p_γ→α({S_γ}|{S_α})   p_δ→α({S_δ}|{S_α}) 
          e^-βH̃_α ({S_α
        ,
        S_β
        ,
        S_γ
        ,
        S_δ})  .

The process of joining two neighboring cavity plaquettes (e.g., plaquettes α and β of Fig. <ref>) involves another “free-energy shift”, defined as

    e^- βΔ F_α↔β^(l)≡ Z_α↔β^(l) = ∑_α↔β
     
        p_α→β({S_α}|{S_β})   
        p_β→α({S_β}|{S_α})   
        e^-βH̃_β→α^(1 D)  ,

where the Gibbs' factor H̃_β→α^( D) has been defined in Eq. (<ref>). Similarly, the contribution to the average energy coming from the interactions between two neighboring cavity plaquettes is given by

    ⟨ e_α↔β^(l)⟩ = ∑_α↔β
        p_α→β({S_α}|{S_β})   
        p_β→α({S_β}|{S_α})  H̃_β→α^(1 D) 
        e^-βH̃_β→α^(1 D)/∑_α↔β
        p_α→β({S_α}|{S_β})   
        p_β→α({S_β}|{S_α})   
        e^-βH̃_β→α^(1 D)  .

We recall here the convention adopted for the notation of the summation over the spin degrees of freedom in the expressions above:

    ∑_γ, δ_α   ≡∑_{S_γ},{S_δ}
      S_γ,1 = S_α,2
     S_δ,1 = S_α,3  ,                                 ∑_β, γ, δ_α≡∑_{S_β}, {S_γ},{S_δ}
     
        S_β,1 = S_α,1
    
        S_γ,1 = S_α,2
     S_δ,1 = S_α,3  , 
    
       ∑_α→ (β, γ, δ)    ≡∑_{S_α},{S_β}, {S_γ},{S_δ}
     
        S_β,1 = S_α,1
    
        S_γ,1 = S_α,2
     S_δ,1 = S_α,3 
        = ∑_{S_α}∑_β, γ, δ_α  ,                   ∑_α↔β≡∑_{S_α},{S_β}
     
        S_β,1 = S_α,1  .

The free-energy 
of the system can be obtained by combining the free-energy shifts involved in the different processes, as explained in <cit.>:

    F = ∑_α=1^N_△Δ F_α^(s)
    - ∑_⟨α, β⟩Δ F_α↔β^(l)
    = 1/2∑_⟨α, β⟩( 
    Δ F_α→β^( iter) + Δ F_β→α^( iter)) - 1/2∑_α=1^N_△Δ F_α^(s)  ,

where 
⟨α , β⟩ denotes the sum over the  3 N_△ /2 nearest-neighbors plaquettes on the graph (the last equality simply comes from the fact that Δ F_α^(s) = Δ F_α→β^( iter) + Δ F_α↔β^(l) by construction). The average entropy of the system is then given by ⟨ S ⟩ = β (⟨ E ⟩ - F).
Analogously, the total average energy can be written as

    ⟨ E ⟩ = ∑_α=1^N_△⟨ e _α^(s)⟩ 
    - ∑_⟨α, β⟩⟨ e_α↔β^(l)⟩  .



Equations (<ref>) can be written for arbitrary (large) RRGs and are expected to become exact in the thermodynamic limit.
On each triangle of the RRG one can define three cavity plaquettes by removing one of its three neighbors. Thus, Eqs. (<ref>) represent a set of 32 × 3 × N_△ coupled nonlinear algebraic equations for the 32 marginal probabilities p_α→β({S_α}|{S_β}) associated to the 32 possible configurations of the spins {S_α} on the cavity plaquette α, given the configuration of the spins {S_β}. Once the fixed points of these equations is found, one can compute the marginal probabilities on each plaquette of the graph from Eq. (<ref>), along with the free-energy and all observables.
In the following, we will discuss three specific solutions of the equations in the thermodynamic limit, corresponding to the (RS) homogeneous paramagnet, the (RS) ordered crystalline state, and the (RSB) glassy phase.



 §.§ The paramagnetic phase


The paramagnetic phase is characterized by translational invariance and corresponds to the homogeneous and RS solution of the recursion relations:

    p_α→β({S_α}|{S_β})  =  p_ para({S_α}|{S_β})      ∀α,β   .

 The probabilities p_ para({S_α}|{S_β}) are given by the fixed point of Eqs. (<ref>) which in this limit become a simple system of 32 coupled nonlinear algebraic equations. The free-energy, the energy, and the magnetization (which is identically zero by Z_2 inversion symmetry in the paramagnetic phase, which implies that p_ para({S_α}|{S_β}) = p_ para({-S_α}|{-S_β})), can be easily computed from Eqs. (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>). 
This phase is expected to be stable at high temperature. However, 
the average entropy density ⟨ s ⟩ = β( ⟨ e ⟩ - f) becomes  negative  below  a  certain temperature, T_s=0(D).  This indicates that the homogeneous solution is certainly not appropriate to describe the low temperature region of the phase diagram.



 §.§ The stability of the paramagnetic phase



The manifestation of the failure of the RS solution also shows up via a loss of stability of the RS fixed point, as given by a simple linear analysis.
To describe this instability one needs  to  introduce  a  probability distribution  P [ p⃗
], where p⃗ is a short-hand notation for the 32 marginal probabilities p({S_α}|{S_β}) and P [ p⃗] is defined as the probability density that the 
probabilities p_α→β({S_α}|{S_β}) on the cavity plaquette α 
are equal to p({S_α}|{S_β}).

In the homogeneous phase from Eq. (<ref>) one has that  the probability distributions of the marginal probabilities on the triangular plaquettes must satisfy the following self-consistent equation:

    P [ p⃗ ] = ∫ d P [ p⃗_γ ]   d P [ p⃗_δ ]  δ [ p⃗ - 𝚙⃗ ( p⃗_γ, p⃗_δ )]   ,

where 𝚙⃗ ( p⃗_γ, p⃗_δ ) is a short-hand notation for the r.h.s. term of the recursion relations (<ref>). Close to the homogeneous paramagnetic solution we have, to first order,

    p ({S_α}|{S_β})  ≈  p_ para({S_α}|{S_β}) + δ p  ({S_α}|{S_β})    .

Starting with δp⃗ identically and independently distributed and injecting the expression above into Eq. (<ref>), one has that the deviation of the marginal probabilities from the homogeneous solution evolves under iteration as

    ⟨δp⃗⟩ = 2  .∂ 𝚙⃗ ( p⃗_γ, p⃗_δ )/∂p⃗_γ|_ para⟨δp⃗⟩  .

where ⟨·⟩ refers to the average using the distribution P(p⃗). ∂ 𝚙⃗ / ∂p⃗_γ is actually a 32 × 32 Jacobian matrix. If λ_ max denotes the eigenvalue of largest modulus of that matrix, the stability criterion simply reads 2 |λ_ max| ≤ 1.
When 2 |λ_ max| > 1, the paramagnetic solution is instead unstable with respect to a “modulation” instability, corresponding to a transition to a regime with successive (homogeneous) generations of the tree carrying different values of the marginal probabilities. Such modulation instability is thus a manifestation of an instability toward an ordered phase, which breaks translational invariance.

This instability criterion can also be obtained by studying response functions to a perturbation (which is related to correlations through the fluctuation-dissipation theorem) <cit.>. In this setting the instability is detected by means of the divergence of the linear magnetic susceptibility in the paramagnetic phase, defined as:

    χ = 1/N_△∑_α,β⟨ S_α S_β⟩_c = 1/N_△∑_α,β.
    ∂⟨ S_α⟩/∂ h_β|_h_γ=0  ,

where S_α≡∑_i ∈α S_i is a short-hand notation for the magnetization of the plaquette α and h_β is an external magnetic field conjugated to the magnetization of the plaquette β.
Making use of the homogeneity of the paramagnetic solution and the tree-like structure of the lattice, the susceptibility can be rewritten as

    χ = 1 + 3 ∑_r=1^∞ 2^r-1⟨ S_α_0 S_α_r⟩_c   .

where S_α_0 and S_α_r are two plaquettes taken at distance r on the tree. The series converges provided that lim_r →∞log⟨ S_α_0 S_α_r⟩_c/r < log 2. To evaluate ⟨ S_α S_β⟩_c, we invoke the fluctuation-dissipation relation:

    ⟨ S_α_0 S_α_r⟩_c = 
    .
    ∂⟨ S_α_r⟩/∂ h_α_0|_h_γ=0

where h_α_0 denotes the external magnetic field conjugate to S_α_0. Since h_α_0 is a function of (the components of) P({S_α_0}), we can use the chain rule along the branch of the tree which connects the plaquette α_0 with the plaquette α_r through the plaquettes α_l, l=1, …, r-1:

    ∂⟨ S_α_r⟩/∂ h_α_0 = ∂⟨ S_α_r⟩/∂p⃗_α_r →α_r-1( ∏_l=2^r∂p⃗_α_l →α_l-1/∂p⃗_α_l-1→α_l-2) ∂p⃗_α_1→α_0/∂P⃗_α_0∂P⃗_α_0/∂ h_α_0  .

In the paramagnetic phase, all the intermediate marginal cavity probabilities are equal and the previous equation factorizes, leading again to 2 |λ_ max| ≤ 1.


The maximal eigenvalue λ_ max increases as the temperature is lowered and the linear susceptibility of the paramagnetic phase diverges at a certain temperature, signaling a modulation instability of the paramagnetic phase toward a crystalline phase (see Sec. <ref>) at a temperature T_ mod(D).

One can also look for another kind of instability, namely a spin glass  instability, which manifests itself as a divergence of the non-linear susceptibility <cit.>, which is defined as

    χ_ sg = 1/N_△∑_α,β⟨ S_α S_β⟩_c^2   .

Equivalently, this instability appears as a widening of the variance ⟨ ( δp⃗ )^2 ⟩ under the recursion of Eq. (<ref>). Both approaches lead to a stability criterion 2 λ_ max^2 ≤ 1. 
Note that this condition is always weaker than that for the modulation instability, 2 |λ_ max| ≤ 1, associated to the crystalline order. However, it is the relevant one in the case of glassy phases, characterized by the establishment of long-range amorphous order.

Solving the recursion relations (<ref>) in the paramagnetic phase, we find that the homogeneous solution becomes unstable below a temperature T_ sg (D), at which the spin glass susceptibility diverges (with T_s=0(D) < T_ sg(D) < T_ mod (D)).

This requires either a phase transition before the spin glass local instability is reached <cit.> (as occurs in the mean-field models of fragile glasses, described by a Random First-Order Transition <cit.>), or a continuous (possibly spin glass) transition at T_ sg. We will show below that the latter scenario is the correct one for the DKIAFM.
In order to do this 
in Sec. <ref> we look for a solution of the recursion relations which breaks the replica symmetry, corresponding to a glassy phase where many local minima of the free-energy exist and where the local marginal probabilities fluctuate from a plaquette to another.





 §.§ The crystal phase
 

One can look for a crystalline RS solution, where the marginal probabilities do not fluctuate from  site  to  site,  but  are  different  in  different  sites (break-down of translational invariance). 
The (sixfold degenerate) crystalline state proposed in <cit.> and observed numerically in <cit.> is characterized by a 12-spin unit cell and breaks (twofold) time-reversal symmetry and (threefold) rotation symmetry (see Refs. <cit.> for more details). In order to be able to account for such ordered phase we need to introduce 8 sublattices of triangular plaquettes (see Fig. <ref>), corresponding to 24 different sets of cavity plaquettes. 
The  merging of the 24 cavity plaquettes is  done  taking  into  account  the structure  of  the  crystalline  phase, as explained in the caption of Fig. <ref>. The recursion equations (<ref>) become then a set of 24 × 32 coupled nonlinear algebraic equations for the marginal probabilities on the 24 cavity plaquettes on each sublattice. The solution of these equations appears  discontinuously at  a spinodal  point T_ sp(D), and becomes  thermodynamically  stable when the corresponding free-energy crosses the paramagnetic one, at the melting temperature T_m(D). At that temperature we observe a first-order phase transition characterized by a spontaneous breakdown of the translational, rotational, and spin inversion invariance, accompanied by a discontinuous  jump  of  the energy density and of the entropy density. Decreasing further the temperature, the energy in the crystalline  phase  approaches  quickly  the ground state value, e_ GS≈ -1.6116 (which turns out to be remarkably close to the one found with Monte Carlo simulations of systems of 300 spins, e_ GS≈ -1.515  <cit.>), and  the  entropy quickly approaches zero.

Inspecting the (ground state) spin configuration of Fig. <ref>, it was noticed in <cit.> that one of the three spins of the Kagome triangles are completely polarized (i.e., the bottom spins of sublattices a,b,c,d and the top spins of sublatices A,B,C,D), 
with the state having zero magnetization overall. 
Note that the need to introduce 8 sublattices of triangular plaquettes is due to the fact that the spin pattern on the two non-polarized rows of spins of the Kagome triangles (i.e., along the horizontal bonds in Fig. <ref>) has period four, with three spins S = ∓ 1 followed by one spin S=± 1.
Based on these observations, suitable order parameters for the transition to the ordered state are the sublattice magnetizations:

    m_X = ∑_i ∈ X S_i   ,

with X denoting the 8 different sublattices: X={A,B,C,D,a,b,c,d}. 
These order parameters essentially coincide with the emergent effective charge variables introduced in Ref. <cit.>, derived from the so-called dumbbell picture <cit.>.  The (ground state) spin configuration of Fig. <ref> corresponds to m_a=m_c=m_A=m_D= ± 1 and m_b=m_d=m_B=m_C= ∓ 1. 
Equivalently, one can choose as order parameter m_ pol =  ⟨ S_ pol⟩, the average magnetization of the spins of the Kagome triangles that are completely polarized, as done in <cit.>.
Following this suggestion, we will use |m_ pol| as the order parameter jumping from zero to a finite value at the transition.

We have also looked for other plausible competing ordered phases, which break the translational and rotational symmetries in different ways and have a different unit cells. 
However, such alternative crystalline states turn out to be less favourable (i.e., they have a higher free-energy) compared to the crystalline phase of <cit.>. Yet, if one does not include the dipolar interactions between the second nearest-neighboring plaquettes, the crystalline phase depicted in Fig. <ref> disappears (i.e., no physically relevant fixed point of the recursion relations is found corresponding to the sublattice structure of Fig. <ref>), and another completely different (fourfold degenerate) ordered phase emerges.
This observation highlights the importance of accounting for the dipolar interaction as accurately as possible, in order to recover the correct description of the ordered phase <cit.>.



 §.§ The spin glass phase
 

The  paramagnetic  phase  is  metastable  below T_m (D), corresponding to a supercooled regime. However, as mentioned 
above, the predicted entropy density becomes negative as the temperature is lowered below T_s=0(D), implying that this solution does not describe well the low temperature region. 
Moreover, the homogeneous solutions becomes unstable below a certain temperature T_ sg (D)>T_s=0(D), at which the spin glass susceptibility diverges. 
The “entropy crisis” and the spin glass instability are manifestations of the appearence of a huge number of metastable glassy states. The RS approach  fails  because  it  does  not  take  into  account  the  existence of several local minima of the free-energy.
This requires either a phase transition before the spin glass local instability is reached (as in the case of lattice models for fragile glasses in the mean-field limit <cit.> described by a Random First-Order Transition <cit.>), or a continuous spin glass transition at T_ sg.
In order to understand which of these two possible scenarios is the correct one for the DKIAFM, we have to look for a solution of the recursion relations which breaks the replica symmetry, corresponding to a glassy phase where many local minima of the free-energy exist and where the local marginal probabilities fluctuate from a plaquette to another.
We thus need to perform a statistical treatment of sets of solutions of Eq. (<ref>). The simplest setting which allows to proceed further in this direction is provided by a one-step RSB ansatz, which starts from the assumption that exponentially many (in N_△) solutions of the recursion relations exist. 
More precisely, we assume that the number N (f) of solutions with a given free-energy density f on graphs of size N_△ is N (f) ∼exp [ N_△Σ (f)], where Σ(f) ≥ 0 is called the configurational entropy (or complexity) and is supposed to be an increasing and concave function of the free-energy f. This is a strong hypothesis which is justified by its self-consistency. Under these assumptions, one can show that the 1RSB self-consistent equation for the probability distribution of the marginal cavity probabilities becomes <cit.>

    P_m [ p⃗ ] ∝∫ d P_m [ p⃗_γ ]   d P_m [ p⃗_δ ]  δ [ p⃗ - 𝚙⃗ ( p⃗_γ, p⃗_δ )]
      e^-β m Δ F^( iter)  ,

where 𝚙⃗ ( p⃗_γ, p⃗_δ ) is a short-hand notation for the r.h.s. term of the recursion relations (<ref>) and Δ F^( iter) is the free-energy shift involved in the iteration process defined in Eq. (<ref>) via the normalization of the cavity marginal probabilities.
The probability distribution depends on the parameter m which is the  breakpoint  in  Parisi's order parameter function at the 1RSB level <cit.>, and is defined as m = (1/β) ∂Σ / ∂ f (all the details of the calculation can be found in Refs. <cit.>).
Similarly to Eq. (<ref>), the 1RSB free-energy density functional is given by

    ϕ(m) = Δϕ^(s)(m) - 3/2Δϕ^(l) (m)   ,

with

    e^-β m Δϕ^(s)    = ∫ d P_m [ p⃗_β ]   d P_m [ p⃗_γ ]   d P_m [ p⃗_δ ]   e^-β m Δ F^(s)  , 
    
    e^-β m Δϕ^(l)    = ∫ d P_m [ p⃗_α ]   d P_m [ p⃗_β ]   e^-β m Δ F^(l)  ,

where the free-energy shifts have been defined in Sec. <ref>.
The other relevant thermodynamic observables, such as, e.g., the average energy, can be obtained in a similar fashion <cit.>.
The parameter m is fixed by the maximization of the free-energy functional with respect to it <cit.>, which allows to recover the complexity as a Legendre transform of ϕ(m):

    m  ϕ (m) = m f - 1/βΣ (f)   .

The RS high-temperature homogeneous description of the phase is recovered by taking 
P_m(p⃗) = δ (p⃗ - p⃗_ para) and m=1 <cit.>. 

Since Eq. (<ref>) is a functional relation, an analytical treatment is not possible in general. Yet the self-consistent equation can be efficiently solved numerically with arbitrary precision using a population dynamics algorithm (for all technical details see <cit.>). For high values of the temperature (T>T_ sg (D)) we recover the paramagnetic solution. Lowering the temperature, a nontrivial solution of the 1RSB equation   appears continuously exactly at T_ sg (D). Right below T_ sg the probability distribution P_m (p⃗) acquires an infinitesimal widening of the variance ⟨ (δp⃗)^2 ⟩. 
This scenario corresponds to a continuous transition to a spin glass phase at the temperature at which the spin glass susceptibility diverges. 
The order parameter of the spin glass transition is the Edwards-Anderson order parameter, q_ EA = (1/N) ∑_i ⟨ S_i ⟩^2, which vanishes linearly as q_ EA∼ (T_ sg - T) for T → T_ sg^- <cit.>.

As it is well-known, the low-temperature spin glass 
phase should be described by full RSB <cit.>. However, any new level of RSB will require considering a more sophisticated situation, namely a distribution over the probability distribution of the previous level. For instance, the two-step RSB will be written as a distribution Q[ P[p⃗]] over distributions P(p⃗). Describing with this formalism a finite connectivity system with full RSB is therefore too complicated, and we will limit ourselves to the 1RSB Ansatz. 
Moreover, since solving the  self-consistent functional equation (<ref>) via  population  dynamics  is  quite computationally demanding, we did not perform the maximization of the free-energy functional (<ref>) with respect to m. For these reasons, our approach only provides an approximate description of the equilibrium properties of the spin glass phase and we have not pushed the 1RSB calculations far below T_ sg (essentially we only consider few values of the temperature in the vicinity of the critical point).





§ PHASE DIAGRAM AND THERMODYNAMIC BEHAVIOR
 

In this Section we discuss the main results found within the mean-field treatment of the DKIAFM described in the previous sections. In order to compare with the numerical results of the Monte Carlo simulations of <cit.>, we start by fixing the parameter D to 1, as in Refs. <cit.>, and measure several observables such as the average energy density ⟨ e ⟩ = ⟨ E ⟩/N_△, the (intensive) specific heat c = ∂⟨ e ⟩ / ∂ T, the magnetization of the polarized spin in one of the sixfold degenerate ground state configurations |m_ pol|, and the average entropy density ⟨ s ⟩ = ⟨ S ⟩/N_△, as a function of the temperature T in the paramagnetic, crystal, and 1RSB glass solutions of the recursive cavity equations. The results are shown in Fig. <ref>.
At high temperature the system is found in the paramagnetic phase. Upon lowering the temperature, a first-order transition to the crystalline phase proposed in Refs. <cit.> (see Fig. <ref>) occurs at T_m. The order parameter |m_ pol| presents a finite jump at T_m, where the average energy and entropy densities also display an abrupt decrease.
The transition to the ordered state turns out to be very weakly first-order, in the sense that the spinodal point of the crystalline solution, T_ sp≈ 0.16   K, is very close to the transition temperature T_m ≈ 0.1566   K where the free-energies of the paramagnetic phase and the crystal phase cross. These temperatures are also numerically close to the modulation instability temperature T_ mod≈ 0.147   K. Since the specific heat of the crystal solution diverges at the spinodal point, the vicinity of T_m and T_ sp results in a very large jump (of about a factor 3) of the intensive specific heat at the transition. This feature might explain the deviations observed in the numerical simulations of the expected scaling of the peak of the (extensive) specific heat as C_ max∝ N_△ <cit.>.


Although approximate, our approach accounts remarkably well for the numerical results of Ref. <cit.>. As expected, the transition temperature is overestimated by the mean-field approximation (by about a factor 3). Yet, the temperature dependencies of the specific heat, the energy, and the magnetization are, also at a quantitative level, very similar to the ones found in Ref. <cit.> (recall that the energy, entropy, and specific heat per spin are obtained by multiplying the energy, entropy, and specific heat per plaquette by a factor 2/3).

The  paramagnetic  phase  is  metastable  below T_m, corresponding to a supercooled regime.
If one keeps lowering the temperature within the supercooled phase, the spin glass susceptibility grows and diverges at T_ sg≈ 0.1   K, where a continuous transition to a spin glass phase takes place. Although the spin glass phase is presumably described by full RSB (at least at the mean-field level), our approach only allows one to perform an approximate 1RSB Ansatz for the low temperature glassy phase. Moreover, solving the self-consistent functional equation (<ref>) via population dynamics is computationally heavy. For these reasons, we did not push the calculations of the thermodynamic observables too deep into the spin glass phase, and only solved the equations for few points close to the critical temperature.



In Fig. <ref> we plot the phase diagram of the DKIAFM, showing the position of the different phases when varying the temperature T and the dipolar interaction D (J is fixed to J=0.5   K).
The effect of varying the dipolar interaction turns out to be particularly simple. In fact, we find that the phase boundaries, as well as all the characteristic temperature scales, vary linearly with D:

    T_m    ≈ 0.1566   D   ,
    
    T_ sg   ≈ 0.1   D   , 
    
    T_ sp   ≈ 0.16   D   , 
    
    T_ mod   ≈ 0.147   D   , 
    
    T_ s=0   ≈ 0.0713   D   .

As expected, in the limit D=0 the paramagnetic phase is stable at all temperatures and corresponds to the only solution of the recursion relations. This is due to the fact that for D=0 the system is much less frustrated and has a highly (i.e., extensively) degenerate ground state (i.e., ⟨ s ⟩ approaches a finite value in the T → 0 limit), since each plaquette has a sixfold degenerate ground state which corresponds to the ice rule (two +1 and one -1 spins or two -1 and one +1 spins). In particular, for D=0 the model reduces to the nearest-neighbor Kagome spin ice model of Wills, Ballou, and Lacroix <cit.>, for which a Pauling estimate yields the entropy s_ GS = (3/2) log[2 (3/4)^2/3] ≈ 0.75225, while our mean-field approximation yields s_ GS≈ 0.75204. 
When the dipolar interactions are turned on (D>0), such degeneracy is lifted, and a specific crystalline ground state structure emerges. The minimization of the local interactions produces a much stronger geometric frustration for T ≲ D and gives rise to the emergence of a spin glass phase at low temperatures, characterized by an extremely rough free-energy landscape (at least at the mean-field level).
The fact that all the relevant temperature scales of the problem show an apparent linear dependence on D is precisely due to the fact that the relevant energy scale is the energy difference between the ground state and the first excited states, which goes linearly to zero with D.



§ CONCLUSIONS
 

In this paper we have developed an analytical mean-field treatment for the equilibrium properties of the DKIAFM introduced in <cit.> and studied numerically in <cit.>. 
Our mean-field approach is based  on  a  cluster  variational Bethe-Peierls formalism <cit.> and on the cavity method <cit.>, and consists in studying the model on a sparse random tree-like graph of triangular Kagome plaquettes, and cutting-off the dipolar interaction beyond the second nearest-neighbor plaquettes (i.e., the 5th nearest-neighbor spins).
Our results essentially confirm and support the observations reported in Ref. <cit.>,
which were obtained using Monte Carlo simulations of relatively small system (ranging from 48 to 300 spins), and  might  be  affected  by  both  strong  finite-size  effects
and by the difficulty of reaching thermal equilibrium in a reliable fashion due to strong metastability effects.

The summary of our results is the following.
Upon decreasing the temperature we first find a transition to a sixfold degenerate crystal state which breaks time reversal, translation,  and  rotation  symmetry  as  the
one proposed in <cit.>. Such transition is indeed discontinuous, as suggested in <cit.>, although its first-order character turns out to be extremely weak, which might explain the strong finite-size effects observed in the finite-size scaling of the numerical data of the specific heat.
When the system is supercooled below the first-order
transition, we find that the paramagnetic state becomes
unstable below a temperature at which the spin glass susceptibility diverges and a continuous spin glass transition takes place at the mean-field level. 

On the one hand, the results presented here support and clarify
the numerical findings of Ref. <cit.>. On the other hand, they provide a first step to  bridge  the  gap  between  the  slow  dynamics  observed in geometrically frustrated magnetic systems and the mean-field theory of glassy systems formulated in terms of
rough free-energy landscape.  We believe that this analysis is of particular interest, especially in the light of the experimental relevance of the model, which could be potentially  realized  in  several  realistic  set  ups,  including
colloidal crystals <cit.>, artificial nanomagnetic arrays <cit.>, polar molecules <cit.>, atomic gases with large magnetic dipole moments <cit.>, and layered bulk Kagome materials <cit.>.

Some comments are now in order.

The lower-critical dimension of the spin glass transition 
is expected  to  be  
d_L ≈ 2.5 <cit.> (at  least  in  the  case
of  short-range  interactions).  Hence  we do not expect a genuine spin glass phase for the DKIAFM. Yet in 2d the spin glass amorphous order can establish over very large (although not infinite) length scales and the spin glass susceptibility can become very large
(although  not  infinite)  at  low  temperature  due  to  the
vestiges of the transition.  Indeed, there are plenty of experimental studies using thin films that at sufficiently low temperatures 
behave as the 3d counterparts. See, e.g. <cit.> for a 
very recent reference and <cit.>  
for a more classical ones. The situation is similar concerning numerical simulations <cit.>.

Concerning the dynamics, very early Monte Carlo simulations of the 3d Edwards-Anderson model suggested that the spin auto-correlation function, close but above the expected critical temperature, decays as a stretched exponential <cit.>.
Therefore, although the model is expected to have a conventional second order phase transition with critical slowing down and algebraic decay of correlation functions, 
for the system sizes and time-scales 
accessed in this paper, the time-delayed correlations 
were satisfactorily fitted by such an anomalous
form, with a stretching exponent decaying with decreasing temperature.
Just a bit later, in <cit.> the conventional critical slowing down was recovered. 
A stretched exponential relaxation
of the self-correlation in the DKIAFM was reported in <cit.>.
However, our results suggest that at sufficiently large time and length scales this behavior might be replaced by a conventional power law decay for control parameters in the critical region.

The case J=0 and D>0 <cit.> has been left over by the present investigation and might be an interesting subject for future studies. Possibly, the most interesting questions would be the  investigation of  how the properties of the model are  affected  by  quantum  fluctuations.


We would like to thank C. Castelnovo for enlightening and helpful discussions.
L. F. Cugliandolo and M. Tarzia are members of the Institut Universitaire de France.
This work is supported by “Investissements d'Avenir" LabEx PALM (ANR-10-LABX-0039-PALM) (EquiDystant project, L. Foini).


99

Toulouse G. Toulouse, Commun. Phys. 2, 115 (1977).

Balents2010 L. Balents, Nature 464, 199 (2010).

Moessner2006 R. Moessner and A. P. Ramírez, Phys. Today, No. 2, 59, 24 (2006). 

Jaubert2008 L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. 100, 067207 (2008).

Lieb1972 E. H. Lieb and F. Y. Wu, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic, New York, 1972), Vol. 1, Chap. 8, p. 331.

Youngblood1980 R. Youngblood, J. D. Axe, and B. M. McCoy, Phys. Rev. B 21, 5212 (1980).

Henley2010 C. Henley, Annu. Rev. Condens. Matter Phys. 1, 179 (2010).

Wolynes2012 P. G. Wolynes and V. Lubchenko, Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications (JohnWiley & Sons, 2012).

Berthier2011  L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).

Cavagna2009  A. Cavagna, Phys. Rep. 476, 51 (2009).

Gotze1991 W. Götze, Liquids, Freezing and the Glass Transition, edited by J. P. Hansen, D. Levesque, J. Zinn-Justin, Les Houches. Session LI, 1989 (North-Holland, Amsterdam, 1991).

Debenedetti2006 P.  G.  Debenedetti, Metastable  Liquids (Princeton  University Press,  Princeton,  1996).

Tarjus2011 G.  Tarjus,  in Dynamical  Heterogeneities  in  Glasses,  Colloids,  and  Granular  Media,  edited  by  L.  Berthier,  G.  Biroli, J.-P.  Bouchaud,  L.  Cipelletti,  and  W.  van  Saarloos  (Oxford University Press, New York, 2011), Chap. 2.

Marinari-etal94

E. Marinari, G. Parisi and F. Ritort, 
J. Phys. A: Math. Gen. 27, 7615 (1994).

ibid 7647 (1994).

BouchaudMezard

J.-P. Bouchaud and M. Mézard, 
J. Phys. I 4, 1109 (1994).

Cugliandolo-etal95

L. F. Cugliandolo, J. Kurchan, G. Parisi and F. Ritort, Phys. Rev. Lett. 74,  1012 (1995).

Bouchaud96

J.-P. Bouchaud, L. F. Cugliandolo, J. Kurchan, and M. Mézard,
Physica A 226,  243 (1996).

Grousson
M. Grousson, G. Tarjus, and P. Viot, Phys. Rev. E 65, 065103 (2002); 
J. Phys.: Condens. Matter 14, 1617 (2002).

Tarjus2005  G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot,

J. Phys.: Condens. Matter 17, R1143 (2005).

Westfahl2001 J. Schmalian and P. G. Wolynes, Phys. Rev. Lett. 85, 836 (2000); H. Westfahl, J. Schmalian, and P. G. Wolynes, Phys. Rev. B 64, 174203 (2001).

Biroli2002 G. Biroli and M. Mézard, Phys. Rev. Lett. 88, 025501 (2002). 

Ciamarra2003 M. P. Ciamarra, M. Tarzia, A. de Candia, and A. Coniglio, Phys. Rev. E 67, 057105 (2003); M. P. Ciamarra, M Tarzia, A de Candia, A Coniglio,
Phys. Rev. E 68, 066111 (2003).

Weigt2003 M. Weigt and A. K. Hartmann, Europhys. Lett. 62, 533 (2003).

Rivoire2004 O. Rivoire, G. Biroli, O. C. Martin, and M. Mézard,  Eur. Phys. J. B 37, 55 (2004).

Tarzia2007 M. Tarzia,
Journal of Statistical Mechanics: Theory and Experiment, P01010 (2007).

krzakala2008 F. Krzakala and L. Zdeborová, Europhys. Lett. 81, 57005 (2008). L. Zdeborová and F. Krzakala, Phys. Rev. E 76, 031131 (2007).

Kurchan2012 J.  Kurchan,  G.  Parisi,  and  F.  Zamponi, J.  Stat.  Mech. (2012) P10012;  J.  Kurchan,  G.  Parisi,  P.  Urbani,  and  F.  Zamponi, J.   Phys.   Chem.  B 117, 12979 (2013);   P.   Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, J. Stat. Mech. (2014) P10009; Annu.  Rev.  Condens.  Matter  Phys. 8, 265 (2017); T. Maimbourg,  J. Kurchan,  and F. Zamponi, Phys. Rev. Lett. 116, 015902 (2016).

Ritort2003
F. Ritort and P. Sollich, Adv. Phys. 52, 219 (2003).

Chandler2010 D. Chandler and J. P. Garrahan, Annu. Rev. Phys. Chem. 61, 191 (2010).

Garrahan2000
J. P. Garrahan and M. E. J. Newman,
Phys. Rev. E 62, 7670	(2000). 

plaquette J.  P.  Garrahan, Phys.  Rev.  E 89,  030301  (2014); R. M. Turner,  R. L. Jack, and J. P. Garrahan, Phys.  Rev.  E 92,  022115  (2015); R.  L.  Jack  and  J.  P.  Garrahan, Phys.  Rev.  Lett. 116,  055702  (2016).

Sethna1992 J. D. Shore, M. Holzer, and J. P. Sethna, Phys. Rev. B 46, 11376 (1992).

Cannas2
S. A. Cannas, M. F. Michelon, D. A. Stariolo
and F. A. Tamarit,

Phys. Rev. E 78, 051602 (2008).

Cannas
O. Osenda, F. A. Tamarit, and S. A. Cannas, Phys. Rev. E 80, 021114 (2009).

Kivelson1995
D. Kivelson, S. A. Kivelson, X. L. Zhao, Z. Nussinov, and G. Tarjus,  Physica A 
219, 27 (1995).

TarjusAvoided 

I. Esterlis, S. A. Kivelson, G. Tarjus, 
Phys. Rev. B 96, 144305 (2017).

cepas2012
O. Cépas and B. Canals, Phys. Rev. B 86, 024434 (2012).

cepas2014
O. Cépas, Phys. Rev. B 90, 064404 (2014).

Dobro2015 S. Mahmoudian, L. Rademaker, A. Ralko, S. Fratini, and
V. Dobrosavljevic, Phys. Rev. Lett. 115, 025701 (2015).

Budrikis2012
 Z. Budrikis, K. L. Livesey, J. P. Morgan, J. Akerman, A. Stein, S. Langridge, C. H. Marrows, and R. L. Stamps,
 New J. Phys. 14, 035014 (2012).
 
Levis2012 D. Levis and L. F. Cugliandolo, Europhys. Lett. 97, 30 002 (2012). D. Levis and L. F. Cugliandolo, 

Phys. Rev. B 87, 214302 (2013). 

Harris1997 M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Phys. Rev. Lett. 79, 2554 (1997).

VBG M. Tarzia and G. Biroli, Europhys. Lett. 82, 67008 (2008).

Hamp2018  J. Hamp, R. Moessner, and C. Castelnovo, Phys. Rev. B 98, 144439 (2018).

Chioar2016 I. A. Chioar, N. Rougemaille, and B. Canals, Phys. Rev. B
93, 214410 (2016).

Mezard1987 M. Mézard, G. Parisi, and M. A. Virasoro, Spin Glass Theory and Beyond (World Scientific, Singapore, 1987).

Kirkpatrick1989 T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys. Rev.
A 40, 1045 (1989).

Lubchenko2007 V. Lubchenko and P. G. Wolynes, Annu. Rev. Phys. Chem. 58,
235 (2007).

Mezard2001 M. Mézard and G. Parisi, Eur. Phys. J. B 20, 217 (2001).

Ogielski1985

A. T. Ogielski, 
Phys. Rev. B 32, 7384 (1985)

OgielskiMorgenstern1985

A. T. Ogielski and I. Morgenstern, 
Phys. Rev. Lett. 54, 928 (1985).

Franz1994 S. Franz, G. Parisi, M. A. Virasoro, Journal de Physique I 4, 1657 (1994); V. Astuti, S. Franz, and G. Parisi, 
Journal of Physics A: Mathematical and Theoretical 52, 294001 (2019).

Mattson1992
J. Mattsson, P. Granberg, P. Nordblad, L. Lundgren, R.Loloee, R.Stubi, J. Bass, J. A. Cowen, 
J. Magn. Magn. Mat. 104–107, 1623 (1992).


Gucchhait2017
S. Guchhait and R. L. Orbach, Phys. Rev. Lett. 118, 157203 (2017).

Fernandez2019

L. A. Fernández, E. Marinari, V. Martín-Mayor, G. Parisi 
and J. J. Ruiz-Lorenzo,  J. Phys. A: Math. Theor. 52,  224002
(2019).

Wills2002 A. S. Wills, R. Ballou,  and C. Lacroix, Phys. Rev. B 66, 144407 (2002).

remark Note that the fact that the mean-field value of the ground state energy is found to be lower than the 2d value is due to the fact that the dipolar interactions are cut-off beyond the  second  nearest-neighbour plaquettes in the mean-field calculations.

Han2008 Y. Han, Y. Shokef, A. M. Alsayed, P. Yunker, T. C. Luben-sky,  and A. G. Yodh, Nature 456, 898 (2008).

Zhou2017 D.   Zhou,   F.   Wang,   B.   Li,   X.   Lou, and   Y.   Han, Phys. Rev. X 7, 021030 (2017).

Chioar2014 I. A. Chioar, N. Rougemaille, A. Grimm, O. Fruchart,
E. Wagner, M. Hehn, D. Lacour, F. Montaigne, and B. Canals, Phys. Rev. B 90, 064411 (2014).

Nisoli2013 C. Nisoli, R. Moessner,  and P. Schiffer, Rev. Mod. Phys. 85, 1473 (2013).

Ni2008 K.-K.  Ni,  S.  Ospelkaus,  M.  H.  G.  de  Miranda,  A.  Pe'er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin,  and J. Ye, Science 322, 231 (2008).

Griesmaier2005 A.  Griesmaier,  J.  Werner,  S.  Hensler,  J.  Stuhler,    and T. Pfau, Phys. Rev. Lett. 94, 160401 (2005).

Scheie2016 A. Scheie, M. Sanders, J. Krizan, Y. Qiu, R. J. Cava,  and C. Broholm, Phys. Rev. B 93, 180407 (2016).

Paddison2016 J.  Paddison,  H.  Ong,  J.  Hamp,  P.  Mukherjee,  X.  Bai,
M.  Tucker,  N.  Butch,  C.  Castelnovo,  M.  Mourigal,   and S. Dutton, Nat. Commun. 7, 1 (2016).

Dun2017 Z. L. Dun, J. Trinh, M. Lee, E. S. Choi, K. Li, Y. F. Hu,Y. X. Wang, N. Blanc, A. P. Ramirez,   and H. D. Zhou, Phys. Rev. B 95, 104439 (2017).

Takagi1993 T. Takagi and M. Mekata, J. Phys. Soc. Japan 62, 3943 (1993).

Levis2013
D. Levis, L. F. Cugliandolo, L. Foini, and M. Tarzia,
 
    Phys. Rev. Lett. 110, 207206 (2013) 

Foini2013
 L. Foini, D. Levis, M. Tarzia and L. F. Cugliandolo,
  
   J. Stat. Mech.  P02026 (2013).

Cugliandolo2015
   L. F. Cugliandolo, G. Gonnella and A. Pelizzola,

 J. Stat. Mech. P06008 (2015). 
 
Pelizzola2005 

E. N. M Cirillo, G. Gonnella, D. A. Johnston, and A. Pelizzola, 
Phys. Lett. A 226, 59 (1997);
A. Pelizzola, J. Phys. A: Math. Gen. 38,
R309 (2005).

RRG The properties of random-regular graphs have been extensively studied. 
For a review see N. C. Wormald, Models of random-regular graphs, in Surveys in Combinatorics, J. D. Lamb and D. A. Preece, eds., London Mathematical Society Lecture Note Series 276, 239 (1999).

Castelnovo2008 C.  Castelnovo,  R.  Moessner,   and  S.  L.  Sondhi,  Nature 451, 42 (2008).

cutoff P. Andriushchenko, 
Journal of Magnetism and Magnetic Materials 476, 284 (2019).

remark2 In particular, it can be shown that m= 1 always gives back the paramagnetic RS solution.

Tierno1

A. Le Cunuder, I. Frerot, A. Ortiz-Ambriz, A. Ortíz-Ambriz,
and P. Tierno, Phys. Rev. B 99, 140405 (2019).

Tierno2

J. Loehr, A. Ortiz-Ambriz, A. Ortíz-Ambriz, and P. Tierno, 
Phys. Rev. Lett. 117, 168001 (2016). 

Tierno3

A. Ortíz-Ambriz and P. Tierno, 
Nat. Comm. 7, 10575 (2016).

Colloidal

A. Libal, C. J. O. Reichhardt-Olson and C. Reichhardt, 
New J. Phys. 17, 103010 (2015).

Heyderman1
A. Farhan, A. Kleibert, P. M. Derlet, L. Anghinolfi, A.  Balan, R. V. Chopdekar, M. Wyss, S.  Gliga, F. Nolting,  L. J. Heyderman,

Phys. Rev. Lett. 111, 057204 (2013);  

Phys. Rev. B 89, 214405 (2014).

Luning

O. Sendetskyi, V. Scagnoli, N. Leo, L.  Anghinolfi, A. Alberca, J.  Luning, U. Staub, P. M. Derlet, L. J.  Heyderman, 
Phys. Rev. B 99, 214430 (2019). 




