Exact numerical solution of the fully-connected classical and quantum Heisenberg spin glass

,

While the physics of classical and quantum Ising spin glasses has been rather thoroughly understood, glasses of Heisenberg (vector) spins have remained a difficult and largely unsolved problem, including especially its quantum version, which governs the local moments in randomly doped, strongly correlated materials.
The approach of Sachdev and Ye [1] (SY), which takes a double limit of a fully connected exchange model with SU(2) spins promoted to SU(M ≫ 1), has attracted a lot of attention, as it exhibits very interesting features in the M = ∞ limit.In a fermionic representation of the SU(M → ∞) generators, the spins do not freeze, but remain in a spin liquid state [1] similarly to the related SYK model of randomly coupled Majorana fermions [2].When expanding around M = ∞, it was conjectured that a transition into a glassy phase occurs at an exponentially small critical temperature log(T g ) ∝ − √ M [3], and a recent study has established that this phase displays full replica symmetry breaking (RSB) [4].A glass phase with similar features was found to occur in the transverse field Ising model [5,6] and the related multicomponent quantum rotor models [7][8][9].In contrast, for bosonic SU(M ) representations, a spin glass phase occurs even for M = ∞, albeit with one-step RSB [3,10].
There is a fundamental difference between glasses displaying full RSB and one-step RSB.The energy landscape of the former features marginally stable local minima with a gapless spectrum of excitations, while the latter display a fully stable, gapped ground state, lying far below the manifold of marginally stable excited states that trap the dynamics.Marginal stability is key to the dynamics and the physics of avalanches in full RSB systems [11].Which of these two scenarios applies to the physically relevant Heisenberg SU(2), S = 1/2 spin glass is an open question.Exact diagonalization of small fully-connected systems is limited by finite size effects [12,13].Surprisingly, not much is known about the classical (large S) limit of the Heisenberg mean-field spin glass either, apart from the replica-symmetric analysis of Refs.[14,15].
The equilibrium quantum spin dynamics of the SY spin-liquid [1] is similar to that of a marginal Fermi liquid [16].It was realized early on [17] that this opens a new perspective on 'strange' metals, culminating in recent models in which disorder and strong interactions conspire to prevent the emergence of quasiparticles and lead to T -linear resistivity down to the lowest temperatures [2,18].It has recently been shown that such a behavior is found in the quantum critical region around the melting point of a metallic SU(2) Heisenberg spin glass [19,20].However, it is not known how the quantum dynamics change within the metallic spin glass phase itself.
In this Letter, we answer these open questions and present a solution of the fully connected spin-1/2 Heisenberg spin glass throughout its ordered phase.In the quantum case, the model reduces to an impurity problem for the spin dynamics coupled to Parisi's equations for the spin glass order.The Heisenberg glass has a full RSB solution, whose structure, however, differs crucially from an Ising glass.The insulating quantum glass displays the universal spin dynamics found in all solvable mean-field quantum glasses with full RSB in the limit T → 0, but deviations from it persist to surprisingly low temperatures.Upon doping, we find a metallic spin glass with unexpectedly slow spin dynamics, close to the SY dynamics that dominate the quantum critical melting point of the spin glass.

