Equilibrium dynamics of infinite-range quantum spin glasses in a field

We determine the low-energy spectrum and Parisi replica symmetry breaking function for the spin glass phase of the quantum Ising model with infinite-range random exchange interactions and transverse and longitudinal ($h$) fields. We show that, for all $h$, the spin glass state has full replica symmetry breaking, and the local spin spectrum is gapless with a spectral density which vanishes linearly with frequency. These results are obtained using an action functional$\unicode{x2014}$argued to yield exact results at low frequencies$\unicode{x2014}$that expands in powers of a spin glass order parameter, which is is bilocal in time, and a matrix in replica space. We also present the exact solution of the infinite-range spherical quantum $p$-rotor model at nonzero $h$: here, the spin glass state has one-step replica symmetry breaking, and gaplessness only appears after imposition of an additional marginal stability condition. Possible connections to experiments on random arrays of trapped Rydberg atoms are noted.

Modern advances in the development and control of programmable quantum simulators [1][2][3][4] have led to remarkable implementations of the idea of solving classical optimization problems by quantum tunnelling [5,6].In a recent experiment, for instance, Ebadi et al. [1] used a two-dimensional Rydberg atom array to investigate quantum optimization algorithms and demonstrated a superlinear quantum speedup in finding exact solutions.
In such a setup, each atom can be either in the atomic ground state or in a highly excited Rydberg state (with a large principal quantum number), thus realizing a quantum two-level system, which can be represented by the eigenstates of the Pauli operator Z i on atom i.The laser-induced Rabi flipping is then described by the operator g i X i , while the laser detuning is given by ∆ i Z i .The long-ranged van der Waals interactions between two atoms is active only when both are in the Rydberg state, so an interaction J ij between atoms i and j, J ij (1 + Z i )(1 + Z j )/4, provides a route to implementing pairwise constraints on the optimization problems [1,[7][8][9][10].The atoms, which are trapped in optical tweezers, can be arranged in arbitrary geometries, and positioning them on randomly site-diluted lattices-as in Ref. 1-introduces an element of spatial disorder in the Hamiltonian.
Inspired by the success of random, infinite-range models in understanding classical optimization problems [11], we will examine here the equilibrium dynamics of the infinite-range Ising spin glass in a field, h, with the Hamiltonian where i, j = 1 . . .N denote the lattice sites, and the interactions between them, J ij , are taken as independent random numbers drawn from the probability distribution This model has been much studied in the quantum spin glass literature [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26], but only a few results have been obtained for the model with a nonzero field, h ̸ = 0 [18,25,26], which is an essential ingredient in the optimization toolbox of Rydberg quantum simulators.Here, we shall provide exact results for the equilibrium long-time dynamics of the N = ∞ model with h nonzero.Our results here are a prelude to the study of the experimental case where the couplings are time-dependent, but nonequilibrium results will not be presented here.
The equilibrium solution for the N = ∞ quantum Ising model with independent random interactions requires the self-consistent solution of a (0+1)-dimensional replicated Ising model with long-range interactions [25,26].Such a model is not exactly solvable, so a full solution at all times requires Monte Carlo simulations.However, it was argued in Ref. 18 that a Landau-theory-like strategy of expanding the quantum action functional in powers of an appropriately subtracted Z autocorrelation function (which is also a matrix in replica space) can provide the exact form of the long-time correlations in the N = ∞ model.
Here, we shall implement this strategy for the h ̸ = 0 model, and obtain the low-frequency dynamic spin spectrum, along with the Parisi spin glass order parameter with full replica symmetry breaking.
Given the difficulty in obtaining the exact solution of the N = ∞ quantum Ising model at all times, we will also study a cousin of the quantum Ising model which has been the focus of some attention in the literature, but only at zero field: this is the 'spherical quantum p-rotor model' [20,[27][28][29][30][31].In the literature, this model has been referred to as a 'p-spin' model rather than a 'p-rotor' model.While the distinction between spins and rotors is not important for classical systems, it is crucial for quantum systems.'Spin' usually refers to a quantum degree of freedom whose components do not commute with each other, while the components of a rotor all commute.Accordingly, we will use the 'p-rotor' terminology in this paper.
As an aside, we also note studies of the Heisenberg spin glass model [32][33][34][35][36][37][38][39][40][41][42], which will not be considered in the present paper.In the Heisenberg model, the states on each site i have a twofold degeneracy, unlike both the Ising and rotor models considered here.An important consequence is that there is no trivial paramagnetic state in the Heisenberg model.In contrast, the Ising and rotor models have a gapped paramagnet with a nondegenerate ground state at large g.
We define the spherical quantum p-rotor model here by an imaginary time (τ ) path integral over continuous real rotor/spin coordinates σ i , which are the analogs of the discrete Z i = ±1 Ising spins.The partition function of this model at an inverse temperature β = 1/T , and in a field h is with repeated indices summed over, and σi (τ ) ≡ dσ i /dτ .Importantly, in order to keep the energy finite, we equip the rotors with a spherical constraint i.e., the rotor degrees of freedom lie on an N -dimensional sphere of radius √ N .These rotors interact with p-rotor couplings J i 1 ...ip , and all of these couplings are taken to be independent random variables with distribution where the width of the distribution over couplings is set by the scale J.The factor of N p−1 ensures that the Hamiltonian is of order N , and accordingly, enforces an extensive scaling of the energy and free energy.
Owing to the global nature of the constraint in (1.4), the solution of the spherical p-rotor model is simpler than that of the Ising model with a local constraint Z i = ±1 on every site i.On the other hand, the global constraint also makes the p-rotor model a less physical generalization of the Rydberg experiments.
The solution of the p-rotor model requires analysis of a closed set of self-consistent Schwinger-Dyson equations, in which the self energies are written as polynomials of the Green's functions.Such Schwinger-Dyson equations have been numerically studied earlier at h = 0 and in imaginary time and frequency [27,29]; here, we will extend the numerical solution to h ̸ = 0 and obtain the full dynamic rotor spectrum by direct solution in real-frequency space.We will also compute the Parisi spin glass order parameter and find that as in the spherical classical p-rotor model [43][44][45][46][47] (obtained by taking the g = 0 limit of (1.3)), there is only one-step replica symmetry breaking.

A. Main results
Our results are expressed in terms of a field Q ab (τ 1 , τ 2 ) which is bilocal in imaginary time and a matrix in replica space with indices a, b = 1, . . ., n.This field is related to the spin/rotor autocorrelation functions via where a, b are replica indices.Time-translational symmetry requires that the N = ∞ solution take the form where ν n is a bosonic Matsubara frequency.We choose the following ansatz for Q ab : The replica off-diagonal terms in Q ab are chosen to be time-independent, because there is no correlation between the time evolution of the spins in distinct replicas [18].We have parameterized these off-diagonal terms in terms of a Parisi matrix q ab ; as in the classical spin glass theory [48], this matrix has an ultrametric structure which is characterized by the Parisi function q(u), with q(1) ≡ q EA , the Edwards-Anderson order parameter at T = 0. We have included an additive factor of βq EA in the replica diagonal term of (1.8) for convenience.We will find that this ensures the solution for Q r (τ ) vanishes as τ → ∞ at T = 0. Also, the diagonal components q aa do not appear in the above, and we use this freedom to choose q aa = 0.
For the Ising model studied in Sec.II, we find that low-frequency (ω) spectrum in the replica-symmetrybreaking phase with h ̸ = 0 is the same as that obtained earlier [18,22,26] at h = 0: Note that a gapless spectrum for h ̸ = 0 is more surprising than at h = 0 because there is no Z i → −Z i symmetry that can be broken in the spin glass phase.The replica symmetry breaking is characterized by < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 N 5 K 5 4 t 0 R 7 6 B q I a Y w u < l a t e x i t s h a 1 _ b a s e 6 4 = " x X G h e 7 3 j R K t s 5 3 g 8 X q 9 9 L R v 7 x < l a t e x i t s h a 1 _ b a s e 6 4 = " z M m D t Y / I j r P 2 9 O r C M 5 b t p f 6 U i u 0  < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 N 5 K 5 4 t 0 R 7 6 B q I a Y w u < l a t e x i t s h a 1 _ b a s e 6 4 = " x X G h e 7 3 j R K t s 5 3 g 8 X q 9 9 L R v 7 The Parisi spin glass order parameter for the spherical quantum p-rotor model for p ≥ 3.
the Parisi function q(u) shown in Fig. 1.This function has the same form as that found for the classical Ising model with g = 0 [49].The function q(u) is nonanalytic at u = x, and the value of x vanishes linearly with T as T → 0 as (see Eq. (2.27)) where y and κ are couplings of order unity in the Landau action (see Eq. (2.8)).Thus, although the replica symmetry breaking vanishes at T = 0, it is nevertheless important that the T -and h-dependent structure in Fig. 1 be included to obtain the gapless spectrum at T = 0 in (1.9).The smaller u plateau in Fig. 1 is at [18,49] The results for the spherical quantum p-rotor model, p ≥ 3, studied in Sec.III, do have differences from the Ising model tied to the one-step replica symmetry breaking in the former, shown in Fig. 2, which is in turn tied to the nonlocal rotor constraint in (1.4).This q(u) is characterized by a break-point at u = x, and the large-N saddle-point equations leave the value of x undetermined.The mathematical solution of the classical problem at g = 0 [47] implies we should compute the free energy F, and then use the solution of ∂F/∂x = 0 to determine x.At such an x, and indeed at all generic values of x, the spectrum of Q r (ω) turns out to have a gap.Earlier work [27,28,50] has advocated use of a 'marginal stability' condition, which leads to precisely the value of x for which the spectrum is gapless.We find that such a gapless spectrum obeys (1.9), and show a plot of Im Q r (ω) for all frequencies in Fig. 3.The lower plateau at q 0 in Fig. 2 vanishes as h → 0 as q 0 ∼ h 2 , unlike (1.11) for the Ising model.The value of x still vanishes linearly with T , as in the Ising model.We note that the p = 2 spherical model has no replica symmetry breaking, and the rotor spectrum is always gapped for h ̸ = 0.

II. ISING MODEL A. Landau action
We where ν n is a Matsubara frequency, ∆ is a small energy gap, and Λ is a high cutoff frequency.In the spin glass phase of the spherical model, the solution is gapless and replica symmetric: where q EA is the Edwards-Anderson order parameter.
We shall be interested in extending this solution beyond the spherical limit to the Ising model-at low frequency scales in the vicinity of the quantum critical point, where |ν n |, ∆, T ≪ Λ, q EA is small-and allow for replica symmetry breaking.In this regime, we can approximate the paramagnetic solution of the spherical model by where A, B are positive constants, while the spin glass solution is We now notice that if we shift Q by a frequency-independent and replica-diagonal constant, then the shifted Q ab (iν n ) becomes small at the relevant low-frequency scales on both sides of the quantum critical point.This makes the shifted Q ab (iν n ) a suitable field in which to carry out the Landau expansion of the Ising model.In terms of the original bilocal field, the shift is and the constant A characterizes nonuniversal, short-time physics not of interest to us.
Working with this shifted field, the important low-order terms in the Landau expansion for the action of the Ising model are [18] A where h is the longitudinal field, χ hb is a background contribution to the linear spin susceptibility, and the h = 0 Landau action is see Appendix A of Ref. 18 for a derivation of Eq. (2.8) from (1.1).Here, r is the parameter which tunes across the spin glass transition at h = 0: it is analogous to the coupling g in the Ising Hamiltonian in (1.1).The cubic term κ is analogous to the cubic term in Parisi's original theory of classical spin glasses [48,52].The U term only involves a single replica, and is a quantum self-interaction of the soft Ising order parameter.The y term is analogous to a quartic term in the classical case [48], where it is the term responsible for full replica symmetry breaking in the spin glass phase; however, the dynamic quantum effects of y are weak, and can be treated perturbatively.
Note that (2.8) does not contain the allowed quadratic term dτ We have removed such a term by exploiting the shift in (2.6).We will see below that such a choice is equivalent to the requirement for the validity of the Landau theory that the full function Q ab (iν n ) is small near the quantum critical point.This shift strategy is analogous to that followed for the critical theory of the Yang-Lee edge singularity [53,54].
Our analysis of the physics of the action A h will examine the behavior as a function of the tuning parameter across the spin glass transition r, the longitudinal field h, and the temperature T .We will keep the couplings κ and U of order unity and obtain results for small y.We find that for h > 0 there are contributions nonanalytic in y that are important to include; however, analytic corrections in integer powers of y are not crucial, and these are relegated to Appendix A.
The analysis of thermodynamic properties in the spin glass phase in Ref. 18 was carried out with a vanishing coefficient of the quartic term, y = 0: in this case the order parameter has replica symmetry.
Here, we will extend the solution to small y ̸ = 0, and show that the solution has broken replica symmetry.
We also show that the spin fluctuation spectrum remains gapless both at y = 0 and y ̸ = 0.

B. Free energy
Inserting the time-translational symmetric ansatz of (1.7) into (2.7) and (2.8), we obtain the free energy where the h-independent free energy is (2.10) We now insert (1.8) into (2.10) and obtain the free energy which consists of two pieces: a 'spin glass' and a 'quantum' component.
The first 'spin glass' component in (2.11) is given by where the h-independent terms are with The terms involving q ab in (2.13) are identical to those the Landau theory of the classical spin glass in Section 3.4 of Ref. 48.
The second 'quantum' component in (2.11) is h-independent (2.15) We will solve the saddle-point equations for F sg exactly, while those for F Q can be solved order-by-order in y.We will present the y 0 solution below, while the y 1 solution is presented in Appendix A. It turns out that the solution of F sg contains terms nonperturbative in y for h ̸ = 0, so it is important to treat the spin glass terms exactly.

C. Saddle-point equations
The saddle-point equations are most easily determined by taking the derivative of (2.10) with respect to Q ab (iν n ): The replica off-diagonal equation of (2.16) is while the replica diagonal part gives (2.18)

D. Zero-field limit
In this section, we first rederive the results obtained earlier [18,36] at h = 0, which we will then contrast with the case for a nonzero field later in Sec.II E.

Quantum paramagnet
In the paramagnetic phase, with q EA = 0 and q ab = 0, the spin glass free energy is F sg = 0.At h = 0, F sg,h = 0 wherefore only the quantum component survives, and Eq.(2.18) simplifies to an equation for This can be solved iteratively in powers of y.At order y 0 , the solution agrees with the form in (2.3) where the gap ∆ is given by the solution of and we have introduced a frequency cutoff to make the frequency summation finite.We will also need the cutoff at higher orders in y, but the low-frequency dynamics should remain cutoff-independent.The q(u) < l a t e x i t s h a 1 _ b a s e 6 4 = " x X G h e 7 3 j R K t s 5 3 g 8 X q 9 9 L R v 7   critical point r = r c0 is obtained by setting ∆ = 0, and is given by The order y 1 corrections to the saddle point appear in Appendix A 1 a.
a. Free energy: As mentioned above, for the case of the paramagnet, we have F sg = 0 for the spin glass component of the free energy in (2.11).The free energy F Q in (2.11) at order y 0 is where ∆ is the solution of (2.21).The order y 1 correction to F Q appears in Appendix A 1 a.

Spin glass
We begin by solving (2.17) at h = 0, which is the same equation of state as that obtained for the classical spin glass.In terms of the Parisi function q(u) characterizing the n → 0 limit of the replica matrix q ab , we can write (2.17) at h = 0 as (see Appendix E) This is the same as Eq.(3.74) in Ref.48 for the classical case.A replica-symmetric solution with q(u) = constant is possible, but this is not the preferred solution at h = 0 (as in the classical case, and, as we will see in Sec.II E, for the quantum h ̸ = 0 case).So, we only consider the replica-symmetry-breaking solution here, as shown in Fig. 4: with EA , x = 2yq EA /(βκ) . (2.27) It can now be verified that the term proportional to δ νn,0 in (2.18) vanishes We can therefore conclude that the complete saddle-point equations for Q r (iν n ) and q(u) reduce to the following three equations for Q r (iν n ), q EA and x: ) with α > 0.Then, at T = 0, the terms of order Q 2 r and Q 3 r in (2.29) will have imaginary parts which vanish faster than |ω| α as |ω| → 0. Collecting all terms of order |ω| α in the imaginary part of (2.29), we obtain the condition which is automatically satisfied from (2.30).So, the spin dynamics are gapless at all values of y in the spin glass phase with h = 0. We will see in the explicit solution below that the exponent α = 1, so that ImQ r (ω → 0) ∼ ω.
b. Solution: The equations (2.29-2.31)can be solved iteratively in powers of y.At order y 0 , we have The order y 1 corrections to the saddle point appear in Appendix A 1 b.
c. Free energy: Inserting the solution (2.25)-(2.27)back into (2.13), and using the relation (see we obtain This is the full expression for F sg in (2.11), valid to all orders in y in terms of the exact q EA .However, q EA is known only to order y 1 in (A8).
As in the paramagnet, we will only determine F Q in (2.11) order-by-order in y.At order y 0 , the value (2.37) The order y 1 correction to F Q appears in Appendix A 1 b.

Phase diagram
At order y 0 , the system is the spin glass phase for r < r 0c , where r c0 is given by Eq. (2.22).In the classical limit, T ≫ Λ, we need only consider the ν n = 0 term, so the phase boundary is at r = 0.In the quantum limit T ≪ Λ, we can evaluate the summation over Matsubara frequencies using the identity for odd spectral functions ImD(Ω) = −ImD(−Ω).Then, we obtain the phase boundary at (2.39)

E. Nonzero field
Having reviewed the results in the absence of a longitudinal field above, we now address the h ̸ = 0 case and determine its phase diagram.

Replica-symmetric solution
Due to the longitudinal magnetic field, there is an average moment on each site, so q EA is always nonzero; therefore, the replica-symmetric solution is nothing but the paramagnet.For q ab = q EA (1 − δ ab ), Eq. (2.17), together with (2.14), reads Resultantly, the term proportional to δ νn,0 in (2.18) vanishes so that the equation for Q r (iν n ) remains unchanged from that in (2.29).The gapless condition in (2.33) is not satisfied by (2.40), even though q EA ̸ = 0. Thus, the solution for Q r (iν) will have a gap.
At order y 0 , the solution for Q r0 (iν n ) from (2.18) has the same form as (2.20), but the equation for the gap ∆ in (2.21) is now modified to [18] The solutions to Eq. (2.41) are shown in Fig. 5.The order y 1 corrections to the saddle point appear in Appendix A 2 a.
a. Free energy: Inserting the solution for Q r (0) in (2.40) into (2.12)we obtain the spin glass component of the free energy (2.43) Note that we need to insert the solution for q EA in (A15) into (2.43).
For the quantum component, the order y 0 contribution to the free energy in (2.15) is which reduces to (2.23) at h = 0.The order y 1 correction to F Q appears in Appendix A 2 a.
Gap as a function of (r,h)

Replica symmetry breaking
From (2.17), the equation (2.24) for the Parisi function is modified to The solution of (2.45) is modified from (2.25) to that shown in Fig. 1 q q EA u/x, (q h /q EA )x < u < x q EA , x < u < 1 . (2.46) The function q(u) is nonanalytic at u = (q h /q EA )x and u = x, but has no discontinuities.The values of q EA and x are unchanged from those in (2.27), and the value of q h is in (1.11).The result in (1.11) is the origin of the nonanalytic dependence on y.
It can now be verified that the term proportional to δ νn,0 in (2.18) vanishes, and the equations for Q r (iν n ), x, q EA remain unchanged from those in (2.29)-(2.31).Indeed, the only change from the h = 0 solution in Sec.II D 2 is in the form of q(u) for x < (q h /q EA )x in (2.46).
Hence, the gapless condition in (2.33) is now satisfied, and the solutions for Q r (iν n ) and q EA remain unchanged from that for the gapless spectrum in Sec.II D 2 b.
a. Free energy: Inserting the solution (2.46), (1.11) into (2.12),we obtain the extension of (2.36) to nonzero h As the introduction of h does not modify the values of Q r (iν n ) and q EA , the free energy F Q remains the same as that in Sec.II D 2 c.

Phase diagram
To obtain the phase boundary between the replica-symmetric and replica-symmetry-breaking phases, we compare their free energies at leading nontrivial order in y.
For the gapped replica-symmetric phase we have the free energy given by (2.42)-(2.44) where the gap ∆ is given by the solution of (2.41), For the gapless replica-symmetry-breaking phase we have from (2.37) and (2.47) The summations in (2.41), (2.48), and (2.49) can be evaluated for T, ∆ ≪ Λ using identities obtained from (2.38): e Ω/T − 1 . (2.50) The phase diagram so obtained in shown in Fig. 6; note that due to the n → 0 limit, we have to choose the phase with the maximum free energy [49].We observe that the extent of the RSB phase in parameter space, which occurs for r < 0, shrinks with increasing h and as y is reduced.Note that ∆ remains finite, albeit small, at the transition point for h > 0 (it vanishes at the transition for h = 0), indicating a first-order transition.

III. QUANTUM SPHERICAL p-ROTOR MODEL
We now turn to the analysis of the p-rotor model described in Eq. (1.3).The classical infinite-range spin glass with p-spin interactions, but without any spherical constraints, was originally introduced by Derrida [55,56].Gross and Mézard [43] first studied the generalization of this model to nonzero magnetic fields, obtaining an exact solution for p → ∞.The classical spherical model [46] can also be solved for any finite p, including in the presence of an external field [45,57].Quantum extensions of these models, by adding a noncommuting transverse field, have been investigated both with (for p = 3) [27,29] and without (for p → ∞) [58] the supplementary spherical constraint.However, the problem of the quantum model in a longitudinal field, which we examine next, has remained unexplored so far.

A. Effective action
To begin, we derive the effective action for the quantum spherical p-rotor model, which will then form the basis for our subsequent saddle-point calculations.In the path integral, the spherical constraint (1.4) can be enforced using the exponential representation of the Dirac delta function which is then inserted into Eq.( 1.3) at the expense of introducing an auxiliary field z(τ ).
At this stage, since the disorder is quenched, we average the free energy log Z[J i 1 ...ip ] over disorderusing the replica trick-instead of the partition function Z[J i 1 ...ip ] itself (which would correspond to the annealed average).The replicated partition function is given by where a = 1, . . ., n is the replica index.Now, we can perform the disorder average by evaluating simple Gaussian integrals to find Note that the overlap between two different replicas of the system where the overline denotes an average over disorder realizations, appears naturally in Eq. (3.3).We then insert the identity into the partition function, which, in terms of the collective variables Q ab , λ ab , reads Here, we have replaced bilinears of σ i with Q ab wherever possible and collected in the first two lines all the remaining terms involving the σ i fields, which we will integrate out in the next step.Using, for concreteness, the convention for n-dimensional vectors v, j and an n × n matrix M, we obtain with the matrix G defined in replica and imaginary-time space as The replicated effective action can thus be extracted as One of the advantages of the replica approach we are now positioned to harness is that since the effective action (3.10) is proportional to N , the saddle-point approximation is exact in the limit N → ∞.

B. Saddle-point equations
We will only be interested in the saddle point of the theory (3.10), in which case we can assume that all fields are time-translation invariant and transform to Matsubara frequency space, as in (1.7).Then, we have We now employ a slight change of notation to obtain expressions similar to those in Ref. 51, defining in terms of which, the effective action is Now, we use the identity (E21) to absorb the h 2 term into the log det as Then, the saddle-point equation with respect to Σ ab (ω n ) is Inserting this back into (3.14),we obtain an effective action just for Q ab and λ: From this, the saddle-point equations for Q ab and λ can be read off as with the Fourier transform defined by and similarly for Σ ab (ω n ).We numerically study these equations for the replica-symmetric and onestep replica-symmetry-breaking cases for specific values of p.In the discussion below, we focus on the p = 3 model since, as we show below, it exhibits an interesting and nontrivial phase diagram due to the competition between these solutions.
For completeness, the solution of the p = 2 model is presented in Appendix C. Note that with p = 2 and h = 0, Eq. (3.17) reduce to (33.32) and (33.33) in [51].Moreover, for p = 4, h = 0 the equations (3.17) are very similar to Eqs. (S50-S53) in Ref. 59, the main difference being that ω 2 n is replaced by iω n .

C. Replica-symmetric solution
In this section, we numerically determine the replica-symmetric (RS) solution of the p = 3 model as a function of transverse and longitudinal fields (g and h, respectively), while choosing the temperature to be close to zero.
FIG. 7: Behavior of the order parameter q for the replica-symmetric solution of the quantum spherical p = 3 model as a function of the transverse field g for different values of longitudinal fields h (with J = 1, T = 0.005).The phase transition between the replica-symmetric and paramagnetic solutions occurs only for h = 0 and is continuous.
The ansatz (1.8) for the replica overlap functions in the RS case necessarily has a replica off-diagonal component since, for nonzero h, the static moment is always nonzero.Then, using the identities in Appendix E, the equations in (3.17) become (for p = 3) We solve these equations self-consistently in the imaginary-frequency domain for the order parameter q as a function of the transverse field g at low temperatures for several values of the longitudinal field, as shown in Fig. 7.In this model, the phase transition between the RS phase (q ̸ = 0) and the paramagnet (PM, q = 0) occurs only at zero longitudinal field.With even a slight increase in h, the replica-symmetric solution becomes equivalent to the paramagnet for all values of g, as in Sec.II E 1.
Regardless of the value of the longitudinal field, we can show that the spectrum of the RS solution is gapped.To do so, we compute the dynamic spin susceptibility of the system, which requires analytic continuation of the equations (3.20)-(3.24) to real frequencies.Taking the zero-temperature limit, we obtain equations of the form ) where the spectral representation of the regular component of the replica overlap and Σ ′ r (ω) and Σ ′′ r (ω) are the real and imaginary parts of Σ r (iω n ), respectively.The spectral function is related to the retarded Green's function as ρ(ω) = Im Q r (ω)/π.Taking the z → ∞ limit above, the associated sum rule becomes +∞ 0 dω ωρ(ω) = g/2.
We note that the equations above, strictly speaking, permit several solutions, one of which has a discontinuity at zero frequency given by ρ(ω = 0) = 1/(2πJ √ 6q).To avoid this, we initialize the system with a simple step function with a small gap and let the equations converge (defined as when the error reaches 10 −20 ).As a final step, we also check that the sum rule is satisfied.
The spectral functions obtained in this fashion are presented in Fig. 8.For h = 0, we observe that as g increases, the gap saturates to a finite value and eventually, after the transition point, the system becomes a quantum-disordered paramagnet, as shown in Appendix D. However, at a finite field, the gap increases with increasing g without a visible saturation since the replica-symmetric solution and the paramagnet are one and the same.

D. One-step replica symmetry breaking
In this section, we now allow for nontrivial structure in the replica off-diagonal space.In general, replacing Eq. (3.19), we can write The matrices q ab and ϱ ab are characterized by Parisi functions q(u) and ϱ(u) with u ∈ [0, 1].For the diagonal components, we make the same ansatz as in (3.19) We consider the one-step replica-symmetry-breaking (RSB) ansatz, for which (see Fig. 2) Using the identities in Appendix E, the equations in (3.17) now become ) ) ) ) where we have simplified the equations by changing variables and introducing m 2 = λ − Σ r (ω n = 0).
There are eight unknowns, namely, m 2 , Q r , Σ r , q 1 , q 0 , ϱ 1 , ϱ 0 , and βx, in the equations above.All but one of these variables can be determined by solving the seven equations (3.35)- (3.41).However, the value of the breakpoint βx is undetermined and can be fixed either by demanding a gapless solution, as in Appendix 3 of Ref. 59, or by requiring the free energy to be stationary with respect to x.The additional equation then makes the system complete and the solution can be obtained by solving the set of equations self-consistently.In the following, we derive this extra equation in both cases: by imposing the gapless constraint, and by using the free energy.

