Strongly coupled quantum phonon fluid in a solvable model

We study a model of a large number of strongly coupled phonons that can be viewed as a bosonic variant of the Sachdev-Ye-Kitaev model. We determine the phase diagram of the model which consists of a glass phase and a disordered phase, with a first-order phase transition separating them. We compute the specific heat of the disordered phase, with which we diagnose the high-temperature crossover to the classical limit. We further study the real-time dynamics of the disordered phase, where we identify three dynamical regimes as a function of temperature. Low temperatures are associated with a semiclassical regime, where the phonons can be described as long-lived normal modes. High-temperatures are associated with the classical limit of the model. For a large region in parameter space, we identify an intermediate-temperatures regime, where the phonon lifetime is of the order of the Planckian time scale $\hbar/k_B T$.


I. INTRODUCTION
Is there a fundamental limit to how fast can a quantum many-body system relax back to equilibrium? This basic question arises frequently in the interpretation of experiments in condensed matter systems [1][2][3][4][5][6]. In particular, transport coefficients can be related to relaxation times of electrical and thermal currents [7][8][9][10][11][12]. Often, the current relaxation times are tied to extrinsic mechanisms, such as the momentum loss rate due to impurity scattering. However, in setups where the bottleneck for current relaxation is the intrinsic thermalization time in the system, it has been proposed that the relaxation time has to obey a fundamental 'Planckian' bound, τ th ≥ α kB T , where α is an unknown constant of order unity [13].
Evidence for this intriguing idea comes both from solvable models, such as holographic theories and systems near quantum critical points [2,[14][15][16][17][18][19][20], and from experiments [3][4][5]. Much attention has been devoted to electrical transport in 'strange metals', where the Planckian bound is a natural way to explain the linear dependence of the resistivity on temperature. To test this hypothesis, the transport lifetime can be estimated using the Drude formula for the resistivity: ρ = m * ne 2 τ , in materials where the electronic effective mass m * and density n are well known. This procedure indeed yields τ ∼ kB T , with a coefficient of order one, in a host of different materials in regimes where ρ ∼ T [3][4][5].
Surprisingly, it was recently noted that in a wide class of insulating compounds at high temperature, similar physics may be at play [21][22][23]. In these systems (for example, complex oxides like SrTiO 3 ), the thermal current is carried by lattice vibrations. Around room temperature and above, the thermal conductivity κ is approximately inversely proportional to temperature [24]. Defining the thermal transport lifetime as τ = D th /v 2 ph , where D th is the thermal diffusivity (directly measured or obtained via the Einstein relation κ = cD th , where c is the specific heat) and v ph is a characteristic phonon velocity, operationally defined as the averaged speed of sound, gives again τ = α kB T where α is found to be in the range 1-3 in a host of different poorly thermally conducting materials. In contrast, in good thermal conductors α is much larger. For example, diamond (which exhibits κ ∼ 1/T over a range of temperatures explained as a result of phonon umklapp), α is found to be ∼ 50 [23].
This observation is particularly counter-intuitive since, at such elevated temperatures, one would naively expect the lattice dynamics to be essentially classical. In a classical analysis of the lattice dynamics, any time scale must be proportional to √ M , where M is the ion mass. A simple dimensional analysis argument implies that τ ph ∼ M/m where m is the electron mass [25]. The observation of a short transport lifetime has the intriguing implication that in a broad temperature regime, the system should be thought of as a quantum-mechanical 'fluid' of strongly coupled lattice vibrations, rather than in terms of individual long-lived phonon excitations. This view is supported by the fact that the estimated phonon mean free path is very short, of the order of the lattice spacing or less [23]. Furthermore, at least some of the complex oxide materials that show α ∼ 1 have very high frequency optical phonon branches, far exceeding room temperature, potentially explaining why the lattice dynamics is not fully classical even at room temperature and above.
These intriguing observations call for a theoretical framework where the crossover from classical to quantum dynamics of lattice vibrations can be investigated. Considering a generic system of coupled non-linear oscillators, we expect that at sufficiently low temperatures, the system is always described in terms of low-energy, longlived normal modes (or phonon quasi-particles). Conversely, at sufficiently high temperature, the dynamics is expected to become classical. The quantum relaxation time τ ∼ /k B T may thus appear only at intermediate temperature scales. 1 This regime is the most difficult to analyze theoretically.
In this work, we propose a simple model of strongly interacting phonons, that can be used to address the above questions. The model can be viewed as a bosonic variant of the Sachdev-Ye-Kitaev (SYK) model [27][28][29] of N degrees of freedom coupled via an all-to-all, random interaction; similarly to SYK, the model is solvable in the large N limit. The real-time dynamics of the model is found to follow the trends described above, with long relaxation times associated with long-lived phonon modes at low T , a crossover to classical nonlinear dynamics with τ ∼ √ M at high T , and a broad intermediate T regime where the lifetime is of the order of the Planckian time scale /k B T . Some representative results for the phonon lifetime τ ph scaled by /k B T as a function of temperature are shown in Fig. 1a,b.
This paper is organized as follows. In Section II we introduce the model and identify the relevant energy scales. In Section III we discuss some of its thermodynamic properties. We map out the phase diagram of the model, and in addition we discuss the specific heat in the disordered phase. The dynamics of the model, and in particular the phonon lifetime as a probe to identify different dynamical regimes, are discussed in Section IV. In Section V we discuss a generalized version of the model with multiple phonon branches. In Section VI we comment on several differences between the fermionic SYK model and its bosonic variant. Details on the imaginary-and realtime derivations are given respectively in Appendices A and B, and Appendix C contains details on the numerical methods.

II. MODEL
We consider a system of N coupled nonlinear oscillators in d = 0 + 1 space-time dimensions, described by the Hamiltonian Here, φ i is the displacement of the ith mode, π i is the conjugate momentum (such that [φ i , π j ] = i δ ij ), M is the mass, and Ω i is the frequency in the absence of nonlinearity. We will begin with the case where the frequencies are all the same, Ω i = Ω 0 , considering a more general situation later. The cubic couplingsṽ ijk are chosen to be independent random Gaussian variables, each satisfying v ijk = 0 andṽ 2 ijk = 2ṽ 2 (no sum), where (·) denotes averaging over realizations ofṽ ijk . The quartic interactioñ u > 0 stabilizes the system (whenũ = 0, the energy is not bounded from below, due to the cubic term). The Hamiltonian in Eqn. 1 is similar to the spherical p-spin model [30], studied in the context of quantum spin glasses; the  figure. (b) presents the high-T behavior of the phonon lifetime for the same set of parameters as in (a). This regime is associated with the approach to the classical limit, where the linearity of T τ ph / suggests that the phonon lifetime becomes independent of T (see Eqn. 10). (c) shows the phonon spectral function A(ω) for four increasing temperatures. The dynamical crossover from the semiclassical regime (with T /Ωv = 0.25) to the 'Planckian' regime is demonstrated by the significant broadening of the spectral peak. Here, u/Ω 3 v = 1.4, and Ω * denotes the T → 0 renormalized phonon frequency (see Eqn. 9). difference is that in the p-spin model, 1 N i φ 2 i is constrained to unity, whereas here this quantity is unconstrained (theũ term implements a 'soft constraint' on the magnitude of 1 N i φ 2 i ). In addition to the energy scale Ω 0 , we can define two energy scales associated with the non-linear terms by dimensional analysis. Theṽ term defines an energy scale Ω v = 2 /M 3/5ṽ 2/5 , while theũ term is associated with the scale Ω u = 2 /M 2/3ũ 1/3 . We set = k B = 1 henceforth, unless stated otherwise. These energy scales become apparent in Lagrangian formulation under the where the rescaled couplings are given by u = Ω 3 u and v 2 ijk ≡ 2v 2 = 2Ω 5 v , such that the energy scales of the system are given by Ω 0 , v 2/5 and u 1/3 . In this work, we focus on the strong coupling regime of the model, where Ω 0 ∼ v 2/5 ∼ u 1/3 . We use the rescaled Lagrangian formulation henceforth.

III. THERMODYNAMICS AND PHASE DIAGRAM
In this section, we discuss some of the thermodynamic properties of the model. We study its phase diagram, and compute the specific heat of the disordered phase as a function of temperature. The main interest of this work is the existence of a 'phonon fluid' regime -a dynamical crossover region in parameter space where phonons are not well-defined quasiparticles. However, in this region, it might be thermodynamically favorable for the system to realize an ordered (glassy) phase that masks the 'phonon liquid' behavior. It is therefore crucial to map the phase diagram of the model and find the boundaries of the glass phase. In addition, we study the specific heat of the disordered phase that can serve as a simple thermodynamical diagnostic for the crossover from the quantum mechanical 'phonon fluid' to the classical regime.

A. Phase Diagram
Let us consider the simple version of the model, where Ω i = Ω 0 for all i = 1, ..., N , which we dub the 'singlebranch' (SB) model. We study its equilibrium phase diagram as a function of temperature T and frequency Ω 0 , for fixed values of v, u and M that satisfy Ω v ∼ Ω u .
We compute the disorder-averaged free-energy density within the framework of the replica method, where it is given by Here, f is the free-energy per mode, β is the inverse temperature and Z n is the replicated partition function of n ∈ N replicas of the model, where n is then analytically continued to zero. We proceed by introducing bilocal fields G αβ (τ, and enforce this identity with the Lagrange multiplier fields Π αβ (τ, τ ′ ), where α, β = 1, ..., n are the replica indices and imaginary-time arguments are denoted by τ, τ ′ ∈ [0, β]. This allows us to express (see Appendix A) the disorder-averaged replicated partition function as a functional integral of an effective action, In the limit of N → ∞, this functional integral is controlled by the saddle point of the effective action, δS eff /δA αβ = 0, A = G, Π, leading to a closed set of self-consistent equations. Different phases of the model are characterized by the replica-space structure of G αβ . Namely, a diagonal G αβ corresponds to the disordered phase of the model, while a G αβ with non-zero offdiagonal elements corresponds to an ordered (glassy) phase, where the off-diagonal elements of G αβ are the order parameters. The saddle-point equations for the off-diagonal components of G αβ are identical to those of [30], which implies that only two stable solutions exists in replica space, the diagonal solution and the one-step replica symmetry breaking (1SRSB) solution.
The diagonal solution is given by G αβ (τ, τ ′ ) ≡ G (τ − τ ′ ) δ αβ , and its corresponding self-consistent equations read Here, ω k = 2πk/β, k ∈ Z are the bosonic Matsubara frequencies, hats denote Matsubara-frequency domain functions, and we have used the fact that G is imaginary-time translationally invariant in thermal equilibrium. The one-step replica symmetry breaking (1SRSB) solution is defined as G αβ (τ, τ ′ ) ≡ (g d (τ − τ ′ ) − g EA ) δ αβ + (g EA − g 0 ) ǫ αβ + g 0 , where ǫ αβ = 1 if α and β are in a diagonal block of size m and ǫ αβ = 0 otherwise, and g EA is the so-called Edwards-Anderson order parameter [31,32]. The self-consistent equations and more details on the 1SRSB solution can be found in Appendix A.
One can solve these self-consistent equations numerically by an iterative procedure (see Appendix C). Substituting the solutions back in S eff enables us to obtain the free-energy density of the two phases of the model and thereby to obtain its phase diagram. For a fixed u, we find that the model realizes two phases, a disordered (replica diagonal) phase, and a glass (1SRSB) phase. The transition between the two phases is first order. This transition is associated with a discontinuity in the order parameter g EA while the break-point parameter m < 1. 2 In the limit of T → 0, we find that the break-point parameter m → 0 at the transition, which suggests that the replica symmetry is restored at the T = 0 quantum phase transition, similarly to Refs. [30,33]. The phase diagram in the (Ω 0 , T ) plane for fixed values of u and v is shown in Fig. 2. As can be seen in the figure, the region of the glass phase shrinks upon increasing u.

B. Specific Heat
Consider the specific heat c of the disordered phase in the SB model. The temperature dependence of c can be used as a thermodynamical diagnostic for crossover to the classical limit of the model. At high temperatures, we find that c saturates to a constant value according to an anharmonic variant of the Dulong-Petit law 2 This discontinuity might not be sufficient for a first-order transition if the break-point parameter m → 1 at the transition, because the effective number of degrees of freedom that are involved in the transition is (1 − m)g EA [30]. In our model, we find that 0 < m < 1 at the transition so that it is indeed first-order. (Eqn. 6). Conversely, temperatures for which there is a large variation in the value of c(T ) are associated with a non-classical behavior. We will use this simple diagnostic to further illustrate that the dynamical Planckian regime of the model, discussed in Sec. IV (see also Fig. 1), is indeed of quantum mechanical nature. The computation of the specific heat is done as follows. We first use the disorder-averaged free-energy density to derive an expression for the internal energy U = ∂(βf )/∂β of the disordered phase. Then, the specific heat c = ∂U/∂T can be evaluated numerically, or analytically at the high-and low-temperature limits (see Appendix A 3 for more details). At high temperatures, we find that where b = (Ω 0 /Ω u ) 2 − (Ω v /Ω u ) 5 /3 /2. The high-T limit, where c → 3/4, is a result of the quartic term which becomes dominant as T ≫ Ω u (see further discussion in A 3). At low temperatures the specific heat vanishes exponentially as the system is gapped at T = 0. In Fig. 3, we show the T -dependence of c for a representative set of parameters, which demonstrates that c → 3/4 at high temperatures, and vanishes as T → 0 (see inset). We will further discuss the temperature dependence of c and the correspondence between c and τ ph in Sec. IV.

IV. DYNAMICS
In this section, we discuss the dynamical properties of the disordered phase of the SB model. We identify three dynamical regimes that posses different spectral properties for which the phonon lifetime τ ph is a convenient probe. In particular, we demonstrate the existence of an intermediate-T 'Planckian' regime, where τ ph is of the order of the Planckian time scale τ Pl = /T . Focusing on the disordered phase, where the model is self-averaging, we may use the framework of the Keldysh formalism (see Appendix B) to compute the disorderedaveraged partition function, which can be expressed as a functional integral of an effective Keldysh action. The functional integral is controlled by the saddle point of the effective Keldysh action in the limit of N → ∞. This enables us to obtain a set of self-consistent equations for the well-known retarded, advanced and Keldysh Green's functions, from which one can extract the desired dynamical information. The Keldysh saddle-point equations are given bŷ where the Keldysh component is set according to the fluctuation dissipation theorem (FDT) to enforce thermal equilibrium.
For a wide range of temperatures, Eqs. 7 can be numerically solved by an iterative procedure (see Appendix C). Solving these equations enables us to obtain the spectral function and the phonon lifetime τ ph , defined through the late time behavior of the retarded propagator: G R (t) ∝ e −t/τ ph . Fig. 1a,b shows the dependence of τ ph , normalized by τ Pl = /T , on T for a system with Ω 0 ∼ Ω v ∼ Ω u . The parameters are chosen such that the system is always far from the glass phase. We identify three distinct temperature regimes. At low temperatures, T ≪ Ω 0 , we find that τ ph rises sharply with decreasing T ; this behavior is associated with the finite gap at T = 0. At high temperatures, T ≫ Ω 0 , τ ph becomes temperature independent, and hence T τ ph / ≫ 1. As we will argue below, this is the classical regime, where τ ph ∝ √ M and is independent of . T τ ph / has a shallow minimum at intermediate temperatures.
We refer to this intermediate regime as the 'Planckian regime', where τ ph ∼ /T . Note that within our model, this regime is not parametrically large, but rather extends over a finite range of temperatures around T ∼ Ω 0 . The minimal value of T τ ph / depends on the parameters Ω 0 /Ω v and Ω u /Ω v . As we will show below ( Fig. 4), for an appropriate choice of these parameters (in the vicinity of the glass phase), we find that T τ ph / approaches 1.
At the two limiting cases of high and low temperatures, Eqs. 7 are amenable to analytical approximations, from which one can extract the qualitative behavior of the phonon lifetime. To carry out these approximations, we assume that (a) the phonons are well-defined quasiparticles at these regimes, and that (b) there are only two energy scales in the system, corresponding to the renormalized phonon frequency and the phonon lifetime. We then use these assumptions to make the following ansatz for the retarded Green's function, where (a) means that we assume √ Π 0 ≫ τ −1 ph and we use (b) to ignore higher order terms in the retarded selfenergy. Note that the ansatz in Eqn. 8 contains two unknown quantities: τ ph and Π 0 . Note also that Π 0 is a thermodynamic quantity, that corresponds to the zerofrequency phonon stiffness. Details on the imaginaryand real-time derivations for Π 0 are found in Appendices A 1 a and B 1, respectively, and in Appendix B 2 we give details on the estimation of τ ph .
At low temperatures (T ≪ Ω 0 ), the system is essentially gapped. We then expect the phonons to be exponentially long-lived, since scattering off thermal excitations is exponentially rare. We confirm this expectation and find that At high temperatures (T ≫ Ω u ), we find that Π 0 ≈ √ uT and where r = Ω v /Ω u = v 2/5 /u 1/3 is a dimensionless number of order one (r ≤ 1 in our setting of interest). Interestingly, it appears that the phonon lifetime becomes independent of temperature, which agrees with the numerical data shown in Fig. 1b. Moreover, by reinstating and M in Eqn. 10 (using the definitions of Ω u and Ω v above Eqn. 2), we find that and is independent of . This suggests that for T ≫ Ω u , the model obeys classical dynamics (see also B 2 a). Specifically, the dynamics of underdamped harmonic oscillators, as To see that the aforementioned inequality holds, we use the fact that at high temperatures Π 0 = √ uT (see Eqn. A9). Then, substituting τ ph and Π 0 , we see that the inequality holds only if T ≫ r 10 Ω u . This is indeed the case, because r ≤ 1 and T ≫ Ω u .
Note that the parameters in Fig. 1a,b (with Ω 0 /Ω v = 1.1) are identical to the ones in Fig. 3. It is then interesting to examine the correspondence between the behavior of the specific heat and the phonon lifetime. At low temperatures, we find that c and 1/τ ph vanishes exponentially, as expected due to fact that the system is gapped at T = 0. At high temperatures, the system approaches the classical limit, as can be seen thermodynamically by the fact that c approaches a constant value, and dynamically as τ ph ∝ √ M . Importantly, at intermediate temperatures (T /Ω v ≈ 0.45 to T /Ω v ≈ 1.15), which we referred to as the Planckian regime in Fig. 1a, we find a significant variation in the value of c, which serves as another indication of the quantum mechanical nature of this dynamical regime.

Minimal phonon lifetime
We have thus far demonstrated that the model has an intermediate-temperature dynamical regime where the phonon lifetime is of the order of the Planckian time scale with a numerical coefficient α. It is interesting to ask what is the minimal attainble value of α within our model. For a generic choice of parameters in the stronglycoupled regime, where one finds that the numerical coefficient α in Eqn. 12 is of the order of 10 around the Planckian regime, as demonstrated in Fig. 1. However, as one approaches the vicinity of the glass phase in the (Ω 0 , T ) plane (for fixed u and v), this numerical coefficient tends to decrease. In particular, we find that for sufficiently large values of u, which enable us to approach relatively small values of Ω 0 and T and remain in the disordered phase, α reaches values close to unity, see Fig. 4. Scanning the (Ω 0 , T ) parameter space, we find that α 1, and that α ≈ 1 for regions with Ω 0 < Ω c (u, v) and temperatures slightly above the glass transition, where Ω c (u, v) is the T = 0 glass transition frequency. This is apparent in Fig. 4. The picture for other values of u is similar. Interestingly, we find that α never drops below 1, supporting the conjecture of a universal bound on α.

V. GENERALIZATION TO MULTIPLE PHONON BRANCHES
Considering the simple, SB version of the model enabled us to obtain its phase diagram and specific heat, and served as a convenient platform for the study of its real-time dynamics. This version, however, describes a rather artificial setting in terms of phonons, where we consider N degenerate optical phonon branches (with Ω i = Ω 0 for all i's). In physical insulating compounds, optical phonon branches are typically spread over a finite bandwidth, rather than being degenerate. As a step towards making our model more realistic, we consider a multi-phonon branch generalization, where the Ω i 's satisfy Ω 1 ≤ ... ≤ Ω N . We dub this version the 'multi-branch' (MB) model. The large-N limit is taken such that the distribution of frequencies obeys We consider the Green's function for the ith branch, The generalization of the imaginary-time self-consistent equations in the disordered (replica-diagonal) phase is given bŷ As in the SB case, we first need to determine the phase diagram of the model. The replica analysis for the MB model is more complicated than in the SB case. Instead of calculating the phase diagram explicitly, we use a simple argument to bound the glass phase in the (Ω min , T ) plane. Consider deforming the mode distribution function f (Ω) continuously to that of a SB model with f SB (Ω) = δ(Ω − Ω min ). Such a deformation softens the phonon modes, stabilizing configurations with large equilibrium displacements. We therefore expect that the deformation expands the regime of the glass phase. Indeed, in the SB model, decreasing Ω 0 brings us closer to the glass regime (Fig. 2). Therefore, we assume that for fixed Ω v,u and Ω min > Ω c (u, v) (where Ω c is the location of the T = 0 glass transition in the corresponding SB model), the MB model is in the disordered phase.
We now consider the real-time dynamics of disordered phase of the MB model. Here, the ith phonon branch is characterized by its corresponding spectral function A (i) (ω) and phonon lifetime τ rameters of the MB model 3 .
In Fig. 5a we show the dependence of τ (i) ph /τ Pl on T for 4 out of 10 distinct branches. Note that two of these branches are at the edges of the spectrum, Ω 1 and Ω 10 , and two are near its center, Ω 4 and Ω 7 . Here, the parameters are chosen such that the system is always far from the glass phase, by the assumption stated above. We find that the overall trends are similar to those found in the SB model (Fig. 1a). In particular, we identify the same three dynamical regimes: the low-T semiclassical regime where the lifetime diverges exponentially, the high-T classical regime where the τ In Fig. 5b,c we show the spectral function of the 10 distinct branches for two temperatures, where the significant broadening of the spectral peaks demonstrates the crossover from the semiclassical regime, with T /Ω v = 0.2, to the Planckian regime of the MB model, with T /Ω v = 1.

VI. RELATION TO THE SYK MODEL
It is natural to ask whether our model realize a critical point with an emergent conformal symmetry at low energies, similarly to the SYK model [28,29]. Unfortunately, the answer appears to be no, both in our model, and also in general bosonic variants of the SYK model, as was recently discussed in [34].
We begin by showing why the most naive approach to finding a conformally invariant point is inconsistent, and then we will show that the consistent approach leads to an unstable conformally invariant critical point, which is not the physical solution of the saddle-point equations. Naively, the first step towards an emegernt conformal (time-reparametrization) symmetry at low energies is to bring the self-consistent equations to a reparametrization-invariant form. Observing Eqs. 5 at T = 0, this can be done by neglecting the ω 2 term, while tuning Ω 0 and u to zero 4 . This leads to the following set of 'SYK-like' self-consistent equations: Then, following the analogy with the SYK model, we substitute a scaling ansatz, given by G conf (τ ) ≡ b |τ | −2∆ φ , into the self-consistent equations, and obtain the scaling dimension ∆ φ and the numerical coefficient b. However, the solution obtained by this approach is an inconsistent solution of Eqs. 15. The inconsistency comes from the UV behavior of the self-energy. Namely, the UV piece of the self-energy is non-negligible: . That is, the leading term of Π(iω → 0) is a constant, Π UV , rather than the conformal self-energyΠ conf (iω) with which we started our consistency check. In addition, due to its positive sign, this constant violates the positivity condition thatĜ(iω) > 0, sinceĜ(iω = 0) = −1/Π UV < 0. This behavior originates from the bosonic nature of the degrees of freedom, which manifests itself in the even parity of the Green's function. In the fermionic case, the odd parity of the Green's function enables one to neglect the UV piece of the Green's function safely.
This inconsistency can be fixed by reinstating Ω 0 , and using it as a counter-term, such that the self-energy readŝ We can then fine tune Ω 0 to cancel the non-conformal constant by letting The scaling ansatz is then a consistent solution of the saddle-point equations at low temperatures and long times, which can be verified both analytically and numerically (see A 1 b and C 1 a, and also Section 2 of [35]). Interestingly, there are several warning signs indicating that this scale-invariant solution is not realized in our model. The most obvious one in our setting is the fact that for Ω 0 = Ω conf and u → 0, the system realizes the glass phase associated with the 1SRSB solution of the saddle-point equations. However, even within the disordered phase, it turns out that the conformal solution is unstable.
In Ref. [36], similar conformally-invariant selfconsistent equations were studied in the context of bosonic tensor models (without disorder), with q-body rather than 3-body interactions (q ≥ 4). Notably, [36] found that the scaling dimension of the φ 2 composite operator is complex for the conformal field theory (CFT) associated with this form of self-consistent equations. This violates the unitarity condition and implies that the CFT is an unstable solution of the theory. This is also true in our model [37], and it suggests that there exists another, stable solution of the disordered saddle-point equations with Ω 0 = Ω conf . Indeed, we find that such a solution exists. This solution is gapped at T = 0, and by comparing the free energies of the two solutions, we also find that it is thermodynamically favorable (see further discussion below Eqn. C3 in Appendix C). A similar observation was made in [38].

VII. DISCUSSION AND OUTLOOK
In this work, we have studied a solvable model of N interacting phonons. In the limit N → ∞, the model is solvable for any interaction strength and temperature, allowing us to access the strongly interacting regime. At low temperature and strong interactions, the system undergoes a first order transition into a replica symmetry breaking (glass) phase. Focusing on the dynamics in the replica diagonal (disordered) phase, we find that the system crosses over between three distinct regimes as the temperature increases: a semiclassical regime with longlived quasiparticle (phonon) excitations at low temperatures, a classical regime at high temperatures, and an intermediate strongly-interacting "phonon fluid" regime. In the latter regime, the minimal phonon relaxation time is of the order of the Planckian time scale, τ ph = α /T , with α approaching unity near the transition to the glass phase.
Our work was motivated by measurements of the thermal diffusivity in a broad class of insulating materials, indicating that these systems may indeed be described as a strongly coupled liquid of phonons, with a relaxation time that approaches the Planckian time. Clearly, our model is not meant to realistically model any material; rather, it provides a concrete example of such a strongly interacting quantum regime in a bosonic system. Within our model, this regime is realized over an intermediate temperature range; at sufficiently high temperatures the system always crosses over to a classical regime, at which τ ph ≫ /T and the specific heat approaches its classical limit. Interestingly, we find that in the intermediate quantum regime, the relaxation time never drops below /T , consistent with the notion of a universal "Planckian bound" on thermalization times.
From a theoretical perspective, our model can be viewed as a bosonic variant of the SYK model. However, there are crucial differences between our model and the fermionic SYK model. In particular, we showed that the emergent low-energy conformally invariant saddle-point solution of this model is not realized, as it is found deep inside the glass phase. This is in line with general arguments regarding the low-temperature behavior of bosonic SYK-like models [34] and the presence of operators with complex scaling dimensions at the putative conformally invariant point [36].
Some natural questions remain open. To address the transport properties in the quantum phonon fluid regime, our model needs to be generalized to higher dimensions. This can be done, e.g., by placing a copy of our model on each site of a D−dimensional lattice, along the lines of Ref. [39]. Related to this issue is the absence of acoustic phonon modes in our model. These are protected by Goldstone's theorem, and must remain gapless and longlived even in the presence of strong interactions. Nevertheless, in the strongly coupled "phonon fluid" regime, their contribution to transport may be negligible due to their small phase space.
Another natural question regards the effect of glassiness on the dynamics in our model. In the glass phase, ergodicity is violated and the phase space is fragmented into disconnected clusters. However, intuitively, one may expect the relaxation dynamics within each phase space cluster to be qualitatively similar to that of the disordered phase. We leave a detailed investigation of this question to future studies. S. Sachdev, T. Senthil for useful discussions throughout this work. EB was supported by the European Research Council (ERC) under grant HQMAT (grant no. 817799) and by the US-Israel Binational Science Foundation (BSF).
Here we give a few more details on the replica analysis. Starting from Eqn. 3, the replicated partition function is given by where φ = {φ α i : i = 1, ..., N ; α = 1, ..., n}, the disorder measure is and the action of each replica α = 1, ..., n is To proceed, we integrate over the disorder and introduce composite fields G αβ (τ, τ ′ ) and the and Lagrange multiplier fields enforcing these constraints, mentioned below Eqn. 3. The integration is straightforward and the implementation of the Lagrange multipliers is done with the identity [40] f This identity is enforced for all τ, τ ′ that satisfy (τ, τ ′ ) ∈ [0, β] 2 such that τ > τ ′ . This enforcement avoids the redundancy coming from the even parity of the imaginary-time Green's function.
Then, we integrate out the phonon fields {φ α i } and obtain that where DA = αβ DA αβ with A = G, Π, and the effective action is given by Eqn. 4.

Replica-diagonal solution
The diagonal solution is given in and above Eqn. 5. Here we discuss some other aspects related to the diagonal solution.

a. Instability of the disordered phase
We comment on an instability that arises at intermediate temperatures, where the zero-frequency phonon stiff-nessĜ −1 (iω n = 0) ≡ Π 0 becomes negative. Given that the model is well-defined, in the sense that the energy is bounded from below, the existence of such instability is a first indication that the replica-diagonal solution is insufficient in this region of parameter space, suggesting the existence of a glass phase, as is later confirmed.
At high temperatures, the quartic term is dominant and the zero-frequency phonon stiffness is simply given by inverse of the positive renormalized phonon frequency. Problems may arise at low-to-intermediate temperatures, where the cubic term might become dominant.
Let us focus on intermediate temperatures. We later show that this instability is can avoided at low temperatures as long as Ω 0 is sufficiently large. To proceed, we solve the saddle-point Eqs. 5 for Π 0 . We assume that T is sufficiently large such that the Matsubara summation may be approximated by the ω n = 0 component: and substituting in Π 0 reads The instability is characterized by parameters Ω 0 , u, v and T for which the only real solution of Eqn. A7 is negative. We denote the ratio between the cubic and quartic energy scales as To demonstrate the existence of an instability, we examine the cases where r → 0 and r ≫ 1. We first consider the r → 0 case, which corresponds to setting v → 0 and u > 0. Since the origin of this instability is the cubic term, this case is expected to show no instabilities. Indeed, we find that the real solution to Eqn. A7 for T ≫ u 1/3 is given by where we have assumed for simplicity that Ω 0 ∼ Ω u . Note that this form of Π 0 also holds for T ≫ Ω u in the strongly coupled regime Ω 0 ∼ Ω v ∼ Ω u . Now consider the case of r ≫ 1, which corresponds to setting u → 0 and v > 0. Here we do expect an instability due to the fact that the cubic term becomes dominant at sufficiently high temperatures. Indeed, the only real solution of Eqn. A7 at T ≫ v 2/5 is given by where we have assumed for simplicity that Ω 0 ∼ v 2/5 . The requirement Π 0 > 0 is clearly violated, indicating that the replica-diagonal solution is unstable for u → 0.
For large values of r, the instability exists at intermediate temperatures v 2/5 T rv 2/5 . In general, the instable region in parameter space for which Π 0 < 0 is also a function of Ω 0 , and can be characterized by solving the cubic polynomial in Eqn. A7 and demanding that the solution that is connected to the solution in Eqn. A9 is positive for all T Ω v . However, the full characterization is not needed as the system undergoes a phase transition before it encounters this instability. This can be seen in Fig. 2, where we find that smaller values of u (larger values of r) correspond to a larger regions that realizes the glass phase.
At low temperatures, by approximating the Green's function asĜ(iω n ) −1 ≈ ω 2 n + Π 0 , we find that which gives the following equation for Π 0 : One can check that the limit u → 0 still allows for positive solutions for Π 0 as long as Ω 0 v 2/5 , whereas for intermediate temperatures in the limit of u → 0, the instability exists even for relatively large values of Ω 0 . Interestingly, for generic values of Ω 0 , v and u, one may find multiple real and positive solutions for Π 0 , from which one can understand the existence of multiple saddle points for the same set of parameters as mentioned briefly in A 1 b.

b. Scale-invariant ansatz and relation to SYK
We show that the naive scaling ansatz is not a solution of the saddle-point equations. Recall that we have denoted the scaling ansatz by G conf (τ ) ≡ b/|τ | 2∆ φ . By substituting the scaling ansatz in Eqs. 15, one can check that b −3 = 3v 2 |Γ (−1/3)| Γ (1/3) and ∆ φ = 1/3. We assume that the scaling ansatz is valid up to a short-times cutoff Λ −1 . We then decompose the full two-point function as (A13) where Θ(x) = 1 if x > 0 and zero otherwise, and the short-times piece is given by G UV (τ ) =´| ω|>Λ dω 2π In the fermionic case, the contribution of this short-times piece to the self-energy may be neglected due the odd parity of the fermionic Green's function. Here this is not the case. Instead, it contribute as a constant, which implies that our naive approach, where we took Ω 0 → 0 to obtain a set of 'SYK-like' saddle-point equations (Eqs. 15), is inconsistent.
For completeness let us note that the mapping from T = 0 to T > 0 is identical to the fermionic SYK case (see e.g. [29]) with the scaling dimension ∆ φ = 1/3. At finite T ,Ĝ We have shown that the model may be fine-tuned to a conformally-invariant critical point by a deformation of the bare phonon frequency Ω 0 → Ω conf . As noted in Sec. VI, however, the disordered saddle-point equations admits a different, gapped solution for this set of parameters (u → 0, Ω 0 → Ω conf ) which is thermodynamically favorable.

1-step replica symmetry breaking solution
As shown in [30], the glass phase of the model is described by a 1-step replica symmetry breaking (1SRSB) solution at the level of the saddle-point approximation 5 . The following derivation is largely along the lines of [30]. We begin from the 1SRSB solution, where ǫ αβ = 1 if α and β are in a diagonal block of size m and ǫ αβ = 0 otherwise. Interestingly, the off-diagonal terms of G αβ (τ ) can be shown to be τ -independent [42].
As in [30], the absence of a linear-in-φ i term in H implies that g 0 = 0. To proceed, we preform the following steps. We substitute the the self-energy, obtained by the variation of S eff with respect to Π αβ (τ ), back in S eff . Then we move to Matsubara space. And lastly we substitute the (replica space) eigenvalues and corresponding degeneracies ofĜ αβ (iω n ), which are given in terms of g d (iω n ) ,ĝ EA and m, in the ln det term. These steps will allow us to easily take the limit n → 0, where m is then analytically continued to take real values between zero and one such that m = 1 corresponds to the replica diagonal solution and m = 0 to the replica symmetric solution. We continue by scaling dimensionful quantities with respect to Ω v = v 2/5 , and to lighten the notation, we leave the dimensionless parameters with the same notation. So in practice, we simply set v = 1 wherever it appears and remember that all other parameters are scaled with respect to the appropriate power of v. The saddle-point equations are then given by where the self-energy is given by Π (τ ) = g d (τ ) 2 − ug d (τ ) δ(τ ). Here Eqs. A19,A20 and A21 are obtained by varying S eff with respect toĝ d (iω n ),ĝ EA and m, respectively. Note that we have already eliminated solutions for which g EA = 0 or m = 1 in the saddle-point equations.
Following [30], we define two parameters, y ≡ βg EA /ĝ d (0) and x ≡ my/ (1 − y). Substituting these in Eqn. A20 and Eqn. A21 gives an equation for x, Numerically solving Eqn. A22 gives x = 1.81696. It will be useful to notice that In particular, note thatĝ d (0) is fixed by m and β. We then separate g d and the self-energy to a constant and τ -dependent parts, Substituting Eqn. A24 and Eqn. A25 into Eqn. A19, with the help of the relations in Eqn. A23, the terms proportional to δ ωn,0 cancel and we getĜ (iω n ) −1 = ω 2 n + Ω 2 0 −Π (iω n ). It is important to remember to respect the saddle-point constraint onĝ d (0), which impliesĜ (0) = mβg EA /x. Therefore, for a given m and β, Ω 0 must be chosen selfconsistently such that the constraint is satisfied. That is, we are forced to set Finally, we arrive at a closed set of equations forG and Π,Ĝ (iω n ) = 1 These equations can be solved numerically (see Appendix C) to obtain the free-energy density of the glass phase and construct the phase diagram of the model. In practice we fix m and β, numerically solve forĜ,Π and then extract Ω 0 from Eqn. A26.

Thermodynamic functions
The thermodynamic functions (free energy, internal energy and specific heat) of the model are obtained from the effective action in Eqn. 4. The free energy of model is obtained by substituting the definitions of the replica space solutions into Eqn. 4. Then, we find that the free energy of the 1SRSB solution is given by and for the diagonal solution, which corresponds to taking the limit m → 1 or the limit g EA → 0, we have Here, C is related to the regularization of the ln det term, C = n ln β 2 ω 2 n + Ω 2 0 . It can be shown that C = 2 ln 2 sinh βΩ0 2 [43] (which is simply the free energy of an harmonic oschillator with frequency Ω 0 ). To derive the internal energy of the disordered phase, however, one can use the unregularized form of the free energy, as we show next.
The internal energy U is defined by ∂(βf )/∂β. Due to the saddle-point approximation, the differentiation is non-vanishing only when it acts explicitly on β: To proceed, we substitute the saddle-point Eqs. 5 in Eqn. A29 and obtain that The derivative of the first term is evaluated as where ω n ≡ 2πn. To differentiate the first summand in the second line, we may rescale the intergration variable τ → x ≡ τ /β, similarly to [30]. Finally, we obtain that (A30) Let us use Eqn. A30 to evaluate the specific heat c ≡ ∂U/∂T . We will show the high-temperature specific heat. At low temperatures the evaluation is similar and one finds that the c vanishes exponentially as T → 0 since the system is gapped at T = 0.
For T ≫ Ω u , we use Eqn. A9 to approximate the Green's function as G(τ ) ≈ T /u. Substituting this approximation in Eqn. A30, we obtain the hightemperature specific heat, given in Eqn. 6. The sign of prefactor b determines whether c has a maximum at an intermediate temperature (b > 0) or saturates to its maximum at T /Ω u → ∞ (b < 0).
Moreover, note that the high-temperature limit lim T /Ωu→∞ is not an artifact of the saddle-point approximation. It can be derived by a simple scaling argument: Consider a classical system (with position X and momentum P ) where the highest order term in the potential is quartic, and denote it by uX 4 . At sufficielty high temperature, the quartic term is dominant. The partition function of the system may thus be approximation as By rescaling Hence the free energy is given by βF = ln(Z) = 3 4 ln(T ) + ln(terms independent of T ) (A33) and it follows that c → 3/4 as T → ∞.

Appendix B: Keldysh formalism
We derive the Keldysh action and saddle-point equations of the disordered phase to study the dynamical properties of the phonons, and in particular the phonon lifetime. The derivation below is done in the spirit of [44]. In this formalism we calculate the partition function Z = Tr e −βH U /Tr e −βH where U is the identity real-time evolution operator, evolving forward from time t 0 to t f and backward from t f to t 0 . As usual, the label + (−) denotes the Keldysh forward (backward) contour.
The disorder-averaged Keldysh partition function Z is given by where φ = {φ i : i = 1, ..., N }, the disorder measure was defined in Eqn. A2, and the Keldysh action is Similarly to the replica method, we average over the disorder and introduce composite fields and Lagrange multiplier fields enforcing these constraints Π ss ′ (t, t ′ ). We then integrate out the phonon fields to obtain the disorder-averaged partition function in terms of G and Π, where the effective Keldysh action is given by The saddle-point equation are then obtained by vary-ing the effective Keldysh action with respect to G ss ′ and where we have assumed real-time translation invariance.
To obtain the saddle-point equations in terms of the more conventional retarded, advanced and Keldysh Green's functions, we introduce a Keldysh rotation matrix, which satisfies and The saddle-point Eqs. B6 do not contain the information about the initial thermal equilibrium density matrix.
The generalization to the MB model is given bŷ G aK (ω) = 2i coth βω 2 Im Ĝ aR (ω) , where a, b = 1, ..., N . Note that this approach is equivalent to the Keldysh diagrammatic approach when considering the O(1) diagrams in a 1/N expansion, see e.g. [45].
1. Real-time derivation of Π0 at high T In Section A 1 a, we derived an equation for Π 0 in imaginary-time, where the high-temperature limit allowed us to approximate the Matsubara summation by the zeroth Matsubara frequency. Here we show that Eqn. A7 can be derived exactly by solving the Keldysh equations forΠ R (0).
To begin, we use Eqs. 7 and substitute the Keldysh Green's function inΠ R (0), where we have used the even parity ofĜ K .
To proceed, we expand the coth to leading order, thereby replacing the quantum with its classical version. We later show that this is consistent for T ≫ Ω u . We haveΠ Using the Kramers-Kronig relations, the second summand in Eqn. B12 is given by In the following we show that Then, by substituting the two summands back to Eqn. B12, together with the fact that Π 0 =Ĝ R (0) −1 , we obtain Eqn. A7 from the Keldysh equations. To derive Eqn. B15, we rewrite I as Using the fact that dω 2π we can write I as Expanding the square brackets, the first term iŝ This is because G R (t) = 0 for t < 0, and hence the sign functin can be replaced by unity. On the other hand, using G A (t) = G R (−t), the second term vanishes, sincê where we used the anti-symmetry of the integrand under the exchange of t 1,2 . Hence, we have derived Eqn. B15. We comment on the validity of the substitution of the quantum FDT by its classical version at the end of Section B 2 a.

τ ph at low and high temperatures
Here, we provide a simple derivation for the phonon lifetime in the limit of low and high temperatures in the disordered phase of the SB model. We assume that the system's parameters are chosen such that we are always far from the glass phase. To obtain an equation for the phonon lifetime, we consider the imaginary part of the retarded self-energy. SubstitutingĜ K intoΠ R in Eqn. 7 and taking the imaginary part of both sides reads (B21) Here A(ω) ≡ ImĜ R (ω) is the spectral function. To proceed, we use the ansatz for the retarded propagator, given in Eqn. 8. It is then sufficient to consider the ω → 0 limit of Eqn. B21 in order to extract τ ph . The zeroth order in ω vanishes, and the leading term is given by (B22) We integrate the RHS by parts, use the fact that the spectral function decays at |ω| → ∞ and use the parity of the integrand to obtain In the following, we will substitute the spectral function obtained from our ansatz (Eqn. 8) and solve for τ ph at the high-and low-temperature limits.
a. High T High temperatures can be identified with the classical limit of the model in the following manner. In general, classical systems that obey Newtonian dynamics are invariant under the rescaling of the real-time coordinate t →t = t/ √ M , where M is the mass of the classical degree of freedom. This implies that any time scale, and in particular τ ph , should be proportional to √ M at the classical limit. Furthermore, as this is a property of the equation of motion, we should be able to see this rescaling invariance at the level of the Keldysh saddle-point equations. It is easy to check that the Keldysh equations are invariant under this rescaling in frequency space (where ω →ω = √ M ω) if the temperature satisfies where ω sf is the characteristic frequency beyond which A(ω) becomes negligible. In this limit, the rescaling invariance follows from the fact that the quantum FDT can be replaced with its classical version. In practice, Eqn. B24 holds when T / is much larger than the largest frequency scale in the system. In terms of our ansatz, the classical limit of the model is identified with T ≫ Π 0 (T ). Recall that in Sec. B 1 and Sec. A 1 a we have derived an equation for Π 0 at high temperatures, whose solution is given by Π 0 (T ) ≈ (uT ) 1/2 . Using this form, we may identify the classical limit with temperatures that satisfy T ≫ Ω u . Note also that the criterion in Eqn. B24 is highly incompatible with the Planckian regime, where ω sf ∼ T .
We now proceed to estimate τ ph . In terms of our ansatz, the spectral function reads Substituting Eqn. B25 in Eqn. B23, and using T ≫ ω sf ∼ Π 0 (T ) to expand the csch function to leading order, we have that where in the last line we have used the assumption that τ ph √ Π 0 ≫ 1 to ignore subleading contributions in 1/τ ph √ Π 0 . Finally, we substitute Π 0 (T ) ≈ (uT ) 1/2 into Eqn. B26 and obtain the high-T phonon lifetime as given by Eqn. 10.
As for the consistency of this approximation and the one in section B 1: we have found that Π 0 ≈ (uT ) 1/2 for T ≫ Ω u , which then implies that (see discussion below Eqn. 11) Namely, the spectral function at high-T is peaked around ω ≈ √ Π 0 and its width is much smaller than √ Π 0 In this case, one can indeed replace the quantum FDT by its classical version as the integrands in Eqn.B11 and Eqn.B26 are essentially supported in a frequency interval for which ω/T ≪ 1, as we have assumed in our selfconsistency argument.
b. Low T At low temperatures we expect the phonon lifetime to be exponentially long, τ ph ∼ e √ Π0(T →0)/ηT , where Π 0 (T → 0) is the T = 0 gap and η is some numerical coefficient.
Now, in the limit T ≪ ω sf ∼ Ω 0 , which corresponds to ω ′ β ≫ 1 in the support of the integrand in Eqn. B23, we may approximate csch(ω ′ β/2) ≈ 2e −ω ′ β/2 . Then, Eqn. B23 is given by We proceed by using the ansatz Eqn. B25, (B29) Note that if τ ph √ Π 0 , τ ph T ≫ 1, the leading contribution to the integrand is coming from ω ′ ≈ √ Π 0 . We may therefore approximate the integral similarly to the high temperature limit and obtain Finally, the phonon lifetime is given by Eqn. 9.
Appendix C: Numerical solution of the saddle-point equations Here we give some details on the numerical solution of the self-consistent equations in imaginary-and real-time for the single-branch model. The generalization of these methods to the MB model is straightforward with the use of the summation described in footnote 3. We solve the saddle-point Eqs. 5 iteratively following the method of [11], with small modifications to be specified ahead. We first describe the unmodified algorithm. At the zeroth iteration step we use an initial con-ditionĜ −1 0 (iω n ). The Matsubara frequencies are given by ω n = 2πβn where −P e ≤ n ≤ P e such that the total number of sampling points is P = 2P e + 1. After the jth iteration step, we obtainĜ j , and substitute it intô Π j+1 (iω n ) = v 2 β ω ′ nĜ j (iω n ′ )Ĝ j (iω n − iω n ′ ) − uG j (0) (C1) which we implement using MATLAB's conv_fft2 function package. The convolution function outputs a length 2P − 1 vector from which we take only the P components that contain frequencies in the same window asĜ j . We then update the two-point function aŝ G j+1 (iω n ) = (1−X)Ĝ j (iω n )+X 1 where X ∈ (0, 1).
After each iteration step we monitor the error e j ≡ |G j − G j+1 | 2 and if e j > e j−1 we update X → X/2, as long as X is larger then a minimal updating factor (to ensure convergence). The iteration procedure is terminated when |e j − e j−1 | < ǫ. Usually we start with X = 1/2, set the minimal updating factor to X min = 1/100, and set ǫ = 10 −14 . Within the unmodified algorithm that we are currently discussing, we obtain the Green's function for a given Ω 0 adiabatically by approaching it from above. That is, we start by solving for the Green's function with a large Ω 0 (e.g. Ω 0 /Ω v = 4), with the initial condition of a free phonon. Then we slightly decrease Ω 0 and use the previously obtained solution as the initial condition for the slightly decreased Ω 0 . In general, one may encounter multiple solutions for saddle-point equations as a function of Ω 0 . However, the solutions obtained with the above method are associated with the thermodynamically favorable saddle point of the disordered phase, as we have checked by comparing the free energies of the different solutions.
We proceed to describe the modified algorithm. This modification enables us to approach thermodynamically unfavorable solutions of the disordered saddle-point equations, and conveniently solve the saddle-point equations in the glass phase. Instead of treating Ω 0 as an input (together with v, u and β), we fixĜ(iω = 0) and then ask what Ω 0 corresponds to suchĜ(iω = 0). To implement this we introduce δΠ ≡Ĝ(0) −1 and redefine the iteration step as follows, G j+1 (iω n ) = (1−X)Ĝ j (iω n )+X 1 whereΠ j+1 (iω n ) ≡Π j+1 (iω n )−Π j+1 (0). Once the iteration process converges we can extract its corresponding Ω 2 0 , given byΠ(0) + δΠ. Note that the correspondence is not single valued due to the existence of multiple saddlepoint. The solutions with this approach are always obtained with a free initial condition. No adiabatic tuning is needed.