arXiv:2312.14598v3 [cond-mat.dis-nn] 25 Aug 2024
Model and method.We consider N spins S i with all-to-all interactions, described by the Hamiltonian The S i are either classical Heisenberg spins -vectors constrained to the three-dimensional sphere of radius S = 1/2 -or quantum SU(2) spins S (with S 2 = S(S + 1)).
We use ℏ = k B = 1 and denote by ℓ = 3 the number of spin components.J ij are Gaussian random couplings with zero mean and variance J 2 /N .We obtain the equilibrium solution of this model using the replica method and Parisi's full RSB ansatz [21][22][23][24][25][26], as detailed in the Supplementary Information [27].
In brief, the mean field model is reduced to an effective single spin problem in a random frozen field h with distribution P(h).In the classical case, it is governed by the Hamiltonian H loc (h) = −h • S. In the quantum case, the dynamics of this spin also depend on the self-averaging spin autocorrelation function χ(τ ) = ⟨S(0)S(τ )⟩ − ⟨S⟩ 2 , via the action where τ is Matsubara imaginary time and β = 1/T the inverse temperature.The glass phase is described by an order parameter q(x) (x ∈ [0, 1]).q(x) characterizes the distribution of phase space distances between local minima of the free energy landscape [25].It determines the local field distribution P(h) ≡ P(x = 1, h), which is found by solving Parisi's equations for the magnetization s(x, h) and frozen fields P(x, h) These equations are solved self-consistently, with and, in the quantum case, Such self-consistency conditions are an example of extended dynamical mean-field theory (EDMFT) [5,[28][29][30].The iterative procedure to solve them is summarized in Fig. 1.A similar procedure has recently been implemented for mean field versions of transverse field Ising and quantum Coulomb glasses [6,32].In contrast to the classical case, the single site quantum problem (2) cannot be solved analytically.Its solution is obtained with a "CTSEG" continuous-time quantum Monte Carlo (CTQMC) algorithm [27,31], without fermionic sign problem.The Parisi equations are solved with a high precision numerical method (error bars on q(x) being of order 10 −5 ), using their integral form [33] and filtering methods to suppress numerical instabilities at low T [27].
Glass phase.A spin glass phase appears below a critical temperature T g , where the Edwards-Anderson order parameter q EA ≡ q(x = 1) turns on, as illustrated on Fig. 2a.Our results for T g are consistent with previous analytical work: T g = JS 2 /3 (≈ 0.08J for S 2 = 1/4) in the classical case [34] and T g ≈ J⟨S 2 ⟩/3 √ 3 ≈ 0.15J in the quantum case [22], with q EA (T ) following the prediction of [34] close to T g (Fig. 2a, inset).Note that T g is higher in the quantum case, since the quantum spins are larger (S(S + 1) > S 2 ).
For T ≪ T g , the overlap function q(x) obeys an approximate scaling form q(x, T ) = q EA (T )f (x/T ) + O(T /J) (Fig. 2a and Fig. S13).The function f is determined by solving the Parisi equations (3ab) directly at T = 0, upon changing to the natural variable u = βx.In the quantum case, they require as boundary conditions the zero temperature magnetizations, ⟨S(τ )⟩ S loc (h) , which we approximate by our lowest temperature QMC results (see [27] Fig. S5-S6).The results for q(x) are shown in Figs.2a,b.The maximal possible overlap is q max = ⟨S 2 ⟩/ℓ, corresponding to a product state without fluctuations.As expected, the classical glass approaches this value at T = 0, while quantum fluctuations weakly reduce the order parameter to q EA ≈ 0.81q max (consistent with exact diagonalization in Ref. [12] , while the fit of spin correlations in Ref. [13] likely underestimated the order).Unlike in the Ising case, local stability does not require a pseudogap in the distribution of local fields.Instead P(h = 0) remains finite and P(h) ≈ P(0) + ah (cf.Fig. 2b), since the probability of dangerously small fields is already suppressed by the phase space factor ∼ h ℓ−1 [35][14].This in turn is related to the tail of the overlap function q(x).Indeed, (4) suggests that for x close to 1, q EA − q(x) counts the number of spins that see frozen fields of order T /x, for which one may expect q EA − q(x) ∼ T /x 0 dh h 2 P(h) ∼ (T /x) 3 + O (T /x) 4 , consistent with an apparent power law: ∼ (T /x) α .Numerically we indeed find an apparent power law approach: 1 − q(x)/q EA ∼ 1/(βx) α as T → 0, with α slightly larger than 3.This contrasts with the Ising case [5,25,35], where the linear pseudogap in P(h) leads to α = 2.The lower density of small fields in the Heisenberg case suggests that low-lying metastable states have higher energies.This is consistent with the smaller value of the so-called breakpoint x c , above which q(x) reaches a plateau ([27], Fig. S12).T /x c (T ) can be interpreted as the typical free energy difference between the lowest metastable states [36].While in the Ising case that energy scales linearly with T as x c (T → 0) ≈ 0.5 [35], in Heisenberg glasses it decreases much more slowly as T /x c (T ) ∼ 1/ log(1/T ), corroborating the picture of a rougher energy landscape.These differences will further show in the response to an increasing external field, under which the ground state magnetization increases in random discontinuous steps called 'shocks'.Their size distri-bution ρ(∆m) was numerically found to be very similar to that of field-triggered out-of-equilibrium avalanches in the classical SK model [37,38], a phenomenon attributed to the marginal stability of the full RSB landscape.The equilibrium ρ(∆m) is governed by the asymptotic approach of q(x) to q EA .If it scales as (T /x) α , then ρ(∆m) ∼ 1/(∆m) 2/α for N −1/2 ≪ ∆m ≪ 1 [38].For the classical Ising glass with α = 2, ρ(∆m) ∝ 1/(∆m) exhibits a broad spectrum of avalanches.The larger α ≈ 3 of Heisenberg glasses leads to ρ(∆m) ∝ 1/(∆m) 2/3 , with predominant large scale rearrangements -as numerically observed in avalanches of XY (ℓ = 2) spins.[39] Specific heat.Let us now consider the internal energy per spin.In the classical case, it is given by [27]: In the quantum case, denoting Q(τ ) = χ(τ ) + q EA and τ = τ /β, it reads As shown in Fig. 2c, the classical and quantum internal energies behave very differently as a function of temperature.For classical vector spins, one expects a constant intra-state heat capacity c = (ℓ − 1)/2 as T → 0, each degree of freedom contributing 1/2 by the equipartition theorem.A linear fit to our data yields c ≈ 1.0 ± 0.1, excluding sizeable inter-state contributions, consistently with finite-size simulations [40].Quantum effects gap out the soft degrees of freedom and yield a much weaker temperature dependence of the internal energy (Fig. 2c).A heat capacity proportional to T 3 was predicted for the SU(N ) quantum spin glass in the large N limit, provided that marginally stable rather than equilibrium states are analyzed [41].Our data are compatible with U (T ) ∝ T 4 , but only at the lowest temperatures T /J ≲ 0.02, as finite temperature corrections are substantial.Spin susceptibility.The spin-spin correlator χ(τ ) is shown in Fig. 3a.For large τ it is well described by the conformal scaling form χ(τ ) ≈ χ(β/2)/ [sin(πτ /β)] θ [17], implying that for ω ≳ T , the dissipative part of the susceptibility at real frequencies χ ′′ (ω) scales as ω θ−1 .The exponent θ has significant, slow T -dependence (Fig. 3b).Using a Landau expansion, Ref. [8] predicted that θ(T = 0) = 2, (χ ′′ (ω) ∝ ω).This value was also found in a 1/M expansion of the SU(M ) quantum spin glass [4] and actually holds for all solvable cases of marginally stable states found so far [42].From the limited numerically accessible temperature range, it is hard to unambiguously conclude whether lim T →0 θ(T ) = 2.However, Fig. 3b suggests that this limit is indeed approached, albeit very slowly, as we can fit our data by θ(T ) ≈ 2 − 5.2 × T /J.Interestingly, a similar behavior was found recently in the transverse-field Ising spin glass [6].
Metallic quantum spin glass.In order to study the interplay between electrons and frozen spins arising from doping the spin glass, we use the Hamiltonian: where c iσ and c † iσ are the electronic annihilation and creation operators, and U is the on-site electron-electron interaction.The hopping amplitudes t ij are randomly distributed with variance t 2 /N .We denote the doping by p = ⟨n ↑ + n ↓ ⟩ − 1.This model has previously been solved in the paramagnetic phase [20].Here, we solve the spin glass phase; as in [20], we use U = 4t and J = 0.5t.
Below a critical doping p c , a metallic spin glass appears, i.e. a phase with both a non-zero spin glass order parameter q EA and a non-zero density of states at zero frequency, as illustrated in Fig. 4.This is compatible with the spin glass instability of the paramagnetic solution found in [20].Previous studies of metallic quantum spin glasses have found a local spin susceptibility χ ′′ (ω) ∝ √ ω corresponding to an exponent θ = 3/2 at T = 0 [28,42,43].This can be interpreted as originating from Ohmic damping of the oscillators in the physical picture discussed above.Intriguingly, we find that θ is smaller than unity for the whole range of (low) temperatures investigated.Our data do not appear consistent with θ reaching 3/2 at T = 0 (possibly due to non-Ohmic damping by the non-Fermi-liquid metal) but may be consistent with θ reaching 1.The latter corresponds to the quantum critical dynamics found at the critical point p = p c , which might extend through a large part of the metallic spin-glass phase.
In conclusion, we have solved the Heisenberg quantum spin glass for SU(2) spins.We found that the energy landscape of vector spins differs significantly from the Ising counterpart, resulting in different long time dynamics.Upon decreasing T , the short time quantum dynamics slowly approach the super-universal behavior of marginal mean field glasses, although featuring significantly softer collective modes and a lower freezing temperature than comparable Ising glasses.The marginal Fermi liquid-type quantum dynamics anticipated from the SU(M ≫ 1) approach is essentially absent in the undoped insulating limit of SU (2) spins, but appears to be present in a wide window influenced by the dopinginduced quantum critical point in a metallic regime, that may be relevant for strongly correlated doped materials.
The data and code associated with the paper are publicly available [44,45].We thank Philipp Dumitrescu, Subir Sachdev and Nils Wentzell for fruitful discussions.The Flatiron Institute is a division of the Simons Foundation.N.K. acknowledges support from a Humboldt fellowship.In this section, we describe the analytical procedure that reduces the mean-field (infinite dimensional) Heisenberg spin glass problem to a self-consistent single spin problem, amenable to an exact numerical solution.We make use of the replica method, as introduced by Sherrington and Kirkaptrick in the classical case [1] and by Bray and Moore in the quantum case [2].We then follow Parisi's seminal replica symmetry breaking construction [3,4], which was generalized to the quantum case of the transverse field Sherrington-Kirkpatrick model [5,6], but not to the quantum Heisenberg spin glass.