Gapless condition
Let us first determine the condition for a gapless spectrum in the RSB case; our approach here is similar to the analysis in Ref. 59.For a gapless spectrum, we expect the low-frequency expansion of Q r , for real ω, to have the form for some α > 0, b > 0. From (3.36), we find that the leading singularity in Σ r (ω) is given by the linear-in-Q r (ω) term in (3.36), so However, from (3.39), for low frequencies, we also have Therefore, comparing equations (3.43) and (3.44), we arrive at the constraint on the variable m that has to be satisfied for the gapless solution to exist: It is easy to check that for the special case of p = 2, this expression precisely agrees with Eq. (33.41)It is important to note that upon solving the replica-symmetric equations self-consistently, we did not obtain a gapless solution.One way to understand this is to plug the gapless constraint (3.45) into the equations for the replica-symmetric case (3.20)- (3.24), which expressly demonstrates that the constraint cannot be satisfied for any values of h.Therefore, the replica-symmetric solution always has a gap.

Free energy
Another way to fix x is to determine the saddle point of the free energy with respect to x.To do so, we rewrite the effective action of the p-rotor model (3.16) explicitly for the one-step replica-symmetry-breaking case This action allows one to directly verify that the saddle-point equations are in fact correct: taking the derivative of (3.46) with respect to q 0 and λ, we see that the equations (3.35) and (3.41) indeed hold after a change of variables to m 2 = λ − Σ r (ω n = 0).
To obtain an equation for x, we differentiate the action with respect to the breakpoint and obtain Although such a saddle-point condition has been employed to obtain the equilibrium state for the classical model [47], we will not use (3.47) here.It what follows, we instead use the gapless constraint to determine the breakpoint x.

Spectral functions
As for the replica-symmetric case before, we are interested in the behavior of the spectral functions in the RSB case.To this end, we analytically continue the equations on the imaginary-frequency axis (3. ) ) ) (3.55) For the final equation, as mentioned, we use the gapless constraint derived above in (3.45), and this gapless behavior should also be reflected in the spectral functions.
Therefore, we are allowed to make a linear ansatz We insert this ansatz in Eqs.(3.49)-(3.55)and solve the equations self-consistently for f (ω).It is quite challenging to obtain the self-consistent solution to these equations for all unknowns.Instead, we simplify this task by solving the equations on the imaginary-frequency axis and self-consistently obtain the value of the order parameter q 1 at a temperature close to zero with a high precision.Then, we use this result to obtain the behavior of the spectral functions.To ensure that the converged solution is meaningful, we check the sum rule, which, in this case, is given by We sketch the behavior of the spectral functions Im Q r (ω) in Fig. 3. From the equations above, it is clear that spectral functions are independent of h, βx and q 0 .Only two parameters, namely, q 0 and βx, depend on h.This is similar to the behavior of the Ising model in Sec.II E 2, where the only change from the h = 0 solution in Sec.II D 2 was in the form of q(u) for u < (q h /q EA )x in (2.46).It is also instructive to determine the dependence of q 0 and βx on h in the h → 0 limit, which can be calculated from the equations (3.54) and (3.55).We find that, to the lowest order, the dependence is simply which, as noted in Sec.I A, is in distinction to the h 2/3 scaling observed for the Ising spin glass.FIG.9: (a) Behavior of q 0 and q 1 as a function of g for the spherical quantum p = 3-rotor model in a longitudinal field with T = 0.005 and J = 1.Note that q 1 is independent of h but q 0 is not.(b) Boundary of stability of the RSB phase as determined by setting q 0 (g, h * ) = q 1 (g).For h > h * (g), the necessary condition q 0 ≤ q 1 required for the RSB solution no longer holds.