Model definition
We consider an all-to-all interacting system of N spins S i , described by the Hamiltonian The J ij are independent random variables with distribution G(J ij |J 2 /N ), where G(x|v) is the normalized Gaussian with zero mean and variance v.We will consider both a classical and a quantum version of the model.In the classical case, the S i are vectors constrained to reside on the ℓ-dimensional sphere of radius S. In the quantum case, the S i are operators, forming a representation of spin S of the group SU(2); we denote ℓ = 3 the number of spin components.

Replica trick
Our aim is to compute f , the thermodynamic limit of the free energy per site at inverse temperature β.To this end, we introduce the quantity f n (N ) ≡ (−1/βN ) log(Z n ), where the bar denotes averaging over the disorder.Z n is the partition function of n replicas of the N -spin system, with the same realization of the disorder.Allowing n to take non-integer values through analytical continuation, we find Since the free energy is self-averaging in the thermodynamic limit, We will assume that the order of the limits can be reversed, and first carry out the computation of f n in the thermodynamic limit.

Classical case
We average over the random couplings the partition function of n replicas of the spin system, labeled by the index a: = Here, the indices α, β denote the spin components.We may drop in the following the term proportional to n 2 , since it vanishes in the n → 0 limit.We now perform a Hubbard-Stratonovitch transformation to decouple the sum over spins: with In the thermodynamic limit, the integral in Eq. ( 7) is dominated by the saddle point {Q αβ⋆ ab }: with ∂F/∂Q αβ ab | Q αβ⋆ ab = 0, ∀(a, b, α, β).The latter amounts to