E. Phase diagram
Finally, having assembled all the necessary ingredients, we now compute the phase diagram of the quantum spherical p = 3-rotor model in a longitudinal field at low temperatures.To do so, we first determine the boundary of stability of the RSB solution by self-consistently computing the spin glass order parameter from equations (3.35)-(3.41)and comparing q 0 and q 1 .For the stability of the RSB solution, one has to require q 0 ≤ q 1 [46], and the solution ceases to be physical when q 0 > q 1 .Our numerical results in this regard are displayed in Fig. 9.In particular, Fig. 9(a) presents q 0 and q 1 as a function of g for various values of h; the RSB solution exists only for q 0 that lie below the line q 1 (g).Based on this information, Fig. 9(b) traces, in g-h space, the line h * (g) along which q 0 (g, h * ) = q 1 (g).Since this marks the boundary of stability of the RSB solution, the phase transition between the RSB and RS phases may occur at this line.
To check whether this is indeed the case, we compute the free energies of both the RS and RSB solutions.The free energy in the RS case can be written using (3.16) and we obtain where we inserted a factor of (ω 2 n + Γ 2 )/g into the argument of the logarithm to regularize the behavior at large |ω n |, and subtracted a constant to regulate the frequency summation in the term linear in Q r .This is equivalent to normalizing the path integral with an oscillator of frequency Γ.In order to make the resulting expression independent of Γ, we subtracted the free energy of this oscillator.For the RSB case, we use the expression for the free energy computed above, focusing on p = 3, Here, we use the same regularization as for the free energy of the replica-symmetric solution.In both cases, Γ has to be much smaller than the Matsubara frequency cutoff.
We evaluate the difference in free energies between the RSB and RS solutions to find the parameter regions where each phase is dominant.Following the discussion in Ref. 49, we pick the ground state as the one that has a larger free energy.The resulting phase diagram is shown in Fig. 10; we find that the phase boundary obtained by comparing the free energies coincides (to a precision of |F rsb eff − F rs eff |/(N n) < 10 −5 ) with the line h * (g) demarcating the stability of the RSB solution in Fig. 9(b).