Quantum case
In the quantum case, we use an imaginary time path integral formalism to express the partition function: where the notation S [.] is now taken to mean the path integral [DS(τ )]e −S 0 [S(τ )] [.], where S 0 [S(τ )] is the imaginary time action describing an isolated spin, which we do not need to specify explicitly.We may then proceed as in the classical case: = We may again drop the term proportional to n 2 as it vanishes in the n → 0 limit.The Hubbard-Stratonovitch transformation yields with The saddle point approximation in the thermodynamic limit yields with ∂F/∂Q αβ ab (τ, τ ′ )| Q αβ⋆ ab (τ,τ ′ ) = 0, ∀(a, b, α, β, τ, τ ′ ).The latter amounts to

Interpretation of the Q ab
We have reduced the initial lattice problem to a self-consistent local problem of n coupled replicas in the limit n → 0, which remains to be defined.We now establish the link between the spin-spin correlation functions of the local problem and those of the initial lattice problem.This will allow us to simplify some of the coefficients Q ab that define the local problem.For generality, we present here the reasoning in the quantum case.
Let S be the action corresponding to the initial lattice problem.We supplement it with a source term: Then, in the thermodynamic limit, Since the correlation function on the left hand side vanishes for two different spin components, the Consider now two copies of the initial system, governed by the action By applying the replica method to this action in the same way as above, we find the correspondence As the copies 1 and 2 are independent, the right hand side of the above equation, and therefore the Q αβ⋆ ab , a ̸ = b, cannot depend on τ, τ ′ .Similarly, in the absence of quadrupolar order in the initial system, the Q αβ⋆ ab vanish when a ̸ = b and α ̸ = β.The reasoning is analogous in the classical case.Thanks to the normalization of the classical spins, Eq. ( 8) may be further simplified by integrating over the Q αα aa in Eq. ( 7).Let us summarize our results so far.
In the classical case, the free energy per spin is with In the quantum case, with and We now need to specify an ansatz for the matrix Q ab which allows for analytic continuation to non-integer matrix dimensions.