IV. CONCLUSIONS
Motivated by the quantum simulation of spatially disordered long-ranged interacting systems for solving combinatorial optimization problems [1,60], we have studied the equilibrium dynamic properties of two quantum spin glass models with infinite-range interactions: the Ising model in (1.1), and the spherical p-rotor model in (1.3) and (1.4).In addition to random couplings between spins/rotors, both models have an applied longitudinal field h, and do not break any global symmetry in any phase.A coupling g tunes the strength of the quantum fluctuations, and at large g, we obtain a gapped paramagnetic phase.Our primary interest was in the small-g regime, where replica symmetry breaking leads to a quantum spin glass phase.
For the Ising model, we employed the Landau theory approach of Ref. 18 to show that the spin glass phase has full replica symmetry breaking, with the Parisi function sketched in Fig. 1.We computed the spin autocorrelation function in this spin glass phase, and showed that it is generically gapless, with a spectral density which vanishes linearly with frequency (see (1.9) and Secs.II D 2 b and II E 2).
In contrast, the spherical p-rotor model (p ≥ 3) has a different behavior in the spin glass phase.There is only one-step replica symmetry breaking (see Fig. 2), and the breakpoint x remains undetermined by the saddle-point equations for the matrix spin glass order parameter.Rigorous mathematical work on the classical g = 0 model [47] has shown that the equilibrium value of x is determined by the stationarity of the free energy with respect to x.For this x, we showed that the spin glass state has an energy gap.In contrast, if a marginal stability condition [27,28,50] is used to determine x, we found a gapless spectrum with a linear frequency dependence at small frequency (see Fig. 3), similar to that of the Ising model.
Our results on the nature of the spin glass phase also find relevance in the context of quantum optimization problems since the properties of the gap-or lack thereof-directly inform the feasibility of preparing low-energy states via adiabatic algorithms [61][62][63].
Looking ahead, it would be useful to obtain the results described here without the replica method, using a path integral on the Schwinger-Keldysh contour.This has already been studied for the spherical p-rotor model [28,37,64], where the breakpoint of one-step replica symmetry breaking, x, has been connected to an effective temperature of the aging dynamics.Moreover, the equation for a gapless spectrum in (3.45) is the same as Eq.(6.24) obtained in Ref. 28 from aging dynamics.However, the aging dynamics of the quantum Ising model have not been studied so far, and it would be interesting to relate that to the full replica symmetry breaking characterized by Fig. 1.Such a computation would be the analog of earlier studies [65][66][67][68][69] of aging dynamics in classical spin glasses described by Langevin equations, and we expect that there would be significant differences [69,70] between the Ising model and the spherical model in the quantum problem too.ϱ = gJ 2 q, (C2) (C4) Specifically, the equation for Q r (ω n ) can be written as Then, the order parameter is given by From the equation above, we observe that the condition λ ≥ 2gJ has to be satisfied for the order parameter to be a real number.Note that for h = 0, we find which is consistent with Eq. (33.41) of Ref. 51.For h ̸ = 0, the equation for the order parameter can be formulated as We can simplify the expression above further to where we introduce the notation We note that the order parameter has to be positive, q > 0, so λ > λ * .In addition, we require λ > 2gJ such that λ * > 0.
Finally, we rewrite the correlator Q r (ω n ) in terms of our new notation as The system is gapless in the limit h → 0. (b) The order parameter q as a function of the transverse magnetic field for several values of h, illustrating the absence of a sharp phase transition at any nonzero g and h.
Since λ * > 0, there is always a gap in the spectrum.The gap is given by the value of the frequency at which the square root in the Green's function (analytically continued to the real frequency) changes sign.
This value is simply which is always greater than zero in the presence of a nonzero field h.The resulting behavior of the gap as a function of the longitudinal field is shown in Fig. 12(a).
To find the value of the spin glass order parameter, we need to find λ according to , (C15) The numerical solution for q as a function of g is shown in Fig. 12 of order y 1 in F Q for the Ising model 1.Zero longitudinal field a. Paramagnet b.Spin glass I. INTRODUCTION r F y 5 O y 1 V r 7 I 4 8 n A A h 3 A M H p x D F W 6 h B n U g 8 A D P 8 A p v j n R e n H f n Y 9 a a c 7 K Z f f g D 5 / M H c I K P D g = = < / l a t e x i t > q EA < l a t e x i t s h a 1 _ b a s e 6 4 = " Z p B P a e t w c m o + F v T H k o o b 6 1 5 Q S I U = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 e I 5 g H J E m Y n v c m Q 2 d l 1 Z l Y I S z 7 B i w d F v P p F 3 v w b J 8 k e N F r Q U F R 1 0 9 0 V J I J r 4 7 p f T m F p e W V 1 r b h e 2 t j c 2 t 4 p 7 + 4 1 d

FIG. 3 :
FIG. 3: Behavior of the zero-temperature spectral function ρ(ω) = ImQ r (ω)/π, as determined by Eqs.(3.49)-(3.56),for the one-step replica-symmetry-breaking solution of the p = 3 spherical quantum p-rotor model at different values of the transverse field g (with J = 1).The marginal stability condition (3.45) is imposed.The spectrum is independent of h as long as we are in the replica-symmetry-breaking phase in Fig. 10.
motivate the structure of the Landau action for the Ising model by recalling the solution of the spherical model for the p = 2 case [17], in which the interactions have the same form as in the Ising model, but the local Ising constraint has been replaced by the global constraint in (1.4).At h = 0, the large-N solution in the paramagnetic phase is of the form (see Appendix C, and Chapter 33 in Ref. 51)

5 9 3 1 <
5 m L e u O P n M E f y B 8 / k D f d 2 M v w = = < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " Z w K k a J c d U 3 h H l n h M y V D 4 t d z 5 5 Q M = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H a N r y P R i 0 d I 5 J H A h s w O D Y z M z m 5 m Z o 1 k w x d 4 8 a A x X v 0 k b / 6 N s s j j w c w T G c g g d X U I E 7 q E I d G C A 8 w y u 8 O Q / O i / P u f M x b c 0 4 2 c w h / 4 H z + A O l 5 j Q Y = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " z M m D t Y / I j r P 2 9 O r C M 5 b t p f 6 U i u 0