Recursive construction
We use Parisi's replica symmetry breaking (RSB) ansatz for the matrix Q ab .The matrix Q RSB k at the k th stage of RSB is defined by the choice of two sequences: (q p ) p=0..k and (m p ) p=0..k , the latter being strictly decreasing, with m 0 = n.It is built "from the inside out" according to the following prescription: Here U m is the matrix of size m filled with ones of size, diag r (A) is the matrix with r times the block A on the diagonal, and we set q −1 = 0. We note that this construction places q k instead of 0 on the diagonal, and the corresponding term needs to be subtracted in the end.In the classical case, the local Hamiltonian at the k th stage of RSB is therefore which we may write as In the quantum case, we introduce the notation S ≡ 1 0 dτ S(τ ).Then, the local action at the k th stage of RSB is (30) For simplicity, we denote We now establish a recursion between the local partition functions Z loc at steps p and p + 1 of the recursive construction.We use the quantum notation, but the relation will be identical in the classical case.At step p the local partition function under an external field h reads: We note that h is introduced only as an auxiliary term, and we will set h = 0 in the end.In particular it does not represent an actual external field applied on the initial system: the matrix Q ab would then no longer be rotationally invariant.At the first step (p = k), we find At step p, we may consider separately the effect of the two terms in Here, α labels the blocks of the matrix Q RSB k p .After a Hubbard-Stratonovitch transformation on the squared sum, we find the recurrence relation We now take the simultaneous limit n → 0 and k → ∞, keeping m k = 1.Then, m p becomes a continuous variable x between 0 and 1, and q p → q(x), such that m p → x. q(x) is the Parisi order parameter.The corresponding limit of the recurrence relation (35) may be taken if one uses the formal identity that holds for any sufficiently regular function f .Z RSB k p (h) then becomes a continuous function This may be recast into a partial differential equation: with boundary condition ζ(1, h) = z ∞ (h) as defined below.At this point, a distinction needs to be made between the classical and the quantum case.In the classical case, the "single-replica" part of the local Hamiltonian H k does not depend on the spin.We may hence split it off in the partition function, so that applying Eq. ( 35) at p = 0 yields and the boundary condition is explicitly We may define ϕ(x, h) ≡ (1/x) log ζ(x, h), which satisfies We also anticipate that q(0) = 0, so that the Gaussian in Eq. ( 39) becomes a δ function.
Then, introducing our results into Eq.( 23) the free energy of the classical Heisenberg spin glass is given by where ϕ(x, h) satisfies Eq. ( 41) with ϕ(1, h) = log z Cl ∞ (h).
We have used that, with the RSB ansatz [7], This provides the solution of the classical spin glass problem given a function q(x).
In the quantum case, the single-replica part of the action cannot be simplified and we obtain lim and with Then, introducing our results into Eq.( 25) the free energy of the quantum Heisenberg spin glass is given by where ϕ(x, h) satisfies Eq. ( 41) with This provides the solution of the quantum spin glass problem given q(x) and Q(τ ).

Local observables and two point functions
We now determine the functions q(x) and Q(τ ) by applying the RSB ansatz to the self-consistency conditions in Eqs. ( 24) and ( 26).We go back to the recursive construction in 1.3.1 and consider the observable R a (τ, τ ′ ) = S a (τ )S a (τ ′ ), which is local in replica space.For every τ, τ ′ , we define the quantity R RSB k p (h), which corresponds to the average value of R computed with the matrix : indeed, at a given step, only one of the blocks is perturbed by the insertion of R, while the others are not.Therefore, one has the recurrence relation The situation is slightly different for an observable that is non-local in replica space, in particular the spin-spin correlator ⟨S . (49) The corresponding PDE is with the boundary condition ρ(1, h) τ,τ ′ = ⟨S(τ )S(τ ′ )⟩ ∞,h , where ⟨•⟩ ∞,h is the mean value with respect to H ∞ (h) or S ∞ (h).If we now define βs(x, h) = (1/x)∇ log ζ, we find that it satisfies the same equation as ρ, with the boundary condition s(1, h) = ⟨S⟩ ∞ .Then ρ satisfies Using Eq. ( 48) with p = 0, The spin-spin correlation function ⟨S a S b ⟩ RSB k p becomes a continuous function ν(y, x, h), where y is the continuum analogue of the index p 0 defined above.For y < x, ν satisfies Eq. ( 52), with the boundary condition ν(y = x, x, h) = s 2 (x, h).The self-consistency conditions ( 24) and ( 26) become This completes in principle our solution of the spin glass problem.One only needs to compute the average magnetization and the time correlation function in a single spin problem governed by H ∞ or S ∞ .This yields, through integration of Eq. ( 52), the functions q(x) and Q(τ ), which allow one in turn to compute the free energy, starting from the partition function associated with H ∞ or S ∞ .However, these results are formulated much more conveniently by introducing the probability distribution of the local field.