FIG. 4 :
FIG.4:The Parisi spin glass order parameter for the Ising model at h = 0.

8 FIG. 5 :
FIG.5:The gap ∆ of the replica-symmetric solution of the Ising model as a function of r and h for various temperatures, as determined from (2.41).The parameters U , κ, and Λ have all been set to unity in these calculations.

FIG. 6 :
FIG. 6: Phase diagram of the Ising model in the (r, h) plane for various values of y at T = 0.2, illustrating the replica-symmetric (RS) and replica-symmetry-breaking (RSB) phases.The value of ∆ at each point in parameter space is obtained from Fig. 5, and U, κ, Λ have been set to unity, as previously.

FIG. 8 :
FIG.8: Spectral functions of the quantum spherical p = 3 model for (a) h = 0 and (b) h = 1 at zero temperature (with J = 1).(a) The phase transition between the RS phase and the quantum paramagnet occurs at a finite g > 3, and the gap thus eventually saturates upon increasing g.(b) The replica-symmetric solution persists for all values of g with a gap that grows with increasing g.

FIG. 10 :
FIG. 10: Phase diagram of the spherical quantum p = 3-rotor model in a longitudinal field with T = 0.005 and J = 1.The phase diagram is evaluated by calculating the difference between the free energies of the replica-symmetric (RS, 3.59) and replica-symmetry-breaking (RSB, 3.60) solutions.The orange line on the phase boundary traces the function h * (g) plotted in Fig. 9(b).

)FIG. 12 :
FIG.12:(a) Gap in the spectrum (C13) as a function of the longitudinal field strength h for the quantum spherical p = 2-rotor model, taking g = 1.The system is gapless in the limit h → 0. (b) The order parameter q as a function of the transverse magnetic field for several values of h, illustrating the absence of a sharp phase transition at any nonzero g and h.
in Ref. 51 and is similar to Eq. (4.11) in Ref. 29.