Probability distribution of the local field
We introduce the probability distribution P(h) of the local field acting on the effective impurity problem.Let us define the function P (x, h|x ′ , h ′ ) as the solution of the "backward" eq. ( 52) with boundary condition For any x ′ ≥ x, it has the Markovian property Differentiating with respect to and since P satisfies (52), we obtain Integrating by parts, Since this is true for any h ′′ , P (x, h|x ′ , h ′ ) is found to satisfy the "forward" equation We now rewrite eq. ( 53) as having defined P then satisfies eq. ( 59) with boundary condition P(0, h) = G(h|2J 2 q(0)) = δ(h).Similarly, for the Parisi order parameter we find

Final formulation of the solution
We are now in position to formulate the solution in its final form.
We have reduced the solution of the fully connected Heisenberg spin glass problem to the solution of a single spin problem.It is governed in the classical case by the Hamiltonian H ∞ (h) = −h • S and in the quantum case by the action The Parisi order parameter and local susceptibility are obtained as where ⟨•⟩ ∞,h denotes the mean value with respect to H ∞ (h) or S ∞ (h).P(x, h) and s(x, h) are obtained by solving Once q(x) and Q(τ ) are known, the free energy can be evaluated according to Eqs. ( 42) and (47).
In practice, the local problem is solved trivially in the classical case: and the self-consistency needs only to be established between Eqs. (64) through (66).In the quantum case, the local problem is still a many-body problem, which can only be solved numerically.Moreover, the local action depends self-consistently on the observables q(x) and Q(τ ), so that convergence must be reached for two interconnected self-consistency loops (see Fig. 1 of the main text).

Internal energy
The internal energy per spin U/N can in principle be obtained by numerically differentiating the free energy with respect to temperature.However, it can be obtained in a more straightforward way by carrying out the differentiation analytically before replica symmetry breaking.Going back to Eqs. ( 9) and ( 16), we compute the internal energy as Since the Q ab were introduced as integration variables, they may be arbitrarily rescaled.If we define Qab = β 2 J 2 Q ab , the "complicated" terms in Eq. ( 8) and ( 15) no longer have an explicit β dependence.We then obtain, in the classical case, We then use the property lim n→0 a<b Similarly, in the quantum case, so that

Doped quantum spin glass
Our solution is readily generalized to a "doped" Heisenberg spin glass model.Physically, the randomly-interacting spins are now being carried by electrons, which can hop between the sites of a fully connected lattice and experience an on-site repulsion U .The corresponding Hamiltonian is with Here the c iσ and c † iσ are the electronic annihilation and creation operators, σ =↑, ↓ labels the spin state, n iσ = c † iσ c iσ is the occupation number on site i, U is the on-site electron-electron interaction, µ is the chemical potential and the t ij are hopping amplitudes, randomly distributed with variance t 2 /N .The spins are defined in terms of the fermionic operators as S α i = (1/2) σ,σ ′ c † iσ σ α σσ ′ c iσ ′ , with σ α the Pauli matrices.The solution of this model proceeds along the same steps as above, and here we only outline the main differences.We represent the partition function as a path integral over the Grassman variables (c † iσ (τ ), c iσ (τ )): Note that we now let the imaginary time run between 0 and β.We then compute Z n , with the disorder average now being carried out over both the couplings J ij and the hoppings t ij .The Hubbard-Stratonovitch transformation now introduces integration variables Q ab and ∆ ab , that decouple respectively the spin-spin interaction term and the hopping term.Upon performing the relevant simplifications (for instance, the ∆ a̸ =b vanish), and dropping indices for simplicity, with The saddle point equations read Upon replica symmetry breaking, the original lattice problem is reduced to a self-consistent impurity problem, governed by the action The self-consistency conditions read q(x) = 1 P(x, h) and s(x, h) are obtained by solving Eqs. ( 65) and (66).
We note that we have used in Eq. ( 83) the spin-symmetrized Green's function G(τ ) = (G ↑ (τ ) + G ↓ (τ ))/2.In the absence of an external magnetic field, it is equal to the Green's function of either spin.In the presence of a field h, G ↑ and G ↓ depend on both the magnitude and direction of the h, but G depends only on ||h||.Indeed, let us apply a unitary transformation U to the basis of spin states {|σ⟩}.This defines new basis vectors {|σ ′ ⟩}.We may define creation and annihilation operators for spin up and spin down electrons in this new basis.They are given by The spin up Green's function in the new basis is Similarly, Now, the unitarity of U imposes that The spin-symmetrized Green's function is thus rotationally-invariant even under an external magnetic field.
2 Numerical procedures

Integral equation formulation
The Parisi equations Eqs. ( 65) and (66) are non-linear partial differential equations that need to be solved numerically.We introduce the Green's function G(x, h|x ′ , h ′ ) of their linear part, which satisfies G is simply a Gaussian: If we now treat the non-linear terms as source terms, Eqs. ( 65) and (66) can be formally integrated according to We now use the fact that s(x, h) is always in the direction of h: s(x, h) = s(x, h) ĥ.Then, Hence, we may introduce angle-integrated Green's functions: Eqs. ( 92) and ( 93) then reduce to one-dimensional integral equations: and

Numerical procedure
Our solver code for Eqs. ( 98) and (99), written in Python, is available on Zenodo [12].The variable x is discretized, and Eqs. ( 98) and (99) are solved on successive x points.We typically use 128 to 256 x points in the region where q(x) is non-constant.When q(x) and q(x ′ ) are very close, the Green's functions G(x, h|x ′ , h ′ ) are strongly peaked around h = h ′ .However, the result of the convolution varies much more slowly as a function of h.Therefore, we sample h ′ on a very fine grid (typically 20000 points) and h on a much coarser Chebyshev grid (typically 256 points), and then use barycentric interpolation to obtain the result of the convolution on the fine grid for the next integration step.When the non-linear term is large (typically, when q(x) is large), the numerical integration procedure is unstable.Indeed, at step n, the non-linear term contains a contribution from the n th numerical derivative of the initial condition, which develops an oscillatory instability with a frequency on the order of the Chebyshev grid spacing.As long as the grid is much finer than the structure in P(x, h) and s(x, h), the spurious oscillations may be filtered out.We use Savitzky-Golay filtering, which amounts to splitting the h interval into panels, and fitting P(x, h) or s(x, h) with a second degree polynomial in h on each of the panels.The panel size is determined adaptively, as a fraction of the value of h where the non-linear term has a maximum.With this procedure, and given our level of discretization, the filtered-out part is ≲ 10 −5 of the remainder at every x step, which roughly sets the precision level of our solver.We further check that the normalization ∞ 0 4πh 2 P(x, h)dh = 1 and the sum rule [8] 1 in the classical case are satisfied to a precision better than 10 −4 after self-consistency has been reached.

Solution at T = 0
Upon the change of variable u = βx, the Parisi equations ( 65) and (66) become In this form, they may be solved directly in the limit β → ∞ (T = 0).We then impose the boundary condition for s at a finite u max , ensuring that the value of u max does not affect the solution in the relevant range of u.To obtain the data in Fig. 2b, we used u max = 100/J.

Solution of the quantum impurity problem
We solve the quantum impurity problem defined by the action in Eq. ( 81), and its restriction in Eq. (63), by using a continuous-time Monte Carlo procedure based on a hybridization expansion in the segment picture (CTSEG algorithm) [9,10].Our quantum Monte Carlo code, based on the TRIQS library, is available on GitHub [11].Briefly, the algorithm is based on the double expansion of the hybridization (∆) and spin-spin interaction (Q) terms.A configuration is defined by the order of the expansion and by the choice of points in imaginary time.Several possible updates or "moves" are defined for the configuration, which are accepted or rejected according to the Metropolis prescription.Further details of our implementation will be described elsewhere.

Classical case
In the classical case, the magnetization ⟨S⟩ ∞,h is known analytically and does not depend on the spin glass parameters.Therefore, we start from an initial guess for q(x) and iteratively solve Eqs.(66) and (65) until convergence.As self-consistency is approached, q(x) develops a singularity: q(x) = q(1) for x > x c .Yet, our numerical solution cannot produce a singularity in the mathematical sense.Two procedures can be envisioned for finding the breakpoint x c .
1. We may apply a threshold to dq/dx.Concretely, once a new q(x) has been determined from Eq. ( 64), we enforce that if q(1) − q(x) < ϵ, q(x) = q(1), with typically ϵ = 10 −5 , on the order of our solver's accuracy.This will result in the breakpoint moving as the solution is iterated, and eventually converging.
2. We may iterate the equations until convergence for a range of fixed values of x c , and then choose the x c for which the sum rule in Eq. ( 100) is most accurately satisfied.
We found that procedure 1 converges in a reasonable number of iterations at low temperatures T /J ≤ 0.04, while procedure 2 is required for higher temperatures.At zero temperature, since x c /T logarithmically diverges as T → 0, there is no breakpoint in the finite interval of x/T in which we solve the equations.

Quantum case
In the quantum case, the magnetization ⟨S⟩ ∞,h is obtained by the solving the action in Eq. ( 81) or (63) with quantum Monte Carlo for a range of values of h.Furthermore, the action depends on the spin glass parameters and must be updated until self-consistency is reached.We start from an initial guess for q(x) and Q(τ ) (and ∆(τ ) in the doped case), and determine ⟨S⟩ ∞,h (as well as the spin-spin correlation function, and possibly the electron Green's function) on the impurity, for typically 30 values of h sampled on a Chebyshev grid.The observables at arbitrary values of h are then inferred by barycentric interpolation.Then, given ⟨S⟩ ∞,h , Eqs. (66) and (65) are iterated N times.An arbitrarily large value of N is not necessarily best for rapid convergence of the whole self-consistent procedure; we found empirically that N = 5 is adequate.Procedure 1 (defined above) was sufficient to obtain reasonable convergence of the breakpoint.Once q(x) and P(h) are obtained, we determine the Q(τ ) and ∆(τ ) that define a new impurity action, and repeat the procedure until convergence of q(x), Q(τ ) and ∆(τ ).The Parisi equations (66) and (65) can be solved directly at zero temperature, but not the quantum impurity problem.However, the input for the Parisi equations is only the magnetization, which we find to be almost converged with temperature for βJ > 100 (Fig. S5).We use the magnetization obtained at the lowest accessible temperature βJ = 200 as an input for the zerotemperature Parisi equations.

1 FIG. 4 .
FIG. 4. Doped case.Glass order parameter qEA, density of states at the Fermi level, and scaling exponent θ of the spin-spin correlation function, as a function of doping p at βJ = 50.Dashed lines are guides to the eye.Upper inset: sketch of phase diagram.The spin glass (SG) melts at the quantum critical point (QCP).Lower inset: θ at fixed doping p = 0.1, as a function of temperature.

2 .
RSB k .When k is sufficiently large, a and b are always in different blocks at k-RSB.If at step p they are still in different blocks, then ⟨S a S b ⟩ RSB k p = (⟨S⟩ RSB k p ) Let p 0 by the first step where S a and S b are in the same block.Then for p < p 0 , the spin-spin correlator satisfies the same recurrence relation as a local observable:

Figure S1 :
Figure S1: Parisi order parameter q(x) for the classical Heisenberg spin glass, for a range of temperature values.

Figure S2 :Figure S5 :
Figure S2: Local field probability distribution P(h) for the classical Heisenberg spin glass, for a range of temperature values.

Figure S6 :
Figure S6: Magnetization ⟨S⟩ ∞,h versus temperature, for three magnetic field values.At the lowest accessible temperatures, the magnetization is almost converged to its zero temperature shape.

Figure S11 :
Figure S11: Scaling exponent θ, as a function of temperature in the metallic spin glass phase, for a range of values of the doping p.

Figure S12 :
Figure S12: Breakpoint x c (T ) vs T .xc (T )/T is a characteristic inverse energy scale governing the distribution of free energies of excited states.It scales logarithmically with T in both the classical and the quantum case.Dashed lines are guides to the eye.Inset: definition of the breakpoint x c (q(x > x c ) = q(1)).

Figure S13 :
FigureS13: Difference between q(x)/q EA (T ) and the scaling function f (u = βx), versus temperature.Finite-temperature corrections to the scaling are of order T .

Material for: Exact numerical solution of the fully connected classical and quantum Heisenberg spin glass
Nikita Kavokine, Markus Müller, Antoine Georges and Olivier Parcollet -Contents