On Phases Of Melonic Quantum Mechanics

We explore in detail the properties of two melonic quantum mechanical theories which can be formulated either as fermionic matrix quantum mechanics in the new large $D$ limit, or as disordered models. Both models have a mass parameter $m$ and the transition from the perturbative large $m$ region to the strongly coupled"black-hole"small $m$ region is associated with several interesting phenomena. One model, with ${\rm U}(n)^2$ symmetry and equivalent to complex SYK, has a line of first-order phase transitions terminating, for a strictly positive temperature, at a critical point having non-trivial, non-mean-field critical exponents for standard thermodynamical quantities. Quasi-normal frequencies, as well as Lyapunov exponents associated with out-of-time-ordered four-point functions, are also singular at the critical point, leading to interesting new critical exponents. The other model, with reduced ${\rm U}(n)$ symmetry, has a quantum critical point at strictly zero temperature and positive critical mass $m_*$. For $0<m<m_*$, it flows to a new gapless IR fixed point, for which the standard scale invariance is spontaneously broken by the appearance of distinct scaling dimensions $\Delta_+$ and $\Delta_-$ for the Euclidean two-point function when $t\rightarrow +\infty$ and $t\rightarrow -\infty$ respectively. We provide several detailed and pedagogical derivations, including rigorous proofs or simplified arguments for some results that were already known in the literature.

Melonic theories are models for which the perturbative expansion is dominated, in some well-defined limit, by a special class of Feynman graphs called melonic. These graphs are relevant in a surprisingly wide range of problems. They were discovered for tensor models in the context of the discretized approach to quantum gravity [1] but also appeared implicitly in the study of disordered condensed matter systems [2]. They dominate a suitably defined large d limit of planar graphs [3][4][5] and thus provide an interesting non-perturbative approximation to the large n limit of matrix models. Moreover, melonic graphs can be recursively constructed and explicitly summed over using the technique of Schwinger-Dyson equations, allowing for considerable analytical progress to be made in the study of such models.
Strong interest in melonic theories, mostly fermionic quantum models, has been rekindled in recent times by the Sachdev-Ye-Kitaev (SYK) model [2,6] and various generalizations thereof [7,8], see also [9] and the review [10]. Within a program that has seen the convergence of very diverse approaches, new and interesting relations have been found between condensed-matter disordered systems [11][12][13][14], tensor models [15][16][17] and random matrix theory [18,19]. The emerging general picture is that melonic theories provide a quantum mechanical description of the near-horizon, low energy limit of near-extremal black holes [6,9]. The AdS 2 factor in the near-horizon geometry is at the origin of universal properties that are also found in the infrared limit of the melonic models [20][21][22]. A crucial ingredient is the emergence of an approximate time reparameterization symmetry. This symmetry is spontaneously broken down to SL(2, R), and the associated pseudo-Goldstone modes are governed by a Schwarzian action which matches the gravitational dynamics in nearly AdS 2 spacetime [23,24]. The low-temperature thermodynamical properties of the gravitational system are reproduced, together with some correlation functions and Lyapunov chaos exponents [25].
Melonic quantum models are also extremely interesting in their own right [1,16,[26][27][28]. The class of melonic graphs is situated in between the very simple trees of bubbles which dominate the large n limit of vector models, and the intractable set of planar graphs which are relevant for the matrix models. This "midpoint" has striking properties: the models remain simple enough to be amenable to analytical study; yet they display non-perturbative features typically associated with matrix models, including the emergence of a continuous spectrum and chaos. It is even plausible that the melonic large d approximation to matrix models [3] is able to capture the correct qualitative physics associated with the full sum over planar diagrams.
The aim of this paper is to present an extensive study of two fermionic quantum mechanical models in the melonic limit, motivated by the line of research initiated in [29]. The models can be formulated either as matrix quantum mechanics, or as disordered models. The random Hamiltonians in the disordered formulation read and where the χ i are complex fermions for 1 ≤ i ≤ N , λ ij kl and ξ i jkl are random couplings and m is a mass parameter. The model based on H is often called the "complex SYK model" and has already been considered in the literature, see e.g. [20,21,[30][31][32]. It preserves the fermion number symmetry and its matrix model version, discussed in Sec. 2, has a full U(n) × U(n) × O(d) symmetry. For this reason, we call this model the "U(n) 2 model." The model based on H has to the best of our knowledge never been studied before. It breaks the fermion number symmetry and its matrix model version has a reduced U(n) × O(d) symmetry; we therefore call it the "U(n) model." The Hamiltonians (1.1) and (1.2) look very similar, but this is misleading. We shall see that they yield qualitatively very different physics.
In the natural units for which the quartic coupling strengths are set to one, the models depend on two parameters, the mass m and the temperature T . Both the large T , fixed m and large m, fixed T regions of the (m, T )-plane are weakly coupled. The matrix model formulation makes it clear that the large mass regime describes a situation akin to a system of weakly interacting and well-separated D-particles governed by open string perturbation theory. At small masses, on the other hand, the models are strongly coupled and share properties with black holes. The introduction of the mass parameter thus allows us to study in detail the interesting transition between a weakly coupled "D-brane" regime and a strongly coupled "black-hole" regime [33]. We are going to discuss at length both the physics of Euclidean time (twopoint and four-point functions, thermodynamics and phase diagrams) and real time (two-point and out-of-time-ordered four-point functions, quasi-normal frequencies, chaos exponents).

Salient features and main new results for the U(n) 2 model
A fundamental property of the U(n) 2 model, first observed in [29], is the presence of a line of first order phase transitions in the strongly coupled region of the (m, T )plane, see Fig. 21 in Sec. 4.2. The first order phase transition exists because the Schwinger-Dyson equations resumming the melonic diagrams of the theory have two distinct solutions within a finite region of the (m, T )-plane. One solution is SYK-like and describes black-hole physics, whereas the other is an analytic continuation of the perturbative solution at large mass. The first order phase transition occurs when the free energies of the two solutions match. A crucial feature is that the line of first order phase transitions terminates at a critical point (m = m c , T = T c ). Standard thermodynamical critical exponents were evaluated numerically in [29] and found to be non-mean-field and asymmetric; for instance, at m = m c , different exponents are obtained depending on whether T → T c from above or from below.
In the present paper, we provide full derivations of the results presented in [29] and extend the analysis in two directions.
i) First, we study the phase diagram of the q-generalizations of the model (1.1), for which the quartic interaction is replaced by a fermion number preserving term of even degree q. The physics is similar for all values of q ≥ 4, with a line of first order phase transitions terminating at a critical point (m c ). We find that the values of some of the critical exponents do depend on q, see Table 2 in Sec. 4.2.3. The variations with q seem to be well above numerical error bars, especially for low values of q, providing evidence that the critical points found for different values of q are inequivalent. Moreover, we can also define "multi-q" models for which interactions of various degrees are included simultaneoulsy. We have initiated the study of such models, whose phase diagrams depend on the relative strengths of the couplings, see Fig. 31 in Sec. 4.2.4 for the case where the quartic and sextic couplings are both turned on. We find in this case evidence for a line of continuously varying non-trivial critical points, interpolating between the pure q = 4 and q = 6 cases.
ii) Second, we study real time correlation functions, from which we extract the leading quasi-normal frequency ν qn governing the exponential decay of two-point functions, and the Lyapunov chaos exponent η L governing the exponential growth of the connected piece of out-of-time-ordered four-point functions. We find that the chaos bound of [34] is saturated in the zero-temperature limit for all values of the mass in the black-hole, SYK-like phase. This is true for both the particular four-point function studied in [34], but also for the physical OTOC, for which the proof of the bound presented in [34] does not a apply. An interesting new result comes about when one studies the behaviour of ν qn and η L near the critical point (m c , T c ), see Fig. 44 and 50 in Sec. 5.1.4 and 5.2.6, respectively. We find that both ν qn and η L have a singular behaviour at the critical point. This singular behaviour leads to a new set of critical exponents, tied to the real time physics of dissipation and chaos and characterizing the critical point.

Main results for the U(n) model
In contrast with the situation in the U(n) 2 model, for any strictly positive temperature the U(n) model undergoes a smooth crossover from the weakly coupled large mass regime to the strongly coupled small mass regime. In particular, we find only one physically acceptable solution of the Schwinger-Dyson equations for each set of values of the parameters m and T . Interesting physics thus occurs at zero temperature. The mass gap m eff is found to decrease from its perturbative value m eff m at large mass to zero for a strictly positive critical mass m = m * > 0, see Fig. 33 in Sec. 4.3. Below the quantum critical point at m * the model is in an interesting new gapless phase.
The most novel and important feature of the model has to do with the infrared behaviour of this gapless phase found for 0 < m < m * . Instead of the standard conformally invariant solution, we find that the conformal symmetry SL(2, R), including the usual scale transformations, is spontaneously broken. This occurs because the power-law decay of the zero-temperature Euclidean two-point function is governed by two distinct exponents ∆ + and ∆ − , This non-standard behaviour, with the strict inequality (1.4), is established numerically in Sec. 3.3.2 with a very high degree of confidence.
The low-temperature thermodynamical properties of the exotic gapless phase are similar to the standard behaviour also found in SYK or in the U(n) 2 model. The entropy can be expanded as S/N = σ + cT + o(T ), with the zero-temperature entropy σ being strictly positive for all 0 ≤ m < m * and vanishing at m = m * . The coefficient c varies smoothly in the interval 0 ≤ m < m * and jumps discontinuously to zero at m = m * , see Fig. 35 in Sec. 4.3. On the other hand, we find convincing numerical evidence supporting the fact that the Lyapunov exponent does not saturate the bound η L ≤ 2πT when T → 0 (compare Fig. 47 for the U(n) 2 model with Fig. 48 for the U(n) model in Sec. 5.2.6). Note that there is no inconsistency here. The standard methods [9] to compute η L in the zero-temperature limit rely on the SL(2, R) symmetry; it is a non-trivial open problem to generalize these methods to a phase satisfying (1.3), (1.4) in which the SL(2, R) symmetry is spontaneously broken.

Plan of the paper
In Sec. 2, we introduce and motivate the U(n) 2 and U(n) models in detail. We explain that, to leading melonic order, they can be formulated equivalently either as large d planar matrix quantum mechanics with fundamental degrees of freedom given by d complex fermionic matrices ψ a µ b , 1 ≤ a, b ≤ n, 1 ≤ µ ≤ d (Sec. 2.1), or as large N disordered models with fundamental degrees of freedom χ i , 1 ≤ i ≤ N (Sec. 2.2). We emphasize the matrix formulation because it provides a genuine quantum mechanical theory even when the number of degrees of freedom is finite. We explain that these models are the most general single-trace matrix quantum mechanics based on the variables ψ a µ b , with quartic interaction terms and the additional condition that the fermion number violating two-point functions tr ψ µ ψ µ vanish at leading melonic order. We also comment in Sec. 2.3 on the relation with tensor models and the different large N limits that can be defined in these models; on subtleties associated with the case where the symmetry U(n) 2 × O(d) is broken down to U(n) × O(d), as in the U(n) model; and on gauging.
In Sec. 3 we introduce the Euclidean time technology which is heavily used in the rest of the paper. In Sec. 3.1 we recall, in a detailed and pedagogical way, the fundamental properties of Euclidean two-point functions, the Schwinger-Dyson equations they satisfy, some important properties of the space of solutions of these equations (non-trivial monodromies), their reparameterization invariance and the standard low energy ansatz. We also introduce and briefly discuss the fundamental new low energy ansatz (1.3) in Sec. 3.1.5. In Sec. 3.2 we derive explicit formulas for all the thermodynamical potentials (free energy, energy and entropy). These derivations are well known in the case of the disordered models, for which an effective action can be derived through a straightforward path integral argument based on the introduction of a suitable bilocal auxiliary field. We focus on the case of the matrix (or tensor) models, for which such a simple argument does not exist. We clarify the meaning of the effective actions for these models and provide rigorous and elementary derivations of all the required identities. In Sec. 3.3 we introduce the numerical techniques used to solve the Euclidean Schwinger-Dyson equations and provide extensive tests of their precision by accurately reproducing well-known properties of the IR limit of the U(n) 2 model in the SYK-like phase. We then apply the same methods to the U(n) model and discover that its low energy limit, for low enough mass m, is governed by the generalized ansatz (1.3), establishing unambiguously the inequality (1.4). Finally, in Sec. 3.4 we discuss the Euclidean four-point functions. As usual, their connected pieces are obtained by resumming a geometric series associated with a set of ladder diagrams. We provide a diagrammatic derivation of the associated kernels, together with a fully rigorous and general algebraic proof (see also App. B).
Sec. 4 is devoted to the study of the phase diagrams of the models, i.e. their properties in the (m, T )-plane. A general physical discussion is given in Sec. 4.1, where the main features are explained and motivated in simple terms. In Sec. 4.2, we focus on the U(n) 2 model. After recalling and extending the main results already obtained in [29] (line of first order phase transition, critical point, physics at fixed fermion number, low temperature behaviour), we study the q-generalizations of the model (see also App. A). This includes the derivation of the critical exponents for the qgeneralizations of the critical points in Sec. 4.2.3, and the discussion of multi-q models with their associated domain wall solutions and line of critical points in Sec. 4.2.4. The U(n) model is discussed in 4.3. The main feature of its phase diagram is the vanishing of the mass gap at a quantum critical point m = m * > 0. For 0 < m < m * , the model is in the exotic gapless phase with a spontaneously broken scale symmetry. We find strong evidence that the zero-temperature entropy is non-vanishing in this phase and also compute the energy and heat capacity at low temperature.
Real time physics is discussed in Sec. 5. We begin in Sec. 5.1.1 with a discussion of the general properties of real time two-point functions, including Carlson's theorem. We then explain in Sec. 5.1.2 how the Schwinger-Keldysh formalism can be used to constrain the two-point functions in the melonic limit, and how the resulting real time Kadanov-Baym equations can be solved by an elegant argument based on Carlson's theorem, apparently yet unpublished. This provides a rigorous path from the Euclidean to the real time solution of the models. The numerical algorithm used to solve the real time equations, which is an extension of the algorithm used for the Euclidean equations, is presented in Sec. 5.1.3 and used to compute some spectral functions. We then apply our results on two-point functions to the computation of quasi-normal frequencies in Sec. 5.1.4, finding an interesting singular behaviour at the critical point of the U(n) 2 model, with associated non-trivial critical exponents. The problem of the real time four-point functions is tackled in Sec. 5.2. General properties (symmetries, spectral functions, KMS condition) are reviewed in Sec. 5.2.1 and the Schwinger-Keldysh formalism in the melonic case is set up in Sec. 5.2.2. We then discuss the expected behaviour of the real time four-point functions, which strongly depend on the time ordering of the operators, and introduce the chaos Lyapunov exponent η L associated with the out-of-time-ordered correlators in Sec. 5.2.3. We carefully explain in Sec. 5.2.4 how η L can be extracted by simplifying the exact Schwinger-Keldysh equations written down in Sec. 5.2.2. This yields the so-called real time ladder kernels. In Sec. 5.2.5, we present two methods to extract η L numerically from the ladder kernels. The first method is the most direct and rigorous and does not rely on the assumption of any particular ansatz for the correlators. It amounts to inverting a large matrix of discretized real time ladder kernels. The second method is due to Kitaev [6] and is based on a natural ansatz for the long time behaviour of the correlators. It reduces the problem to the computation of the largest eigenvalue of a one-dimensional kernel operator. Both methods are found to yield perfectly consistent results and are applied in Sec. 5.2.6 to the explicit computation of η L in both models. The most interesting new results we obtain concern the non-saturation of the chaos bound in the gapless phase of the U(n) model, and the singular behaviour of η L near the critical point in the U(n) 2 model.

Presentation of the models
The models we shall study can be formulated in two different ways, either as matrix or disordered models. In the matrix models there are N = n 2 d degrees of freedom, organized as a vector of complex fermionic matrices ψ a µ b with U(n) indices a, b and an O(d) index µ. On the other hand, the disordered model variables form a vector of N complex fermions χ i with i = 1, 2, . . . , N . For the purposes of the present work, we restrict ourselves to leading order in a suitably defined N → ∞ limit [3][4][5]. It is a remarkable fact that in this limit both formulations turn out to be strictly equivalent, although this equivalence breaks down if one were to consider 1/N corrections.
The advantage of the disordered models is that they are a priori technically much simpler. Being vector models, they are straightforwardly solved in the case of annealed disorder using the standard auxiliary field techniques. The main and crucial novelty, compared to the old-fashioned vector models [35], is that the auxiliary field is now bilocal [36]. By examining the leading graphs it is also easy to show that the quenched and annealed disorder averages coincide at leading large N order, see e.g. [37] for a detailed discussion.
The advantage of the matrix models is that they make the link with possible applications in holography much more explicit [3], since they are standard quantum mechanical matrix models with the usual string theory interpretation. Indeed, in retrospect one may see the startling equivalence with the disordered models at leading order in N as an explanation for why disordered modelsà la SYK could be related to black-hole physics.
The leading order equivalence between the two points of view deserves to be carefully justified. This essentially amounts to devising proofs that never use the auxiliary field technique, which is not applicable in the matrix case. This will be explained and exemplified in detail along the rest of the paper.
The models we consider are, in some sense to be clarified below, the most general ones based on complex fermions with quartic interaction terms and some simple additional conditions. In this Section we provide a general presentation and exposition of their basic properties.

Large d matrix models
2.1.1 General setup and the large N limit For convenience, we work with fermionic matrices ψ a µ b normalized in such a way that the canonical anticommutator reads The fermion number operator and fermion number are respectively defined bŷ where we use matrix notation for the U(n) indices. With the normalization (2.1), we have that 0 ≤ Q ≤ 1.
We wish to consider the most general single-trace Hamiltonian invariant under U(n) × O(d) with a quadratic mass term and quartic interactions. It is important to distinguish two classes of interaction terms, according to the cyclic order of the O(d) indices in the traces. We will say that interaction terms with alternating indices as in µνµν are of type I or tetrahedric, while those with repeated consecutive indices such as µµνν are of type II or pillow. 1 For example, tr ψ † µ ψ ν ψ µ ψ ν is tetrahedric whereas tr ψ † µ ψ µ ψ † ν ψ ν and tr ψ † ν ψ µ ψ † µ ψ ν are pillows. Packing up terms of each class in h I and h II respectively, we can write the Hamiltonian as Note that the most general quadratic term depends on a unique mass parameter m, because expressions like tr ψ µ ψ µ identically vanish.
We then consider the large n and large d limits defined in the following way [3][4][5]: i) First, we take n → ∞ keeping d and the couplings in h I and h II fixed. This is the usual 't Hooft limit [38], so for example the free energy has a large n expansion of the form Here F g is the sum over genus g vacuum Feynman graphs. We will focus on the leading order F 0 which selects the planar Feynman graphs.
ii) Next, we take d → ∞ again keeping the couplings in h I and h II fixed. Note that the √ d enhancement of the tetrahedric (type I) terms in (2.3) is absolutely crucial here. Without it, the d → ∞ limit would exist but would be trivially vector-like [3,4]. This second limit generates a 1/ √ d expansion of the sums over planar graphs, so for example for the leading-order free energy Once more, in what follows we focus on the leading contribution, i.e. F 0,0 for the free energy and similarly for other quantities.
In the limiting procedure above, the order in which one takes the limits n → ∞ and d → ∞ is crucial. At fixed n, due to the √ d enhancement of the tetrahedric couplings in (2.3), it is easy to find graphs having arbitrarily high powers of d. The large d limit at fixed n is thus ill-defined. On the other hand, taking the large n limit first restricts us to consider only planar graphs, and these have a power of d which is bounded above. Thus the large d limit of the sum over planar graphs is well defined even with the enhanced couplings.
For conciseness, recalling that N = n 2 d we shall refer to the above successive large n and large d limits as "the large N limit" and call the leading contributions in this limit "the leading large N contributions." In this way, the terminology is also appropriate for the discussion of the disordered models presented in Sec. 2.2, so it is possible to switch back and forth between both equivalent formulations.

Inequivalent models
It is straightforward to check that all the single-trace interaction terms corresponding to a given cyclic ordering of the fields in the trace are equivalent at leading order in the large N limit described above, up to a possible redefinition of the mass parameter m and a c-number shift of the Hamiltonian. For example, Note that in some instances a double-trace term tr ψ µ tr ψ † µ may be produced, but this is also sub-leading at large n. Alternatively, such terms can be altogether eliminated by imposing a tracelessness condition tr ψ µ = 0, which does not play any role at leading order. 2 This remark greatly reduces the number of inequivalent terms we may consider at leading order. For the tetrahedric interactions, noting that tr ψ µ ψ ν ψ µ ψ ν = 0, we are left with three possibilities There are a few more possibilities for pillow interactions, but those that contain the combination ψ µ ψ µ or its complex conjugate cannot contribute to leading order graphs. The simplest proof of this uses the auxiliary field technique, which works for pillow interactions but not for tetrahedric interactions (an explicit example will be given in Sec. 3). Taking this into account, we are also left with the three pillows Treating these three terms by introducing auxiliary fields ψ µ ψ † µ and ψ † µ ψ µ , one may show that their effect is to renormalize the mass and to shift the energy by a simple function of the fermion number Q. These changes are rather trivial and can be taken into account straightforwardly. For this reason, we shall discard the pillow terms in most of the following discussion, except for a few specific remarks.
The action of U(n) L × U(n) R × O(d) is defined in matrix notation by The diagonal subgroup for which U L = U R is simply denoted by U(n), and the fermion number symmetry is the U(1) Q subgroup acting by a phase as ψ µ = e −iα ψ µ .
All the terms we are considering are O(d)-symmetric. The terms h I,1 , h II,1 and h II,2 are U(n) L × U(n) R -symmetric and thus also conserve the fermion number. The terms h I,2 and h II,3 violate U(n) L × U(n) R but preserve the fermion number, and the term h I,3 violates both U(n) L × U(n) R and U(1) Q .
One can also define an anti-unitary time reversal operatorT , acting like the identity on the Fock basis obtained by the action of the creation operators ψ † on the Fock vacuum |0 defined by ψ a µb |0 = 0. The hermiticity of the Hamiltonian implies that the couplings λ, λ , κ and κ must be real, and the corresponding terms are then automatically T-invariant. If complex, the couplings ξ and κ violate T; see, however, the discussion in Sec. 2.1.3.

Generalized melons
There are many ways to depict the Feynman graphs of these models. For example, we can use the standard 't Hooft ribbon graph representation, with an additional middle strand associated to the O(d) index. In the case of U(n) L × U(n) R -invariant models the two strands making up the ribbons are associated with distinct unitary groups and can thus be colored, say in green and red. The faces of the graphs then respect such coloring. In this case a so-called "colored graph representation" is also possible, and is used in the tensor model literature (see e.g. [3,4]). For more general models, in which the U(n) L × U(n) R symmetry can be broken down to U(n), a more economic and convenient representation is as follows.
The propagator is represented by a single oriented line. We choose conventionally the orientation to represent the flow of U(1) Q charge, with ψ † µ creating and ψ µ destroying one unit of fermion number. At the vertices, the cyclic ordering of the lines is important: going around a vertex clockwise corresponds to the positions of operators in a single trace from left to right. Moreover, the line in the propagator may be seen as being associated with the O(d) index. Tetrahedric and pillow interaction vertices are thus distinguished by the fact that the lines cross or not, respectively. The vertices associated with the terms in (2.7) and (2.8) are depicted in Fig. 1.
Following the terminology of [4], the graphs that contribute at leading order in N are called generalized melons. More detailed comments on these leading graphs will be given in Sec. 2.3. Here we limit ourselves to describing them.
When only tetrahedric interactions are used, the leading order graphs contributing to the two-point functions are obtained by inserting on bare propagators the generalized elementary melons depicted in Fig. 2, and repeating this process recursively. Note that on top of the "symmetric" elementary melons (λ, λ), (λ , λ ) and (ξ, ξ * ), there exist mixed structures like (λ, ξ * ), etc. These mixed structures, inserted once, can only contribute to the fermion number violating two-point functions tr ψ µ ψ µ and tr ψ † µ ψ † µ , but multiple insertions do contribute to tr ψ µ ψ † µ as well. When pillow interactions are included, we can also insert recursively on bare propagators the elementary melonic structures depicted in Fig. 3. A generic leading order Feynman graph thus has a tree-like structure built out of bubbles attached to each other by a pillow vertex, on which arbitrary recursive insertions of generalized elementary melons may be added. Leading order graphs contributing to the vacuum amplitude, or free energy, are obtained by gluing the external legs of a graph contributing to tr ψ µ ψ † µ at leading order. A sketch of such a leading vacuum graph is depicted in Fig. 4.
Further simplifications result from restricting to the models for which the fermion number violating two-point functions vanish at leading order, tr ψ µ ψ µ = tr ψ † µ ψ † µ = 0. This yields one model for which we can keep the interactions h I,1 and h I,2 and another model based on h I, 3 . But the structure of the elementary melons implies that, at leading order in N , the only effect of adding h I,2 on top of h I,1 is to renormalize the coupling λ → λ + 2λ . It is thus enough to keep only the term h I,1 . We are left with two basic models, As already noted, pillow interactions may also be added to either of these rather trivially. The models (2.10) and (2.11) are those on which we shall focus in the rest of the paper. In spite of very similar-looking Hamiltonians, their physics, in particular their phase diagrams and infrared behaviours, will turn out to be drastically different.
Note that the model (2.10) has the full U(n) L × U(n) R symmetry and in particular conserves the fermion number. The model (2.11), on the other hand, only has the diagonal U(n) symmetry and violates both U(1) Q and time-reversal invariance (for a general complex coupling ξ). However, at leading order the symmetric structure of the elementary melons implies that time reversal invariance is actually restored, because only the combination ξξ * = |ξ| 2 can enter. Moreover, for the same reasons a Z 4 subgroup of U(1) Q acting as ψ → iψ is also preserved at leading order, explaining why tr ψ µ ψ µ and tr ψ † µ ψ † µ vanish. Because of their symmetry properties, we shall refer to the models with Hamiltonians (2.10) and (2.11) as the U(n) 2 and U(n) models, respectively. Note that the U(n) 2 model is the matrix version of the complex tensor model mentioned in [16]. The U(n) model does not seem to have been introduced before, either in tensor or matrix versions.
There are also interesting generalizations involving higher order interaction terms, similar to the q-body interacting SYK model introduced in [9] and based on the complete interaction bubble studied in [4]. We shall briefly comment on such a generalization for (2.10) in the following sections, see also App. A.

Disordered models
In the disordered model framework, the degrees of freedom are N complex fermions conveniently normalized to satisfy the anti-commutation relations and we consider a general Hamiltonian of the form The coupling κ mimics the pillow interactions of the matrix models. The other couplings satisfy the hermiticity conditions 14) and the symmetry conditions They are drawn randomly from a Gaussian distribution with weight (2.16) The large N limit is defined by taking the real couplings κ, λ, ξ and ζ fixed. In this limit, the so-called self-averaging observables, which include the thermodynamical potentials and the two-and four-point functions studied later in the paper, coincide with their average over disorder. This quenched average is usually performed by introducing replicas, but at leading large N order this turns out to be unnecessary. Indeed, the solution is given by the diagonal replica ansatz, which means that the quenched average coincides with the annealed average that is obtained by directly integrating over the random couplings as if they were dynamical variables themselves. The models can then be solved using the standard auxiliary field technique, as will be briefly reviewed in Sec. 3.
It is easy to check that at leading order in N , and for all the purposes of the present paper, the random Hamiltonian (2.13) is equivalent to (2.10) or (2.11) when turning on the coupling λ or ξ, respectively. The effect of the matrix model pillow terms are also correctly mimicked by the coupling κ in (2.13).
Note that differences do exist between the matrix model and the disordered model formulations. For instance, the coupling ζ in (2.13) has no analogue in the matrix case, because a term like tr ψ µ ψ ν ψ µ ψ ν identically vanishes. We will therefore take ζ = 0 in the rest of this paper. Conversely, some features of the matrix model, such as non-vanishing fermion number violating correlation functions, cannot be reproduced by (2.13) because the annealed average automatically yields an emergent U(N ) symmetry which implies fermion number conservation. Equivalently, the integration over the random couplings can produce "symmetric" elementary melons, similar to the ones depicted in the first two rows of Fig. 2, but cannot yield the other asymmetric structures. To achieve this, one needs to generalize the disordered model. A simple possibility is to consider two species of fermionic vectors χ i andχ i transforming in complex conjugate representations of U(N ), using a Gaussian integration measure over random couplings that includes terms mixing different couplings. An explicit example could be a term λ ijk l χ † i χ † jχ k χ l + ξ i jkl χ † i χ j χ k χ l + H.c. in the Hamiltonian with a contribution of the form λ ijk l ξ l ijk + H.c. in the Gaussian weight for the couplings. This does generate asymmetric melonic structures and non-zero fermion number violating two-point functions of the form χ iχ i , but will not be considered further in this paper.
To summarize, the U(n) 2 matrix model is mimicked at leading large N order by the random Hamiltonian 3 which has already been considered in the literature and is often called the "complex SYK model," see e.g. [20,21,[30][31][32]. As for the U(n) model, it is on the other hand mimicked by and we are not aware of any prior work on it, even though it has obvious similarities with (2.17) and might be expected to also be relevant for condensed-matter purposes.

Discussion and relation to other tensor models
This subsection provides details and comments on the matrix model point of view that are not strictly necessary to follow the main new results of the paper, as presented in the next sections. It may thus be skipped upon first reading. The aim is to provide further motivation for our work, as well as to place it in the wider context of recent developments in tensor models. The reader can refer to [41][42][43] for reviews on the latter.
We shall first clarify the nature of the double large n and large d expansions we are considering, since some confusion seems to exist in the literature about certain terminology. We will discuss models with the full U(n) L × U(n) R symmetry in Sec. 2.3.1, and consider models with only one U(n) symmetry in Sec. 2.3.2. Finally, in Sec. 2.3.3 we elaborate on the string-theory inspired interpretation of the models we are studying, as well as the relation between these and other tensor models considered in the literature.

Degree expansion, index expansion, melonic and non-melonic interactions
We have already mentioned in Sec. 2.1.3 that diagrams for models with U(n) L ×U(n) R symmetry can be "colored". As reviewed in detail for instance in [4], both interaction vertices and Feynman graphs have a representation in terms of regular colored graphs, the so-called "bubbles." To any bubble we can associate a Gurau degree [1,41,44] which is in general a non-negative half integer. The original and standard definition of a melon is to be a regular colored graph of Gurau degree zero.
Since interaction vertices are in one-to-one correspondence with regular colored graphs, they can be classified into two groups: the melonic and non-melonic interactions. In our case, the tetrahedric term h I,1 in (2.7) is non-melonic, whereas the pillow terms h II,1 and h II,2 in (2.8) are melonic. To define a large N limit, one should choose the way the coupling constants scale with N . In the case of matrix models, there is a unique interesting choice first discussed by 't Hooft [38], which yields the genus expansion. In the case of tensor models, the situation is much more complex and interesting. Several inequivalent and non-trivial scalings can a priori be defined. In any given model, the optimal scaling of a coupling constant is the one for which it is impossible to further enhance that coupling and still have a well-defined expansion. Computing the optimal scaling in tensor models is in general a difficult technical problem, see e.g. [45].
There are two general procedures known to define scalings with well-defined large N limits for arbitrary interactions in tensor models. 4 One is the so-called BGR scaling [46], which was used originally in the tensor model literature but can also be straightforwardly generalized to the matrix-tensors by distinguishing two sets of indices with d = n, [4]. This scaling has two fundamental properties: i) The Feynman graphs in the large N expansion of the tensor models are classified according to their Gurau degree. In particular, the leading graphs are the standard melons, whose structure was elucidated in [47]. Subleading orders may be studied using the results in [48].
ii) Only melonic interaction vertices can contribute to leading order graphs. In particular, the tetrahedric terms in (2.7) do not contribute at leading order in this framework. This is due to the fact that the √ d enhancement that we have introduced in (2.3) is not present in the BGR scaling.
The physics of quantum mechanical models based on the BGR scaling is similar to the physics of the large N vector models [35] and in particular cannot describe black-hole physics (for instance, there is no continuous spectrum at large N in this scheme). Note also that, in the BGR scaling, if we distinguish U(n) and O(d) indices the large n and large d limits commute.
The other scaling was introduced in [4], generalizing to any tensor (or matrixtensor) model the scaling used by Carrozza and Tanasa for a O(n) 3 tensor model in their seminal work [49]. The scaling of [4] is actually optimal for a large class of interactions, called maximally single-trace, which generalize the concept of singletrace matrix model interactions to general tensor or matrix-tensor models. This new scaling is the one used in the present paper, and is such that: i) The Feynman graphs in the large N expansion of the tensor models are classified according to their index, which is a non-negative half-integer naturally associated to the colored graph. In particular, the leading graphs have index zero and are called generalized melons in [4]. Note that a graph of Gurau degree zero automatically has index zero, but the converse is not true: there are many more generalized melons than melons, and their full classification is at the moment unknown. Some special cases, which include the interactions used in the present paper, have been treated in [4,49].
ii) Both melonic and some non-melonic interactions can appear in the leading order graphs. In particular, both the tetrahedric h I,1 vertex in (2.7) and the pillows h II,1 and h II,2 in (2.8) contribute at leading order, which makes the physics much richer and more interesting than that of the BGR scaling.
Even though the generalized melons built out of the tetrahedric interaction vertex h I,1 in (2.7) have a strictly positive Gurau degree, they share the basic fundamental recursive structure associated with standard melons: they are built by arbitrary repeated insertions of generalized elementary melons, depicted in Fig. 2, and the sum over such graphs can be evaluated using Schwinger-Dyson equations of "melonic" type. This simple but non-trivial result is known to be valid for the interactions used in the present paper, as well as for some generalizations thereof [4,5]. In more general cases, it is still usually very easy to show that graphs built in such a way will contribute at leading order. However, let us emphasize that the hard part of the proof is to demonstrate that no other graph can contribute; the simple "melonic" recursive structure is not a priori required. This hard part is often overlooked in the modern literature. The full classification of generalized melons remains a very interesting open problem.

U(N )-symmetric models and the tracelessness condition
Almost all the tensor model literature, including [4], relies on the fact that each distinct index of the tensors (or of the matrix-tensors) is associated to a distinct symmetry group. This simplifies the analysis of the Feynman graphs considerably, because different indices of a given tensor cannot mix in a given face.
For a long time it was believed that this simplifying property was necessary to define a consistent large N limit. The fact that this is actually not the case was first pointed out in [3]. A simple argument, based on the possibility to color the faces of the relevant planar graphs with only two colors, implies that the large d limit with the enhanced scaling of [4,49] used in the present paper is consistent and that the graphs contributing at leading order are again the generalized melons. The models containing terms that break U(n) L × U(n) R down to U(n), like (2.11), can thus have a well-defined large N limit.
An obvious question to ask is what happens beyond the planar graphs. To some extent this question may be purely technical or academic, because for the most important physical applications we are aware of, only the planar graphs are needed. Moreover, in the matrix interpretation the higher genus terms are also subleading to all the d 1− /2 , > 0, corrections that already contribute at the planar level. These planar corrections are themselves out of reach with present technology, except for the leading piece in √ d which is easy to include [49]. 5 This being said, it is an extremely interesting problem from the graph-theoretic point of view, and is not limited to the matrix-tensor models. From the purely tensorial point of view, it is very natural to impose symmetry properties on the indices by considering, for instance, completely symmetric or antisymmetric tensors. Such symmetry conditions have the same effect as the U(n) L × U(n) R symmetry breaking terms in the models we study, the result being that distinct indices mix up in faces and the usual tensor model technology based on colored graphs does not apply.
It was very soon realised after [3] that in these models there exists a specific class of non-planar graphs that makes the new large d limit inconsistent, see [5] for explicit examples. Similar graphs spoil the large N limit of other (anti)symmetric tensor models. An elegant way to remove these graphs, which have a very specific structure, is to impose tracelessness conditions [50]. This was shown to be correct in the case of symmetric tensors in [51], using a technical tour de force based on a set of local moves simplifying the graphs [52].
The problem remains open in the case of the matrix-tensors. The additional difficulty is that we do not have complete symmetry between all the indices in this case, and in particular the tracelessness condition is imposed only on the matrix indices. 6 A precise conjecture for the models of interests in the present paper is as follows.
Conjecture: For the models (2.3) for which we impose the condition tr ψ µ = 0, the genus g free energy (and similarly any other observable) has a well-defined large d expansion of the form where η g is an upper bound on the highest possible power of d for a graph of fixed genus g. From [3], we already know that η 0 = 1. Generalizing results proven in [3,4] for the case of U(n) L × U(n) R models at higher genera, for which the tracelessness condition is not required, we may futher conjecture that η g = 1 + g.

Matrix models vs. tensor models, brane picture and gauging
As discussed above, the models we are considering are close cousins of the fermionic tensor models introduced in [15,16]. To make the relation more explicit, one can set d = n and treat the ψ a µ b as a three-index tensor to take the large n limit. In the present context, this procedure seems rather artificial because the indices a, b and µ are associated with different symmetry groups. We believe it is more physical and interesting to keep the matrix model point of view. In the matrix model picture, one can always interpret the indices a, b associated with the unitary symmetry as Chan-Paton factors of oriented open strings ending on some kind of D-particle. The additional orthogonal O(d) symmetry is then naturally interpreted as corresponding to the rotational invariance in the directions transverse to the D-particle worldline. This picture allows interesting physical interpretations of the results, and makes the relation with string theory, holography and black-hole physics much clearer. The large d limit might also have an interesting relation with the large dimension limit of general relativity studied in [53].
At subleading order, the difference between the matrix model and the tensor model interpretations becomes quantitative.
One reason is that the corrections are organized differently in the two pictures. In the matrix models, planar contributions are always dominant over non-planar contributions, even if we consider higher order corrections in 1/ √ d. On the other hand, if d = n as in the tensor model interpretation, a given order in the large n expansion can mix up matrix model graphs having different genera.
Another reason is that the gauging of the global symmetries is done in a different way in the matrix models and in the tensor models. In the former, in line with the open string picture, we gauge only the U(n) or U(n) L × U(n) R symmetries that act on the matrix indices. The O(d) symmetry remains ungauged and is associated with the conservation of angular momentum. One then gets a standard string-like spectrum of operators. On the other hand, in the tensor model picture it is more natural to gauge the full global symmetry group, U(n) 2 × O(n) in our case, or O(n) 3 in well-studied examples [54].
Note, however, that both in large d matrix models and in tensor models the gauging is relevant only at next-to-next-to leading order, because the dimension of the gauge group is of order n 2 . Any contribution coming from the gauging, for instance a gauge-fixing action and a Fadeev-Popov determinant, will thus contribute at order n 2 , whereas the leading and next-to-leading orders are n 2 d and n 2 √ d (or n 3 and n 5/2 in tensor models). For instance, any Hagedorn-like phase transition at weak coupling will be pushed to zero temperature at leading order (in large d matrix models, it is easy to check that the weak coupling Hagedorn temperature is of order 1/ ln d).

Euclidean time picture
We now proceed to study the models (2.10) and (2.11) in the Euclidean formalism, leaving for Sec. 5 the corresponding analysis in real time. Detailed derivations of various formulas are provided for thermodynamic quantities, as well as for two-and four-point functions. The low energy solutions to the relevant Schwinger-Dyson equations are spelled out, and we also discuss the numerical methods used to solve these equations at finite coupling, exemplifying their application.
An important point we discuss is the possibility to have non-trivial monodromies in parameter space for the two-point functions. This means that the Schwinger-Dyson equations can have several distinct and consistent solutions for fixed values of their parameters. This possibility was first discovered in [29] and later also observed in [55] for a related system with two coupled SYK models. It will play a major role in Sec. 4 when we study the phase diagram of the U(n) 2 model (2.10).
An innovative aspect of our discussion is the introduction of a new ansatz for the infrared behaviour of two-point functions. We assume that their power-law decay is governed in general by two distinct exponents ∆ + and ∆ − when the Euclidean time goes to +∞ and −∞ respectively. The usual ansatz used so far in the literature corresponds to ∆ + = ∆ − , but this is not strictly required in certain cases. The possibility to have ∆ + = ∆ − will play a crucial role in elucidating the physics of the U(n) model (2.11).
We also provide a careful justification of the use of effective actions for the standard bilocal fields G(τ 1 , τ 2 ) and Σ(τ 1 , τ 2 ) in the matrix (or tensor) models. Most of the literature has focused on disordered models, in which this effective action can be straightforwardly derived by a path integral argument, see e.g. [56], but this is not applicable in the case of matrix or tensor models.
Finally, we derive a general formulaà la Bethe-Salpeter (3.133) for the kernel used in the computation of the four-point functions, which is applicable in a wide class of models. Its analytic continuation to real time will be performed in Sec. 5, allowing us to study the strength of quantum chaos as diagnosed by the Lyapunov exponent.

General definitions and properties
We define the time-ordered finite temperature Euclidean two-point function by The Euclidean time evolution of any operator is defined as usual by The equations (3.1) are valid for −β < τ < β. The correlators satisfy the KMS condition and can be extended to all Euclidean times by antiperiodicity of period β = 1/T , In particular, we use the Fourier decomposition with Matsubara-Fourier coefficients G k and Matsubara frequencies The tree-level Matsubara coefficients are given by yielding the tree-level two-point function where we have introduced the Fermi distribution function On top of (3.3), the Euclidean two-point function has a set of interesting general properties. 7 To understand these properties, it is very useful to introduce the so-called spectral function, defined by where Z = Tr e −βH (3.10) is the partition function, E p are the eigenvalues of the Hamiltonian, the double sum in (3.9) is over all the eigenstates |p and |q of the Hamiltonian 8 and The two-point function has the following spectral decomposition, which is equivalent to The spectral function ρ is manifestly real and positive. This yields and from (3.12) and (3.13), The moments of the spectral function, govern the UV behaviour of the Matsubara-Fourier coefficients (3.13), and through (3.12) they fix the discontinuities of the two-point function and its derivatives, In particular, the canonical commutation relations (2.1) or (2.12) yield Moreover, as will be demonstrated in the following (see 3.2.2), the first moment ρ 1 is exactly given by a one-loop calculation, and the one-loop contribution is subleading in the N → ∞ limit. At leading large N order, ρ 1 thus coincides with the tree-level resulṫ A last very useful quantity that we need to introduce is the so-called resolvent R, defined by The resolvent is a holomorphic function for Im z > 0 and Im z < 0. Eq. (3.13) is equivalent to The large z expansion of R is governed by the moments (3.17) and in particular When the spectrum of the Hamiltonian is discrete, which is the case in these models at finite N , the only singularities of R are simple poles on the real axis at the Bohr frequencies. At infinite N , the spectrum becomes continuous and R develops a discontinuity across the real axis. The spectral function can be read off from this discontinity as

Schwinger-Dyson equations
Let us introduce the self-energy Σ(τ ), which is an antiperiodic function with Matsubara coefficients This equation is equivalent to where we have introduced the notation δ β for the antiperiodic δ-function The self-energy is computed by resumming the one-particle irreducible (1PI) twopoint graphs from which the full two-point function is obtained by the usual geometric series, The recursive structure of the leading order graphs discussed in Sec. 2.1.3 implies simple self-consistency conditions, called Schwinger-Dyson equations, on the sum of the 1PI two-point graphs. For instance, for the U(n) 2 model (2.10) and the U(n) model (2.11), we get Including contributions from the pillow vertices is very easy. For instance, if we turn on the h II,1 term of (2.8) in the U(n) 2 model, the right-hand side of (3.30) gets an additional contribution −κG(0 − )δ β (τ ) = κQδ β (τ ) due to the possibility of inserting the structure on the left of Fig. 3. The full Schwinger-Dyson equation resulting from this addition is depicted schematically in Fig. 5. The same term would be added if h II,1 were to be turned on in the U(n) model. The term h II,2 produces a contribution −κ G(0 + )δ β (t) = κ (Q − 1)δ β (τ ) and the term h II,3 a contribution (Re κ )(G(0 + ) + G(0 − ))δ β (τ ) = (Re κ )(1 − 2Q)δ β (τ ), consistently with the other structures in Fig. 3.
From the point of view of the disordered models, one does not need to analyze the leading graphs to derive the Schwinger-Dyson equations. One can instead apply a purely algebraic reasoning using the standard auxiliary field method, based on the saddle-point equations of the effective action derived from the path integral.

On the solutions of Schwinger-Dyson equations and monodromies
The Schwinger-Dyson equations are perturbative in nature, since they are associated with the summation of a class of Feynman graphs. However, they also contain nonperturbative information because the sum of Feynman graphs they compute has a finite radius of convergence. Performing the analytic continuation of the perturbative series to strong coupling can yield extremely interesting and non-perturbative physics.
Morally speaking, this is very similar to what happens, for example, in some supersymmetric gauge theories [58], where one resums weakly coupled semiclassical contributions that have a finite radius of convergence to obtain non-trivial results at strong coupling via analytic continuation. In the models studied in [58], these are instanton contributions, not perturbative contributions, but the spirit is essentially the same.
In both the U(n) 2 and U(n) models, the perturbative series in λ 2 or |ξ| 2 can be straightforwardly computed from the Schwinger-Dyson equations. For instance, in the U(n) 2 model we may plug (3.7) into (3.30) to obtain the first non-trivial order of Σ, which is O(λ 2 ). Inserting this approximate solution into (3.29) to get a correction for G k , we find This procedure can be iterated to any desired order in λ. In particular, it leads to Eq. (3.21) since all the perturbative corrections to the tree-level result decay at least as fast as 1/k 3 when k → ∞.
An important feature of the models we are studying is that the perturbative expansion converges in a region of parameter space that includes two physically very different regimes. i) One regime is the high temperature limit at fixed mass and coupling, βm 1 and βλ 1 (in the U(n) 2 model) or β|ξ| 1 (in the U(n) model). This regime is perturbative because the fermions get a large effective thermal mass of order πT . In this limit, the leading order solution is G k = i/ν k , which is equivalent to G(τ ) = 1 2 sign(τ ). At the level of the states, this is an expansion around the maximally entropic state associated with the density matrix 2 −N I.
ii) Another regime is the high mass limit λ/m 1 (in the U(n) 2 model) or |ξ|/m 1 (in the U(n) model). This is the standard perturbative regime in particle physics, whose leading order is governed by the usual harmonic oscillator solution (3.7). In this limit, the models remain perturbative even at zero temperature and the ground state is simply the Fock vacuum |0 .
Overall, there is an effective coupling constant g eff = min(λ/T, λ/m) in the U(n) 2 model, or g eff = min(|ξ|/T, |ξ|/m) in the U(n) model. The perturbative series converges when g eff is small enough.
A crucial property of the solutions of the Schwinger-Dyson equations is that they can undergo non-trivial monodromies in parameter space. This possibility was first realized in [29], and later also found in other disordered models with quadratic terms [55]. In other words, depending on the path one uses to perform the analytic continuation from the weakly coupled region g eff 1 where the perturbative series converges to a point in the strongly coupled region beyond the radius of convergence, one may end up with distinct solutions. These distinct solutions are a priori all physically relevant, in the sense that they satisfy all the general consistency conditions listed in Sec. 3.1.1, namely (3.3), (3.14), (3.15), (3.16), (3.20) and (3.21). This property is illustrated in Fig. 6 and will be abundantly discussed in the following. Low entropy, HO-like perturbative regime High entropy, SYK-like perturbative regime Figure 6: Non-trivial monodromy in parameter space associated with two paths along which the weakly coupled solution at A is analytically continued to B. Units are such that λ = |ξ| = 1 and the dashed line represents the separation between the weakly coupled and the strongly coupled regions. When going first through the high T , high entropy perturbative state and then to the low T , non-perturbative regime (upper path), it is natural to get a solutionà la SYK at point B. When going first through the high mass, low T , low entropy, harmonic oscillator-like perturbative state (lower path), one may obtain a completely different, low entropy solution at B [29].

Reparameterization invariance and standard low energy ansatz
The U(n) 2 model Let us look at the U(n) 2 model in which the term h II,1 is being turned on as well (including other pillow vertices would change the discussion only in a trivial way). The Schwinger-Dyson equations read and we neglect the term iν k and the ultraviolet δ-function contributions on the righthand side of (3.33) and (3.34), we can write the low energy (or long time) limit of the equations in the convenient form These equations are covariant under the transformations Invariance of the equations is achieved for diffeomorphisms such that β = β. The equations are also invariant under the local multiplications where g is an arbitrary non-vanishing function of period β. The emergence of the reparameterization symmetry (3.38)-(3.39) is a remarkable feature of the low energy limit of Schwinger-Dyson equations resumming melonic graphs, consistent with a dual gravitational description of the models [6,9]. The invariance (3.41) is a local Euclidean version of the global fermion number symmetry of (2.10). The fact that a global symmetry becomes local is also beautifully consistent with holographic ideas.
The full low energy reparameterization symmetry (3.38)-(3.39) will always be spontaneously broken. In particular, time translation invariance implies that the physical two-point functions depend only on τ 1 − τ 2 , a property which is of course not invariant under arbitrary reparameterizations. Moreover, Eq. (3.36) involves an integral over all time scales. Since it is a low energy equation, a UV cutoff is typically needed to avoid UV divergences. In practice, the reparameterization symmetry is thus also explicitly broken by this cutoff.
A fundamental property of the models, first discussed in [11], is that some solutions can preserve a SL(2, R) subgroup of reparameterizations, or more generally a SL(2, R) subgroup of the full set of transformations (3.38), (3.39) and (3.41) [21], at low energy. For the U(n) 2 model, the standard zero-temperature ansatz is (3.42) Instead of b + and b − , it is often convenient to use C > 0 and θ defined by Note that (3.15) implies that b ± > 0 and thus |θ| < π/4. Computing the Fourier transformG of G, Eq. (3.42) and (3.43) can be rewritten as In terms of the spectral function introduced in Sec. 3.1.1, this is equivalent to and One low energy parameter remains unconstrained, for example the angle θ. This parameter is naturally related, via a complicated non-perturbative formula, to the mass parameter m in the Hamiltonian, possibly renormalized by the pillow couplings. When m = 0 the model recovers a large N particle-hole symmetry which yields G(−τ ) = −G(τ ), as in the standard SYK model, and thus θ = 0.
Note that (3.36) is equivalent tõ and thus implies in particular thatΣ(0) = 0. Using the first equation in (3.35), we see that the Fourier transform of the self-energy Σ at zero frequency is constrained to satisfyΣ for the SL(2, R)-preserving solution. Satisfying this condition is therefore an indicator that the solution is indeed conformal in the IR.
The scaling form of the two-point function at non-zero temperature can then be obtained from the zero-temperature ansatz by using the mapping τ = tan(πτ /β) and (3.38) to go from a solution defined on the real line (β = +∞) to a solution defined on the interval [−β/2, β/2]. 9 The resulting solution does not satisfy the antiperiodicity G(−β/2) = −G(β/2), but this can be corrected by performing an additional transformation of the form (3.41), with a function g(τ with a damping factor Eq. (3.50) provides a good approximation to the two-point function in the scaling regime βλ 1, with the additional condition |τ | 1/λ which excises the UV region. By taking the Fourier transform, we get the Matsubara-Fourier coefficients , (3.52) which are valid when βλ 1 and |ν k | λ, as well as which is valid when βλ 1 and |ω| λ.
U(n) model A very similar analysis can be performed for the U(n) model using (3.31) instead of (3.30). At low energy, the reparameterization symmetry (3.38), (3.39) is still valid, but (3.41) does not hold anymore, consistently with the fact that the fermion number is no longer conserved. The low energy ansatz (3.42), (3.43) can still be used, and plugging it into the Schwinger-Dyson equations we find that the only possibility with b ± = 0 corresponds to The particular mapping τ = tan(πτ /β) from S 1 β to the real line is singled out by the fact that it is consistent with time translation invariance. We would like to thank Paolo Gregori for providing this argument. Important remark: The existence of a consistent low energy ansatz solving the low energy Schwinger-Dyson equations does not ensure that there exists a physical solution to the full Schwinger-Dyson equations whose low energy limit matches the low energy ansatz. In other words, it may not always be possible to UV complete solutions that seem to be consistent at low energies. Sometimes this occurs for very elementary reasons. For instance, the low energy ansatz (3.42) can a priori solve the low energy Schwinger-Dyson equations for any mass parameter m in the Hamiltonian. But, at least in the perturbative regime m λ, we know that such a solution cannot be UV completed: the perturbative series converges in this regime and the solution must be very near the harmonic oscillator which has a trivial IR behaviour. The lack of a UV completion is more subtle and has important consequences in the context of purely bosonic models [5,59,60]; see also the discussion in Sec. 3.3.2.

A new low energy ansatz
The standard ansatz (3.42) for the long time behaviour of the zero-temperature twopoint function can be generalized to where the exponents ∆ + and ∆ − are a priori distinct. This possibility may seem exotic and, to the best of our knowledge, has never been considered in the literature so far. However, in models for which the particle-hole symmetry is broken it is a plausible and mathematically consistent possibility. One salient new result of the present paper will be to provide strong numerical evidence that this behaviour is actually realized in the U(n) model.
An important new feature of the ∆ + = ∆ − case is that the two terms on the right-hand side of (3.31) scale differently at τ = +∞ and τ = −∞. Assuming, for example, that ∆ + < ∆ − , Eq. (3.31) implies that 10 The τ → +∞ term then dominates in the Fourier transform at low energy. Using the Schwinger-Dyson equations (3.36) and (3.56) as before, and keeping the leading terms at low energy, we get with no further constraints on ∆ − and b − . In passing, we note that the relation (3.49) must also be satisfied by such a solution.
Let us remark that both the standard (3.42) and the more general (3.55) ansätze spontaneously break the reparameterization symmetry (3.38)- (3.39). However, the scale symmetry, which is preserved by (3.42), is spontaneously broken by (3.55) when ∆ + = ∆ − . This is an important qualitative difference between the two classes of solutions.
It is natural to seek a finite temperature solution generalizing (3.50) using the reparamaterization covariance of the Schwinger-Dyson equations. However, the usual method does not work. First, the invariance under local multiplications (3.41), which is crucially used in the U(n) 2 model to obtain a finite temperature two-point function with the correct antiperiodicity property, is no longer valid in the U(n) model. Second, and even more crucially, the mapping τ = tan(πτ /β) yields a two point function that is not time-translation invariant, except if ∆ + = ∆ − = 1/4. 11 Because of these difficulties, we shall limit ourselves to the zero-temperature formula (3.55) when On the practical side, it is interesting to observe that, even though the exponent ∆ + in the case ∆ + < ∆ − and the exponent ∆ in the case of the standard ansatz ∆ = ∆ + = ∆ − are both equal to 1/4, the two situations can be sharply distinguished even if we look only at the τ → +∞ asymptotics of the two-point functions. Indeed, there is a factor of 8 discrepancy for b 4 + between (3.54) and (3.57).

Disordered model formulation
In the case of the disordered models, one can compute the Euclidean finite temperature path integral for the Hamiltonian (2.13) straightforwardly. This is a completely standard calculation and thus we shall be brief.
Since the annealed and quenched averages coincide at leading order, as already mentioned in Sec. 2.2, we first integrate over the disorder with the Gaussian weight (2.16). This produces a vector model with bilocal interaction terms which can be treated in the standard way [35]. We introduce two bilocal auxiliary fields G and Σ, where Σ plays the role of a Lagrange multiplier enforcing the constraint G(τ 1 , τ 2 ) = χ i (τ 1 )χ i (τ 2 ). Integrating out the complex fermions yields the effective action where the operators O Σ and O 0 , which act on antiperiodic functions, are defined by Note that, formally, the tree-level terms − ln(1 + e −βm ) and ln det O 0 in (3.58) cancel each other. However, we prefer to use the rigorous formula (3.58), which does not have UV divergences and takes into account that only the ratio of determinants det O Σ / det O 0 is actually well defined.
The large N limit is governed by the saddle points of (3.58). Varying this effective action with respect to Σ and G, denoting the saddle point values of Σ and G by Σ and G respectively, and using time translation invariance, we get  (3.34) to the case where the couplings λ, ξ and κ are turned on simultaneously. 12 In particular, and not surprisingly, the on-shell auxiliary field G coincides with the Euclidean two-point function.
Moreover, by construction the on-shell value of the effective action, From the free energy one can then obtain the entropy and the energy by using the standard thermodynamical identities, The above results are very powerful: they imply that the knowledge of the Euclidean two-point function only is sufficient to derive all the thermodynamical properties of the model. In practice, we solve the Schwinger-Dyson equations numerically as described in Sec. 3.3, and plug the result into (3.58) to compute F , S and E straightforwardly. Very explicit formulas are given in Sec. 3.2.3 below.

Matrix model formulation
The case of the matrix models is a bit more subtle to deal with, because there is no longer an auxiliary field method available to compute the large N path integral. It is thus not immediately obvious how to compute the free energy. In particular, it is unclear whether an effective action like (3.58) could still be relevant. The aim of this subsection is to present a detailed derivation of the free energy and to carefull justify the use of effective actions like (3.58) in the context of matrix models. 13 We explained in Sec. 2.1.3 that turning on both couplings λ and ξ in the matrix model yields complications due to non-vanishing fermion-number-violating two-point functions. For simplicity and because we won't need the more general results in the rest of this paper, we are then going to focus on the models (2.10) and (2.11). Note however that our method can be generalised straightforwardly to the models in which both couplings are turned on simultaneously. Note also that, as usual, pillow couplings may be turned on quite trivially by using the auxiliary field method.
The derivations for (2.10) and (2.11) are extremely similar. We thus present the derivation in parallel for both models.
Step one The Heisenberg Euclidean equations of motion follow from (3.2) and read In the U(n) 2 model, this yieldṡ Contracting with ψ † µ on both sides and taking the trace we get (3.67) 13 The arguments apply to tensor models as well, see also [61].
In particular,Ġ Let us note en passant that, subtracting the above two equations and using (3.20) and the canonical anticommutation relations to evaluate we getĠ This together with (3.19) yields an exact formula for the first moment of the spectral function, generalising (3.21) at finite n and d.
In the U(n) model, we have 14 Massaging a bit the traces, we geṫ In particular, in the U(n) modelĠ and Eq. (3.21) is valid even at finite n and d.
Step two The partition function Z and free energy F are defined in general by where Tr traces over the Hilbert space of states. Then in the U(n) 2 model and in the U(n) model. From (3.68) and (3.73), neglecting the subleading 1/ √ d term in the U(n) 2 model, we thus geṫ (3.81) In particular, since G is real, in the U(n) model we must have This result is obvious from the structure of the generalized melons in the large N limit, but we see here that it is valid even at finite n and d. In the U(n) model, we can thus rewrite (3.81) aṡ Finally, note that, by definition of the self-energy (3.26), we can also writė From the point of view of the Fourier transform, the first equality in the above equation follows from the fact that both G k and Σ k are proportional to 1/k at large k and thus the series k∈Z+1/2 G k Σ k converges.
Step three Define the functionals which, using the Schwinger-Dyson equations (3.30) or (3.31), as well as the antiperiodicity of G, can be rewritten as Comparing with (3.81), (3.84) and (3.85), we get for the U(n) 2 and U(n) models, respectively.
Step four Integrating (3.92) and fixing the constant of integration by using the fact that at zero coupling both models are simply a collection of N = n 2 d decoupled fermionic oscillators, we obtain (3.63). This is exactly as in the case of the disordered model formulation, where S eff is replaced by S for the U(n) 2 or U(n) models, according to which couplings are turned on.

Useful formulas
Numerically, the solution of the Schwinger-Dyson equations is most naturally encoded in the Matsubara-Fourier coefficients G k (or Σ k ). It is thus very convenient to express all the relevant physical quantities, such as fermion number, free energy, energy and entropy, as convergent sums over these coefficients.

Fermion number The fermion number (2.2) is given in terms of the two-point function
One cannot use immediately the expansion (3.4) to evaluate G(0 ± ) because k G k diverges. This is related to the fact that the function G is discontinuous at τ = 0. To go around this problem, one can use Dirichlet's trick of considering G(τ ) + G(−τ ), which is continuous at τ = 0 and has a convergent Fourier expansion. This yields Using then (3.93), (3.20) and (3.14), we get Subtracting the tree-level result, with tree-level coefficients G (0) k given in (3.6), we also get Formulas (3.95) and (3.96) are of course mathematically equivalent, but (3.96) is more convenient to use in numerical evaluations because the infinite series converges more quickly. This is due to the fact that Re G k and Re(G k −G (0) k ) decay like 1/k 2 and 1/k 4 respectively at large k. We shall repeatedly use a similar rewriting of the formulas in the following, separating the tree-level result and the quantum corrections.
There is also an interesting formula for the fermion number in terms of the spectral density, which is obtained from (3.93) and (3.12), This formula coincides with the formula for a collection of decoupled fermionic harmonic oscillators with a frequency distribution given by the spectral density ρ.
Free energy The free energy is obtained by evaluating the effective actions (3.86) or (3.87) on the solution for the two-point function and using (3.63). It is convenient to use the Schwinger-Dyson equations (3.30) or (3.31) to simplify the result. We get in this way Note that this formula is valid for both the U(n) 2 and the U(n) models.
Energy The energy is simply the expectation value of the Hamiltonian, which is valid for both the U(n) 2 and the U(n) models. The right-hand side of (3.100) can be evaluated using (3.95) and (3.85), yielding Subtracting the tree-level result yields the equivalent formula to be used in the numerical evaluations.
It is interesting to express the energy in terms of the spectral density, using (3.100) and (3.12) together with (3.20) and (3.21). We obtain Note that, unlike the case of the fermion number (3.97), this formula differs slightly from the formula for a collection of decoupled fermionic harmonic oscillators, by the replacement of the frequency factor ω by 1 2 (m + ω).
Entropy The entropy is then given by the standard relation We have also checked explicitly that the result so obtained is consistent with the first relation in (3.64).
Zero temperature in the U(n) 2 model One can derive very interesting exact formulas at zero temperature in the U(n) 2 model [13,20] when the low energy ansatz discussed in Sec. 3.1.4 is valid. The charge at zero temperature can be expressed in terms of the low energy parameter θ defined in (3.43) as Moreover, if we denote by σ the zero-temperature entropy, is the so-called Catalan constant. When m = 0, it is useful to consider σ as a function of the zero-temperature charge Q(T = 0, m) instead of m. One can then derive the elegant formula Using (3.43), (3.105) and (3.107) as an initial condition, 15 and after a little massaging, this equation can be integrated to or more explicitly The above zero-temperature formulas will not play a crucial role in our analysis of the model, but will be used in the next subsection to check the precision of our numerical results.

Numerical solution of Schwinger-Dyson equations
In this subsection, and in all numerical computations except otherwise stated, we always choose "natural units" such that the couplings are set to one, i.e. λ = 1 and |ξ| = 1 .
We introduce a high-frequency cutoff k UV and restrict our attention to the 2k UV Matsubara frequencies labelled by k = −k UV + 1 2 , . . . , k UV − 1 2 . This effectively reduces (3.113) to a system of 4k UV algebraic equations for the coefficients G k and Σ k . We solve this system using a standard iterative procedure which, given an approximate solution (G k through (3.113). The weighting factor α can a priori be an arbitrary real number, and this choice will directly affect the convergence of the algorithm. One may refine it by making α frequency-dependent, or updating it from one iteration to the next as in [9]. For most purposes of the present paper, we found that a fixed, frequency-independent value α = 0.75 is sufficient. We found no other solutions by using other values of α, by introducing a frequency dependence or by updating the weighting factor during the iterations.
The above algorithm is a generalization of the iterative procedure commonly used to solve fixed-point equations of the form f (x) = x, of which (3.113) can be thought of as being a multivariate version. Solutions x * to f (x) = x can be found from the iteration ) as long as the initial guess for the solution is close enough to x * and |f (x * )| < 1. When this condition is satisfied, we say that the fixed point is attractive. If |f (x * )| ≥ 1, on the other hand, the fixed point is repulsive and the method cannot converge. This is where the use of a weighting factor α can be helpful. Indeed, all fixed-points of f (x) are also fixed-points of f α (x) = αf (x) + (1 − α)x, and it is easy to check that we can always find a non-zero value of α such that |f α (x * )| < 1. Thus, a repulsive fixed-point of the original problem f (x) = x may be found as an attractive fixed point of the modified f α (x) = x problem, for some suitably chosen value of α.
An important ingredient of the algorithm is the choice of the initial approximation. It may be taken to be the tree-level coefficients (3.6), G k , or alternatively any other previously found Matsubara-Fourier coefficients deemed to be close to the actual solution. The iteration is then repeated until the equations are satisfied up to a prescribed tolerance ε, i.e.
Since we work with 64-bit processors, we chose the tolerance to be ε = 10 −12 in the units (3.112), safely above machine-precision.
There is no guarantee that this algorithm terminates, but if convergence is achieved in I α (ε) steps then the running time is O(I α (ε) × k UV log k UV ). The k UV log k UV factor in this expression comes from performing two Fast Fourier Transforms (FFTs) in each step, first to obtain k and then to obtain Σ , or more generally O(k UV log q) for the q-body generalizations of (3.113) (see App. A), so it is always subleading because we shall always have q k UV .
We implemented the above algorithm in C++, and also in Mathematica for reference and double-checking purposes. In our C++ implementation, FFTs are calculated using the Cooley-Tukey algorithm [62] optimized to work with k UV being a power of 2. In general, lower temperatures require higher values of k UV to obtain good solutions. The reason is that the highest Matsubara frequency we include must always be much greater than the scale set by the coupling constant, i.e.
The condition (3.116) ensures that at scale ν k UV −1/2 we are well into the UV free asymptotic regime of the models, governed by the expansion (3.18) with (3.20) and (3.21), for which After achieving convergence, one should always check that (3.117) is satisfied, as well as the basic consistency conditions (3.14), (3.15) and (3.16).
In practice, the cutoff k UV also affects the precision with which we can compute thermodynamic quantities through formulas such as (3.96), (3.98) or (3.102). Indeed, introducing a high-frequency cutoff amounts to the truncation of the corresponding series, thus resulting in a systematic error in the computed values. To keep this error under control at temperatures of the order T ∼ 0.1, one finds that k UV = 2 10 is high enough, but for T ∼ 0.001 the cutoff should be increased at least to k UV = 2 16 , and for T ∼ 0.0001 we can reach k UV = 2 21 . See App. C for a more thorough discussion.
We shall use the above algorithm extensively in the following. For instance, when computing the phase diagrams in Sec. 4 we map the solutions that exist at each point in the space parameterized by (m, T ). In order to do this we may begin by finding numerical solutions of the Schwinger-Dyson equations in a region where an approximate solution is known to be valid, e.g. in one of the perturbative regimes T 1 or m 1. The algorithm then converges if we use as initial solution the tree-level coefficients (3.6), G k . Keeping all but one of the parameters fixed, we can then find another solution close to this point using the previously found solution as an initial approximation, successively repeating this process to explore all of parameter space. While doing this, it is very useful to keep track of the thermodynamic properties of the solutions and the convergence of the algorithm, in order to adjust the step sizes accordingly.
One extremely important feature, that we have already mentioned in Sec. 3.1.3, is the possible existence of multiple solutions at a given point (m, T ), due to nontrivial monodromies in the space of solutions. This is signalled by the fact that the process described in the paragraph above can sometimes result in hysteresis curves for the thermodynamic functions. A representative example is depicted in Fig. 7. This means that the solution we find numerically depends not only on the parameters (m, T ) explicitly entering the equations (3.113), but also on the initial approximation G k we used to start the iterative procedure. In this way, the solution "remembers" the specific sequence of steps in parameter space that was taken to reach it, and through it also the perturbative regime where the whole process typically started. This phenomenon will be discussed extensively in Sec. 4.
Note on numerical errors and error bars 16 As discussed in detail in App. C, our numerical results are subject to two main sources of errors: systematic errors due to the introduction of cut-offs or discretizations; and methodological errors related to our fitting methods, etc. We do not have rigorous control over systematic errors, S/N m Figure 7: Hysteresis curve for the entropy of the U(n) 2 model as a function of mass, for fixed temperature T = 0.035 in natural units (data points shown with crosses for illustrative purposes). The red curve is obtained starting from the solution found at m = 0 and gradually increasing the mass, whereas the blue curve is obtained by decreasing the mass starting from the solution found at m = 0.6. In a range of masses around the value m ∼ 0.26, two distinct solutions to the Schwinger-Dyson equations clearly coexist.
which typically may become important at very low temperatures. The data analysis methods are however under better control, and error bars can be associated with them as explained in App. C. The error bars that appear in the following therefore always refer to the methodological errors, but do not take into account the systematic errors that may appear on top of them. This implies that the actual values of some low-temperature observables may lie outside our error bar intervals. This issue is only relevant in some specific, yet important, instances and will be mentioned and discussed explicitly whenever it occurs.

Basic numerical results for the U(n) 2 model
The numerical algorithm described above may be used to study the finite coupling Schwinger-Dyson equations (3.113) for the U(n) 2 model, and assess for example the applicability of the low energy ansatz discussed in Sec. 3.1.4. From previous studies of the complex SYK model [20,21] we can expect that this should be the case for low enough masses, whereas at large mass the system must be in a gapped perturbative phase (see Sec. 4 for more details).
In Fig. 8 we plot some solutions found for m = 0 and m = 1 at various tempera-tures. The m = 0 solutions match with the standard SYK model, and we clearly see that the finite-temperature conformal ansatz (3.50) becomes a very good approximation in its range of validity, which in the natural units (3.112) is when β 1, τ 1 and β − τ 1. The m = 1 solutions display a distinctly different behaviour, with the two-point function decaying exponentially. An energy gap can be immediately extracted from the log-plot of this type of solutions, see Sec. 4. For m = 0 the model loses its particle-hole symmetry and thus G(β − τ ) = G(τ ). In Fig. 9 we plot solutions for various intermediate values of the mass at T = 0.01. The marked change in behaviour around m ∼ 0.3, from power-law decay at low masses to exponential decay at large masses, will be associated in Sec. 4 to a first order phase transition.
A very rough but simple criterion to distinguish the low and large mass solutions is to test the condition (3.49). Working at finite inverse temperature β the lowest frequency we have access to is ν 1/2 = π/β, so that we can use Re Σ 1/2 as an approximation toΣ (0) There are different ways in which we may quantitatively test the low energy ansatz of Sec. 3.1.4, beyond the naked-eye impression provided by Fig. 8 or the very rough testing of (3.49) through Re Σ 1/2 . One possible approach starting from the zerotemperature expression (3.42) consists in extracting the scaling dimension from a The precise interval of time we use in practice may be defined in different ways, but fortunately this choice does not greatly affect the fitted value of the scaling dimension ∆ unless obviously unreasonable. Comparing the results obtained with different choices can serve the dual purpose of making sure that the numerical results make sense, as well as providing an estimation of the error or uncertainty associated to the computed values, see App. C for more details.
One important advantage of the procedure just described is that we can repeat it independently for τ → −∞. We then obtain from each numerical solution not one but two values for the scaling dimension, denoted by ∆ ± for τ → ±∞ as in (3.55). These values could a priori be different, but we find that in the U(n) 2 model they are always very similar, and certainly equivalent within the numerical precision we can achieve. This is in sharp contrast with the case of the U(n) model discussed in the next subsection. Applying the same methods, we will clearly see there that ∆ + = ∆ − , signalling the fact that the behaviour of the U(n) model indeed follows the non-standard ansatz (3.55).
More precisely, the β → ∞ limit of the finite temperature formula (3.50) predicts We thus expect a linear dependence in τ , with a vanishingly small slope in the β → ∞ limit. This is illustrated in Fig. 10 for T = 10 −3 and m = 0.1. We indeed find a linear behaviour within the range of Euclidean times considered, both for τ → +∞ and τ → −∞. The scaling dimensions ∆ ± are then obtained from the τ -independent part of the corresponding linear fits, and seen to be consistent with the expected value given by (3.46), ∆ = 1/4. We can also read the coefficient a from the slope and then derive the angle θ from (3.51). Note that this determination uses a subleading term and thus cannot be very accurate.
Datapoints for τ < 0 Datapoints for τ > 0 This result is consistent with a = a + = a − = 0.04 3 and through (3.51) yields an angle θ = 0.12 9 . A more precise determination of a and θ is obtained in Fig. 11.
We can also use the previously found values for ∆ ± as input to fit the coefficients b ± directly in G(τ ), again taking the numerical solution restricted to the finite intervals considered before. This is illustrated in Fig. 11 for the example with β = 1000 and m = 0.1, and provides in particular a better determination of the parameters a and θ through (3.51). A non-trivial check of the whole procedure is the fact that these satisfy with very good accuracy the relation (3.47). Another complementary approach is to fit ∆ and b ± , with a = 1 2π ln b + b − as before, by using the natural generalization of (3.42) to arbitrary ∆ given in App. A, Eq. (A.7). This ansatz can be trusted as soon as |τ | 1, without the need to restrict to |τ | β as before. It is therefore not surprising that the results obtained in this way are even more precise. Some examples are shown in Fig. 12, for solutions with m = 0.1 at various temperatures. The fitting interval is delimited by vertical dotted lines, but once more the results are fairly robust with respect to its variations. Repeating this process at a given fixed mass for various temperatures, we observe that the scaling dimension varies linearly with the temperature for low enough temperatures. Such a linear dependence could certainly be explained by looking at subleading corrections to (3.42), but we have not tried to do so. We may consider the appearance of such linear behaviour as an indicator that the IR regime has been properly reached, since it is then immediate to extrapolate the finite temperature results to T → 0.
The consistency of the results convincingly demonstrates that the solution to the Schwinger-Dyson equations approaches the low energy ansatz discussed in Sec. 3.1.4. Similar results are obtained for other small masses, with the approximation becoming better and better as the temperature gets lower. Note that the aim of the detailed analysis of some standard cases, as discussed above, is mainly to test our numerical methods before applying the same ideas to the much more interesting U(n) model. In that case the non-standard ansatz (3.55) will play a crucial role, and we do not have a finite temperature expression analogous to (3.50). We will therefore not have the luxury of added precision. However, the consistency observed in the U(n) 2 model between the different methods provides us with a high degree of confidence that the results for the U(n) model that we present in what follows are also very reliable.

Basic numerical results for the U(n) model
We can perform on the U(n) model a numerical analysis similar to the one done in the previous section for the U(n) 2 model. Note that for m = 0 we can use the particle-hole symmetry to look for solutions satisfying G(−τ ) = −G(τ ), exactly as in the U(n) 2 model. The Schwinger-Dyson equations are then equivalent to the standard SYK model [2,6,9], the two-point function following at low temperatures the standard ansatz of Sec. 3.1.4.
For m = 0 the situation is very different already at a qualitative level. To illustrate this, we show in Fig. 13 solutions of the Schwinger-Dyson equations found for various masses at T = 0.01 in natural units. Again we expect the system to develop an energy gap for large enough values of the mass, and indeed we see that the two-point function is exponentially decaying already for m = 0.5. However, we observe in this case a gradual, smooth transition into this regime, contrasting with the sharp change of behaviour seen for the U(n) 2 model, e.g. in Fig. 9. Figure 13: Numerical solution of the Schwinger-Dyson equations for the U(n) model at β = 100 and |ξ| = 1, for various values of the mass m. The loss of the particle-hole symmetry is signalled by the asymmetric shape of G(τ ) about τ = β/2. The plot is consistent with a smooth change in the long time behaviour of the two-point function, from power-law to exponential decay, to be contrasted with that of the U(n) 2 model shown in Fig. 9.
Further evidence in favor of a smooth transition into the perturbative regime is provided by the constraint in (3.49). Note that this condition is equally valid for the standard and non-standard ansätze discussed in Sec. To elucidate in more details the nature of the gapless regime, we can try to extract scaling dimensions from the numerical solutions at finite temperature, as done before for the U(n) 2 model. Having fixed the mass to some small value, a scaling dimension growing unboundedly as the temperature decreases would be a telltale sign that the IR behaviour of the two-point function is trivial. On the other hand, a linear dependence of ∆ with the temperature as T → 0, as seen in the U(n) 2 model, would be strong evidence in favor of some non-trivial scaling in the IR.
A representative case with m = 0.1 at T = 10 −3 in natural units is illustrated in Fig. 14. An exponential decay of the two-point function is clearly ruled out, the results being consistent with a power-law decay both when τ → +∞ and when τ → −∞. However, as already anticipated the exponents ∆ ± relevant at τ → ±∞ must now be distinct, in sharp contrast with the situation observed in the U(n) 2 model in Sec. 3.3.1 . In other words, the general ansatz (3.55) appears to be valid, but the standard ansatz (3.42) is unambiguously violated. More precisely, as in Fig. 10 for the U(n) 2 model we find an approximately linear behaviour of d log |G(τ )| d log |τ | when 1 |τ | β and τ > 0. This allows for the determination of an exponent ∆ + with satisfactory precision, the result being consistent with the expected value of 1/4, see Eq.(3.57). On the other hand, linearity is not as good when τ < 0, suggesting that the inverse temperature β = 1000 is not high enough to properly reach the deep IR regime. This yields a significant uncertainty in the determination of ∆ − , but in any case the equality ∆ − = ∆ + is ruled out. 17 We have also proceeded to extract the coefficients b ± using the zero-temperature ansatz (3.55), as shown for the U(n) 2 model in Fig. 11. This yields a rather imprecise determination b + = 0.64 5 and b − = 0.46 7 , again supporting the interpretation that the deep IR regime has not been fully reached.
Further evidence in favour of the generalized ansatz ∆ + = ∆ − is found when analyzing the temperature dependence of the various quantites involved. In Fig. 15 we illustrate the typical situation at low enough masses by plotting the results obtained for m = 0.1. The scaling dimension ∆ + extracted for τ → +∞ is seen to develop a very clear linear behaviour at low temperatures. This allows us to confidently extrapolate its value to T → 0, obtaining ∆ + = 0.251 1 . On the other hand, the linear regime does not seem to be reached for ∆ − . Taking into account the sign of the small non-zero second derivative in the plot, our results should then be interpreted as providing a lower bound on ∆ − , establishing unambiguously that ∆ − > ∆ + . A similar analysis is made in Fig. 16 for the coefficients b ± . The extrapolation of b + , which is observed to be quite reliable, is of particular interest. Indeed, it provides an additional criterion to distinguish the standard ansatz ∆ + = ∆ − , for which b + = (1/4π) 1/4 0.53 from Eq. (3.54), from the generalized ansatz ∆ + = ∆ − , for which b + = (2/π) 1/4 0.89 from Eq. (3.57). We find in this case b + = 0.82 4 , a result clearly in favour of ∆ + = ∆ − and consistent with our previous discussion of the scaling dimensions.
In conclusion, we have discovered a new IR behaviour in the U(n) model, inconsistent with the standard ansatz (3.42) discussed in Sec. 3.1.4 and used so far in the literature. It is however in good agreement with the more general, non-standard ansatz of Sec. 3.1.5, at least to the degree in which we can trust the numerical methods involved. More details, in particular a more global picture of the resulting phase diagram, will be provided in Sec. 4, whereas the consequences of this new IR behaviour in the real time physics shall be presented in Sec. 5.
(3.121) The four-point function also satisfies the reality and symmetry conditions (3.124) The last condition above follows from time-reversal invariance, which is always valid in the leading large N melonic limit we are interested in. 18 It is natural to decompose φ into a disconnected and a connected piece, Once more the linear behaviour is observed to be more accurately followed for τ → +∞ than for τ → −∞, the determination of b − consequently being more imprecise.
Our goal is to compute the connected piece F at leading order in the large N limit. We first present a diagrammatic analysis in the next subsection and then a rigorous and very general algebraic proof in Sec. 3.4.3.

Diagrammatic analysis
Graphs contributing to F at leading order are easy to find. One considers ladders of the form depicted in Fig. 17, with rungs of various types. These ladder graphs can be easily resummed because they form a geometric series. If we introduce a convenient matrix notation for functions of four variables, with a pair of variables forming a "matrix index" and matrix multiplication rule we can write the geometric series for F in the form The function F 0 (τ 1 , τ 2 , τ 3 , τ 4 ) = −G(τ 1 , τ 4 )G(τ 3 , τ 2 ) (3.128) corresponds to the contribution with no rung inserted, whereas the kernel K is built from all the possible rung structures that can appear in the ladder. Figure 17: Ladder graph contributing to the connected four-point function F. For better readability, we use a simple oriented line to represent the full two-point function (instead of a line with a shaded disk as in Fig. 5). The shaded rectangles can be any of the rung structures depicted in Fig. 18 for the U(n) 2 model, or in Fig. 19 for the U(n) model. In the U(n) 2 model, the kernel is given by The first term on the right-hand side of (3.129) is associated with the left and central rungs in Fig. 18, while the second term corresponds to the rightmost rung depicted in that figure.
In the U(n) model, the kernel is similarly given by The first term on the right-hand side of (3.130) is associated with the three rungs in the first row of Fig. 19; the second term with the three rungs in the second row of Fig. 19; and the third term with the six rungs in the last two rows of Fig. 19. Figure 19: Rungs from which we can build the ladder diagrams contributing to the leading order of the connected four-point function of the U(n) model. For better readability, we use a simple oriented line to represent the full two-point function.
Even though it is completely straightforward to check that the ladder diagrams do contribute at leading order to F, it is less immediate to prove diagrammatically that they are the only diagrams having this property. This part of the proof is often omitted in the literature. In the next subsection, we provide a very general algebraic method that permits the rigorous computation of the connected four-point function kernels of a large class of models. Applied to the U(n) 2 and U(n) models this method yields (3.129) and (3.130), respectively, thus proving that the ladder diagrams considered above are indeed the only ones that can contribute at leading large N order.

Algebraic proof
We shall use the language of matrix models to describe the proof of the formulas (3.129) and (3.130) for the four-point function kernels, but the disordered model formulation can be treated in exactly the same way. Only some very general features of the models are needed, namely: i) The large N limit is dominated by generalized melonic graphs such that the Schwinger-Dyson equations for the two-point function G and self-energy Σ coincide with the saddle-point equations of an effective action functional S eff , of the general form ii) The partition function Z is given at leading order by the on-shell value of the effective action, where S eff = S eff (G, Σ) as in (3.62).
The Euclidean four-point function (3.120) is then given by (3.125), (3.127) and (3.128), with a kernel defined by the general formula To prove this formula, it is convenient to introduce the source formalism. The partition function Z has a standard path integral representation, with matrix Grassmannian integration variables ψ µ andψ µ and an action (3.134) The term S int represents the interactions in the Hamiltonian. We can modify this action by adding a source term for the two-point function. The new action is Equivalently, we can express G and the connected piece F of the four-point function defined in (3.125) in terms of the free energy, so that Because the source term in (3.135) is quadratic in ψ andψ, its addition only modifies the expression of the propagator, without changing the large n and large d diagrammatics of the models. The results (3.131) and (3.132), which are a priori assumed to be valid without the source term, can thus be immediately generalized to take it into account. This amounts to replacing the mass term mδ β (τ 1 − τ 2 ) in O Σ and O 0 defined in (3.59) by mδ β (τ 1 − τ 2 ) + J(τ 1 , τ 2 ). The form (3.131) of the effective action then implies that this modification yields a new, source-dependent effective actionŜ eff [G, Σ; J] given bŷ For the particular case of the models we are discussing, one then finds the explicit form of the kernels (3.129) and (3.130) by plugging in (3.133) the relevant functionals for s[G], namely s U(n) 2 and s U(n) given by

Phase diagrams
In this section, unless explicitly stated otherwise, we continue using for convenience the natural units defined in (3.112). Having set λ = 1 (for the U(n) 2 model) or |ξ| = 1 (for the U(n) model), we then have two dimensionless parameters, m and T , and we are going to study in detail the properties of the models at large N in the (m, T )-plane.
We can assume m ≥ 0 without loss of generality. Indeed, the sign of m may be changed by permuting "particles" and "holes," i.e. by exchanging the role of the operators ψ and ψ † (in the matrix models) or χ and χ † (in the disordered models).

General discussion
A qualitative picture of the physics to be expected in the (m, T ) plane is illustrated in Fig. 20. Our goal in this subsection is to elaborate on this picture and elucidate many important features, without resorting to difficult calculations or numerical analysis. Full details are then provided in the following subsections 4.2 and 4.3.
Zero mass When the mass strictly vanishes, the models have particle-hole symmetry. Assuming that this symmetry is not spontaneously broken, a possibility that we shall not investigate here, we can seek a solution to the Schwinger-Dyson equations such that G(−τ ) = −G(τ ), and this solution matches the standard solution for the original SYK model studied extensively in the literature [6,9,[63][64][65]. The only difference with SYK is the doubling of the number of degrees of freedom, because we use complex and not Majorana fermions. This introduces trivial factors of two, for example in the formula for the entropy. The physics at m = 0 is thus familiar. The models flow from a trivial perturbative behaviour at high temperature, governed by The physics at small mass is SYK-like, with a flow from a perturbative regime at high T to a strongly coupled conformal phase at low T . The physics at large mass is harmonicoscillator-like, weakly coupled for all values of the temperature. Our main goal below is to elucidate and study the consequences of the non-trivial transition from the weakly coupled low-temperature region at large mass to the strongly coupled lowtemperature conformal region at small mass.
the standard high T perturbation theory for fermionic systems around the maximally entropic state with density matrix 2 −N I, to a non-trivial, strongly coupled, conformal behaviour at low temperature, reminiscent of extremal black-hole physics. In particular: the zero-temperature Euclidean two-point function G has a power-law decay |τ | −2∆ at large times, with ∆ = 1/4; there exists a non-vanishing zero-temperature entropy given by (3.107); and the chaos exponent, as read-off from real time out-oftime-ordered four-point functions, is 2πT at low temperatures, saturating the bound discussed in [34].
Small non-zero mass When the mass is small but non-zero, three qualitatively distinct behaviours could a priori be expected.
The first possibility is that a mass gap is immediately created, as would be required if the physics were weakly coupled. This is what happens, for instance, when two copies of the SYK model are coupled quadratically as in [55,66].
The second possibility is that the physics remains gapless when the mass is turned on. This is what happens in both of the models we consider, over a finite interval of masses 0 ≤ m ≤ m * , where m * is a strictly positive critical mass beyond which a mass gap is created. The nature of the deformation induced by the mass term is, however, qualitatively different in the U(n) 2 and the U(n) models. This was mentioned in Sec. 3.3.1 and 3.3.2 and will be thoroughly discussed in Sec. 4.2 and 4.3. More precisely, in the U(n) 2 model the mass term flows to a marginal deformation in the IR, in the sense that the IR physics smoothly depends on m. In the U(n) model, on the other hand, the mass operator is relevant at strictly m = 0, in the sense that the deep IR physics is governed by two disconnected fixed points at m = 0 and m > 0. It becomes marginal for all 0 < m < m * .
The third possibility is that the mass term is irrelevant in the IR. This behaviour would actually be quite natural in the U(n) model, since the standard low energy ansatz (3.42), (3.43) has no free parameter in this case, see Eq. (3.54). However, our numerical results rule out this possibility in favor of the most general behaviour described in Sec. 3.1.5.
Large mass At very large masses the models are weakly coupled for all values of the temperature, and are governed by the standard perturbation theory around a set of N decoupled harmonic oscillators. The physics is then very simple. There is a mass gap m eff > 0 of order m, with the zero-temperature Euclidean two-point function G decaying as e −m eff τ at large times; the zero-temperature entropy vanishes, the ground state being the unique Fock vacuum; and the chaos exponent is o(T ) at low temperatures.
In the U(n) 2 model, the fermion number is conserved and the mass term can then also be interpreted as (minus) a chemical potential for the fermion number. Under these circumstances, it is straightforward to derive, for instance from the spectral decomposition formula (3.12), that in a phase with a mass gap the mass dependence of the zero-temperature two-point function is exactly given by a simple exponential term. More precisely, for m 1 and m 2 any two masses in the gapped regime where the ground state is the Fock vacuum. This very general result is realized in the simplest possible way in the U(n) 2 model. Indeed, the tree-level zero-temperature two-point function, is actually an exact solution of the zero-temperature Schwinger-Dyson equations for any value of the coupling, simply because the right-hand side of (3.30) involves the product of G at opposite times.
In the U(n) model this simplification of the large mass physics does not occur. The fermion number is not conserved and the mass is not a chemical potential. Even though the theory is still gapped and the ground state is the Fock vacuum, the simple relation (4.1) is not valid. Indeed the zero-temperature two-point function can get non-trivial corrections, due to the presence of the G(τ ) 3 term on the right-hand side of (3.31). In particular, the physical mass gap m eff , which is computed in Sec. 4.3, is a non-trivial function of m.
The transition from weak to strong coupling The above discussion implies that, when the mass is decreased from large to small values, the behaviour of the models must drastically change from a gapped perturbative regime in which the ground state is the Fock vacuum to a non-perturbative regime with no gap. Therefore, there exists a critical mass m * below which the mass gap closes, m eff > 0 for m > m * , m eff = 0 for m < m * . (4. 3) The zero-temperature transition at m = m * is reflected in the zero-temperature energy, which vanishes in the Fock vacuum at m > m * and must be strictly negative at m < m * . The transition also comes along with the creation of a non-vanishing zero-temperature entropy. This follows from (3.110) for the U(n) 2 model and will be convincingly demonstrated numerically for the U(n) model in Sec. 4.3.
However, in spite of these similarities between the two models, we can infer from our discussion of the high mass regime that the transition at m = m * must be qualitatively different in the U(n) 2 and U(n) models.
In the U(n) 2 model, because the zero-temperature two-point function matches with the tree-level result (4.2) in the gapped regime, the mass gap is exactly m eff = m for all m > m * . The transition from m > m * to m < m * must then be discontinuous, with in particular a discontinuous jump of the mass gap from m * to zero. Such a first order phase transition was indeed found in [29]. Its existence is made possible by the non-trivial monodromies in the space of solutions of the Schwinger-Dyson equations that we have already mentioned in Sec. 3.1.3. The first order transition is not limited to zero temperature, but extends up to a critical point at finite temperature [29]. More details are given in Sec. 4.2.
In the U(n) model, the zero-temperature two-point function gets non-trivial corrections at m > m * and the mass gap is a non-trivial function m eff (m). We shall provide convincing numerical evidence in Sec. 4.3 that m eff smoothly vanishes when m → m * from above. We thus get in this case a large N version of a quantum critical point. In particular, the first order phase transition is replaced by a smooth crossover at finite temperature.
Most of the rest of this paper will focus on providing more details on the transition from weak to strong coupling in both models, using the standard thermodynamical Euclidean point of view in Sec. 4.2 and 4.3 and a more modern and innovative real time point of view in Sec. 5.
The gravitational collapse interpretation There is an interesting qualitative interpretation of the transition from weak to strong coupling when the mass is varied at a fixed (low) temperature in the models we are studying, based on the compelling analogy between black-hole physics and the strongly coupled conformal phase on the one hand and the matrix model formulation of the models on the other hand.

As already mentioned in Sec. 2.3.3, quantum mechanical matrix models can always be interpreted in terms of open strings and D-particles. In this point of view, the operators ψ †a
µ b create open strings whose endpoints are attached to the D-particles labelled by the U(n) indices a and b. Since the mass of a string is proportional to its length, the mass parameter m represents the distance between the D-particles joined by the string.
At large mass, we thus have an assembly of D-particles which are very far away from each other and interact very weakly. Their dynamics is governed by standard matrix model perturbation theory. In the usual terminology, this is the regime where the "brane picture" is reliable. As the mass is decreased, the D-particles are brought closer and closer together and the matrix model becomes more and more strongly coupled. Eventually, we expect the system to undergo gravitational collapse and form a black hole. It then enters a regime in which the "brane picture" is no longer suitable, and should be replaced by a "gravity picture" in which the branes are effectively described by a non-trivial gravitational background containing a black hole.
It is very natural to identify this phenomenon with the expected non-trivial transition at m = m * discussed above. This is perfectly consistent with the appearance of a macroscopic zero-temperature entropy for m < m * , signalling the formation of a black-hole horizon, as well as with other properties of the strongly coupled regime of both models we consider.
Of course, the matrix models under consideration are still quite far from being truly realistic models of black holes. For instance, they do not contain fundamental bosonic matrix degrees of freedom, which are most naturally associated with transverse dimensions of space-time. Moreover, the strongly coupled phases of the models are "black holes" only in the same sense as the SYK model is one: there is a very suggestive matching of various physical properties at the qualitative level (entropy, quasi-normal behaviour, continuum spectrum, chaos) and some non-trivial quantitative agreement (saturation of the bound on the Lyapunov exponent). However, it is also clear that they can at most represent "string size" black holes for which a purely Einsteinian description cannot be fully adequate. In particular, the models contain an infinite tower of states in the bulk that cannot be decoupled from the low mass modes [9]. Moreover, the strongly coupled phase of the U(n) model is quite exotic and its interpretation as a "black hole" is then even more questionable than in the case of the U(n) 2 model (see Sec. 4.3).
Keeping in mind these important caveats, we nevertheless expect that the transition from weak to strong coupling in these models does capture important features of the transition between the "brane description" and the "gravitational description" of a collection of D-branes, i.e. of the phenomenon of gravitational collapse. To the best of our knowledge, such a description has never before been given in a fully quantum mechanical framework, even in a very simplified set-up such as the one presented here.

Phase diagram of the U(n) 2 model
The phase diagram of the U(n) 2 model was first described in [29], and is depicted in Fig. 21. We will briefly review here its main features to offer a self-contained discussion, and present additional details about its construction, generalizations and most important properties. The phase diagram's structure is governed by the existence of a critical temperature T c below which two distinct and consistent solutions to the Schwinger-Dyson equations may coexist. The critical temperature itself can be determined numerically using binary search, taking as division criterion the appearence of the hysteresis phenomenon already discussed in Sec. 3.1.3 and exemplified in Fig. 7. The resulting value T c = 0.067 1 is small but finite, and is compatible with theoretical back-of-the-envelope estimates. Having determined the critical temperature, the phase diagram of Fig. 21 was constructed by taking ∼ 130 fixed temperatures in the range T ∈ [10 −3 , T c ], and sweeping the (m, T )-plane from small to large masses and vice-versa, as explained in Sec. 3.3.
One of the solutions found below the critical temperature T c is the analytic continuation of the weakly coupled solution at large masses. We shall call this solution "harmonic oscillator-like", or simply HO-like solution for short. It can be analytically continued down to a minimum mass m − (T ), but ceases to exist for any mass m < m − (T ). As discussed in Sec. 4.1, at strictly zero temperature this solution matches with the tree-level result, which solves the Schwinger-Dyson equations for all masses, so that m − (T = 0) = 0. At finite temperatures, however, we always find m − (T ) > 0.
A different solution existing for T < T c is the SYK solution at m = 0, which in turn can be analytically continued to higher masses up to the maximum value m + (T ). This solution is conformal in the IR and retains many of the familiar properties of the SYK model, such as non-vanishing zero-temperature entropy, so we shall call it the SYK-like solution.
The curves m − (T ) and m + (T ) are drawn in Fig. 21 with dashed lines. For any fixed temperature T < T c , there is an interval of masses m − (T ) < m < m + (T ) for which both solutions coexist. Their free energies coincide at a transition mass m t (T ), the HO-like and SYK-like solutions being more stable for m > m t (T ) and m < m t (T ) respectively. The system thus undergoes a phase transition when crossing the line m t (T ), the jump from one solution to the other being discontinuous. The corresponding phase transition is therefore of first order. This is illustrated in Fig. 22, where the thermodynamical potentials are plotted as a function of mass for temperatures T < T c , T = T c and T > T c .
As we have mentioned at the end of Sec. 4.1, the first order phase transition observed by decreasing m at fixed T , for T < T c , is reminiscent of the phenomenon of gravitational collapse. On the other hand, if we start from very low T in the mass interval m * < m < m c and increase T , then we also cross the line of first order phase transitions. In this case, the phenomenon is somehow reminiscent of the Hawking-Page phase transition, where the thermal AdS space becomes thermodynamically unstable compared with the large Schwarzschild-AdS black hole solution [67].
A very important feature of the phase diagram is that the mass interval in which the two solutions coexist shrinks as the temperature is increased, and eventually reduces to the point m = m c when T = T c . At this point, the two solutions coincide and the transition becomes of second order. Remarkable properties of this critical point and its generalizations will be discussed below in Sec. 4

.2.3 and 4.2.4. For
T > T c there is only one solution and the continuation from large to small mass is completely smooth.
In [55], Maldacena and Qi studied a system consisting of two copies of the SYK model coupled via a quadratic "mass" term. The resulting phase diagrams, depicted in Fig. 11 and 18 of Ref. [55], have both similarities and differences with the one described here.
An important similarity is that there are two competing solutions in a region of the (m, T )-plane, with an associated first order phase transition. Moreover, it seems that On the other hand, a crucial difference is that the mass parameter m in Maldacena-Qi is relevant: a mass gap is created as soon as the mass is turned on. All the solutions for the two-point functions in their model are thus gapped when m = 0, in parallel with our HO-like solutions; in other words, they do not have SYK-like solutions with a non-trivial IR behaviour when m = 0.
Another interesting aspect of the quadratically coupled SYK model is that it can be fruitfully studied in the large q limit, where q is the order of the interactions. In this limit a third branch of solutions is found, and allows to interpolate smoothly between the two other types of solutions responsible for the existence of the first order phase transition that we have discussed above. Solutions in this third branch have a negative specific heat. Their existence is irrelevant in the canonical ensemble, but Maldacena and Qi argued that they could smoothen the transition in the microcanonical ensemble. 19 Note that this novel class of solutions was not found by solving numerically the Schwinger-Dyson equations at finite q, so it might be an artifact of the large q limit, even though this seems unlikely.
The U(n) 2 model we are describing here can also be generalized to arbitrary qbody interactions, see App. A and the results in Sec. 4.2.3 and 4.2.4. However, the large q limit, when carried out in a standard way along the lines of [9] and [55], turns out to be rather trivial. In particular, the phase transition is never seen. The reason for this seems to be that the critical masses m * or m c , or the critical temperature T c , become exponentially small with q at large q, instead of being governed by a simple power law. For this reason, we shall not discuss the large q limit any further here. Note also that we have not been able to find numerically any different solution with negative specific heat at finite q either, although we can of course not guarantee it does not exist.

The (m, Q) plane and physics at fixed Q
The phase diagram in the (m, Q) plane is depicted in Fig. 23. One notable feature is the existence of a large forbidden region, shaded in grey. Another notable feature is that for temperatures T < T c there is an interval of charges that can only be realized if the two phases coexist (vertical dotted lines on the figure). 20 In the context of most condensed matter work, it is natural to interpret the mass term as minus a chemical potential, which is adjusted to fix the fermion number Q to any desired value. This is not the point of view adopted in the present paper, where we always consider the mass to be the most natural variable to be fixed. An interesting consequence of our analysis is that for any given temperature T < T c , there is a maximum mass m + (T ), associated with the dashed line on the right in Fig.  21, above which the system cannot exist in the SYK-like phase.
The existence of this maximum mass, and thus of a minimum chemical potential −m + (T ), is a very well established fact. Independently of the numerical results, we can provide the following simple argument. Let us assume that the m + (T ) bound does not exist, so that the SYK-like solution can then be analytically continued anywhere in the (m, T ) plane. Starting from the grey-shaded region in Fig. 21 where both the SYK and HO-like solutions exist, and performing the analytic continuation of the two solutions to the high-temperature regime along a path that goes around the critical point (m c , T c ) on its right, we deduce that two distinct solutions must exist in the weakly coupled, high T at fixed m regime. Moreover, these two distinct solutions must have the same UV behaviour G k ∼ i/ν k given by (3.18), since this behaviour cannot change under analytic continuation. But this is impossible. At high temperatures, all the Matsubara frequencies are in the UV regime ν k λ and thus any solution with the standard UV behaviour will necessarily be very near the perturbative solution G(τ ) 1 2 sign(τ ). A unique series expansion for any such solution can be generated from the Schwinger-Dyson equations; this series coincides, of course, with the high-temperature perturbative series. It converges and yields a unique solution, thus proving that the bound m + (T ) necessarily exists.
The numerical results suggest that for the SYK-like phase the lines Q(m) at fixed temperature T < T c all terminate at a strictly positive charge Q min (T ) > 0. This minimum possible charge Q min (T ) corresponds to the endpoints of the upper dashed lines in Fig. 23, and can be computed at any given temperature evaluating Q(m) at m = m + (T ). For instance, we find at zero temperature Q min (T = 0) = 0.150 2 . Note that Q min is of course less than the value of the charge below which the solution becomes thermodynamically unstable; at T = 0, this latter charge corresponds to the upper left endpoint of the black line in Fig. 23 and is found to be approximately Q * = Q(T = 0, m = m * ) = 0.243 1 . Our numerical results seem to indicate that Q min (T ) is a decreasing function of T which goes to the value Q c = 0.050 1 at the critical point when T = T c .
However, it is also conceivable that our numerical algorithm may fail to converge for a small interval of masses near m + (T = 0), for which the slope |∂Q/∂m| is large; one may then imagine that the zero-temperature line Q(m) (the red dashed line in Fig. 23) might extend all the way down to zero charge. On the other hand, at non-zero temperature it seems very unnatural to obtain solutions with zero charge and finite chemical potential. Then Q min (T ) should be strictly positive, at least for all T > 0, while going to Q c when T = T c (since at the critical point the SYK-like solution merges with the HO-like solution).

The physics at zero temperature
It is certainly very interesting to study the phase diagram at zero temperature, but the numerical algorithm used to solve the Schwinger-Dyson equations operates inherently at finite-temperature. Of course, the identification of the two distinct solutions in the U(n) 2 model as SYK-and HO-like already constitutes a statement about the low-temperature regime. However, this identification is possible using only very elementary data, such as the value of the entropy (see Fig. 7) or the rough criterion associated to the constraint (3.49). Indeed, as shown in Fig. 24 the value of Re Σ 1/2 at a small but finite temperature T = 10 −3 is approximately equal to the mass for the analytic continuation of the m = 0 solution, whereas it essentially vanishes for the analytic continuation of the large mass solution.
In order to provide further details about the zero-temperature physics, a careful analysis involving more delicate extrapolations is required. We have already shown in Fig. 12 an example for how such a procedure may be carried out for some small but non-vanishing mass, its results matching the expectations coming from the lowtemperature ansatz, (3.46) and (3.47). A more global picture is provided in Fig. 25 by systematically repeating that process for various masses. We clearly see that the agreement is very good all along the range of masses for which the SYK-like phase actually exists.
A further non-trivial check of the extrapolation procedure is made possible by the exact T = 0 results for the SYK-like phase reviewed in Sec. 3.2.3. In Fig. 26 we plot the T → 0 extrapolations of the charge and entropy as functions of the IR parameter θ, along with the exact results (3.105) and (3.111) for comparison. The agreement between numerical and theoretical values is very good for masses up to m = m * and a bit higher, but starts exhibiting deviations in the entropy for larger masses. We can attribute these discrepancies to numerical artifacts becoming more and more relevant as the mass increases, the entropy being very sensitive to the truncation implemented  by the frequency-cutoff of the numerical algorithm (see App. C). In any case, these demonstrations of the precision of the numerical extrapolation used to implement the T → 0 limit give us confidence that our results are very robust at least up to m ∼ m * . Moving on, in Fig. 27, we plot the zero-temperature energy, entropy and specific heat divided by T , as functions of the mass, for the SYK-like solution and the most stable solution. Note that all these quantities are trivially vanishing for the HO-like solution, since it coincides with the tree-level result at zero temperature. The thermodynamical potentials can be expanded in terms of , σ and c at low temperature as The results are consistent with the general discussion in Sec. 4.1. The ground state jumps from the Fock vacuum at m > m * = 0.213 2 to a highly degenerate state at m < m * . The heat capacity c grows from its c(m = 0) = 0.80 1 value matching the result reported for SYK [9], and experiences a finite jump to zero at m = m * . While the extrapolation for this quantity should certainly experience similar numerical artifacts to those present in the entropy (see Fig. 26b), it seems plausible that c will diverge when m → m + (T = 0). We also plot in Fig. 27d the effective mass obtained from the exponential decay of the HO-like solution, extrapolated to zero temperature, finding perfect agreement with the m eff = m prediction. This result will serve as a baseline reference when studying the U(n) model in the following.
Remark : There is no contradiction between the fact that the mass gap m eff jumps discontinuously from m eff = m for m > m * to m eff = 0 for m < m * (Fig. 27d)  fact that the zero-temperature energy is a continuous function of m at m = m * (Fig.  27b). This is a large N effect: the eigenstates of the Hamiltonian responsible for the small negative energy when m is just below m * do not contribute to the large N two-point function when m > m * .

The q-generalizations and critical exponents
The U(n) 2 model admits a straightforward generalization described in App. A, for which the quartic interaction term in the Hamiltonian is replaced by a sum of terms of even degrees. In this subsection we restrict our attention to the case where only one such term is included, with degree q ≥ 4. The qualitative physics of the models m T q = 4 q = 6 q = 8 q = 10 q = 12  Table 1, along with the charge at the critical point Q The q-body interacting generalizations of the SYK model and closely related models have been considered in the literature (see e.g. [9,55,68]) because the large q limit can be solved analytically and captures interesting aspects of the physics of the models at finite q. However, as we have already mentioned in our comments on the Maldacena-Qi model in Sec. 4.2, the large q limit of the U(n) 2 model turns out to be rather trivial when defined in the standard way. In particular it is not able to capture the line of first order phase transitions or the associated critical point, which are the most interesting features of the model. Trying to understand better the large q behaviour in this model, as well as possible q-body generalizations of the U(n) model, could be an interesting direction for future research.
For our purposes, the main interest of considering the q-generalizations of the U(n) 2 model introduced in App. A is that we can study the properties of the critical point for different values of q. There are various natural observables one may consider; we choose to focus on the following quantities: i) The heat capacity as a function of temperature at m = m (q) c . We shall compute the right and left critical exponents α which govern the singular behaviour of C when T approaches the critical value T iii) The charge Q as a function of temperature at m = m (q) c , with right and left critical exponents q iv) The charge difference ∆Q between the two phases along their coexistence line, as a function of temperature. This quantity is of course defined only for T < T (q) c , and the associated critical exponent β (q) is defined by Note that we could also consider the same quantity as a function of mass, but since m (q) t (T ) is always smooth the associated critical exponent is also β (q) , a fact that we have verified numerically.
v) The entropy S as a function of mass at T = T vi) The charge Q as a function of mass at T = T (q) c , with right and left critical exponentsq To obtain all of these critical exponents we proceed in a way which is very similar to what we did before for the scaling dimension. For example, to extract q ± we can finally fit the Q (q) ± coefficients in (4.10) directly, again using an appropriate fitting window, as shown in Fig. 29b.
The example of this procedure illustrated in Fig. 29 for q = 4 is entirely typical, with other exponents and q-generalizations behaving analogously. The results for q = 4, 6, 8, 10, 12 are listed in Table 2.
Let us emphasize three important properties of these results.
First, the critical exponents that we find are not mean-field. In other words, the critical points are non-trivial and strongly coupled. This may appear as a surprise,    actions that we extremize, (3.86) and its q-body generalization (A.2). The crucial difference with the standard mean-field theory is that these effective actions depend on functions and not simply on a finite number of variables. We see that this functional version of mean-field theory yields a much richer theory of critical phenomena than Landau's, including the possibility for non-trivial critical exponents.
Second, the critical exponents are asymmetric, for instance α (q) etc. This is an unusual property [69]. In ordinary local quantum field theory, standard renormalization group arguments imply scaling forms near the critical points which ensure that the critical exponents must be symmetric (the constants of proportionality, like C (q) ± in (4.6), may however be asymmetric). This general result does not contradict our findings, since the models we are dealing with here have no notion of space or locality.
Third, the critical exponents are q-dependent, especially for lower values of q. We thus obtain an infinite series of inequivalent critical points. We shall argue in the next subsection that these points actually belong to a continuous line of inequivalent critical points.
Clearly, we are here just scratching the surface of a very non-trivial new phenomenon in large N quantum mechanical models. The critical points we are dealing with remain to a large extent very mysterious, and further study is needed to better understand their properties. In particular, we do not know if some version of the renormalization group formalism can be applied in the present context.

Multi-q models, domain walls and a line of critical points
When turning on several couplings λ q at the same time, one expects to find an even richer structure of critical and possibly multi-critical points in the corresponding phase diagrams. A detailed study is outside the scope of the present work, but we want to display here two interesting features.
RG flow and domain walls One way to visualize the renormalization group flow from the asymptotically free regime ω λ q to the conformal region ω λ q at low masses is to plot in log-log scale the imaginary part of the Matsubara coefficients Im G k as a function of the Matsubara frequencies ν k . This is illustrated in Fig. 30a for the standard case of the U(n) 2 model, where only λ 4 = 1 is turned on and we take m = 0 for simplicity. The structure of the RG flow is then clearly visible: we have at large frequencies a −1 slope characteristic of the asymptotic regime (3.18), G k ∼ i/ν k ; and a slope of 2∆ − 1 = −1/2 determined by (3.44) in the conformal region at low frequencies. The transition between the two regimes occurs at frequencies of order ω ∼ λ 4 = 1 and is very visible at the small but finite temperature T = 10 −4 used here. Now consider the case shown in Fig. 30b, where an additional coupling, chosen to be λ 8 = 200 for illustrative purposes, is turned on. There is a deep-UV region ν k λ 8 again governed by asymptotic freedom, as well as a deep-IR region ν k λ 4 governed by the conformal behaviour associated with the most relevant coupling λ 4 . Moreover, between them now lies an intermediate region λ 4 ν k λ 8 , which is approximately conformal with a slope 2∆ − 1 −0.75, corresponding to the dimension associated with the q = 8 interaction term, ∆ 1/8. Note that we are using a model with λ 4 and λ 8 turned on, instead of, for example, λ 4 and λ 6 . This is because the larger difference in the associated exponents ∆ makes the intermediate approximately conformal region easier to visualize when one works at small yet finite temperature and ratio λ 4 /λ 8 , a limitation inherent to the numerical approach. More generally, one may turn on many couplings λ q 1 , . . . , λ qr , organized hierarchically with q i−1 < q i and λ q i−1 λ q i . In the frequency regime λ q i−1 ω λ q i the model is then in an approximate conformal phase dominated by the coupling λ q i and characterized by an exponent ∆ = 1/q i . In a holographic picture, along the lines of [23], this behaviour should be mapped to interesting domain wall solutions interpolating between different AdS 2 spaces. To the best of our knowledge, these solutions have not yet been discussed in the literature.  A line of inequivalent critical points In multi-q models we expect that the phase diagram smoothly depends on the various couplings, so that we have in fact critical lines, surfaces or hyper-surfaces when considering models with two, three or more interaction terms, respectively.
One of the simplest cases to analyze corresponds to the model in which both λ 4 and λ 6 are turned on, with λ 6 = 1 setting the units. Its phase diagram is depicted in Fig. 31, projecting different values of λ 4 to the (m, T )-plane. We see that there is a critical line interpolating between the critical point of the q = 6 model when λ 4 = 0 and the critical point of the q = 4 model when λ 4 λ 6 . This is nicely confirmed by our analysis of the critical exponents. For example, in Fig. 32 we show the entropy exponents s ± , computed as explained for q ± in Fig. 29, as a function of λ 4 . We observe a continuous variation between the values s  Table 2), with the gap s − − s + smoothly widening as λ 4 becomes larger. This provides convincing evidence that we have a continuous line of inequivalent critical points, for which the coupling λ 4 is exactly marginal. Unlike the case of the U(n) 2 model that we have discussed in the previous Section, the thermodynamics of the U(n) model is governed by a unique solution to the Schwinger-Dyson equations that can be analytically continued for all values of m and T without undergoing monodromies. In particular, we do not find first order phase transitions or non-trivial critical points at finite temperature in the U(n) model. All the interesting physics occurs at zero, or near zero, temperature. As discussed in Sec. 4.1, this important difference with the U(n) 2 model is made possible by the fact that the zero-temperature Euclidean two-point function gets non-trivial corrections even at relatively large masses.

Phase diagram of the U(n) model
More precisely, at large enough mass the system is weakly coupled. The groundstate must then coincide with the Fock vacuum |0 defined by the conditions ψ µ |0 = 0, with a mass gap m eff of order m. The spectral function (3.9) then has a support only for ω ≥ m eff > 0 at zero temperature, and the spectral decomposition formula (3.12) implies that In Fig. 33 we plot the mass gap m eff as a function of the mass parameter m, evaluated numerically by matching finite-temperature solutions with the exponential behaviour (4.16) and then extrapolating m eff to T → 0. The essential feature to note here is that the gap closes for a finite value m = m * = 0.44 1 .
A complementary point of view is provided by the constraint (3.49). In Fig. 34 we plot Re Σ 1/2 as a function of the mass at T = 0.01 in natural units. We observe here a clear crossover between the regime at low masses approximately satisfying (3.49), and the regime at large masses for whichΣ(0) 0. The turning point is at m = m * = 0.44 1 , perfectly consistent with the appearance of the effective mass discussed above.
The point m = m * is a new kind of large N , strongly-coupled quantum critical point. For m < m * the zero-temperature Euclidean two-point function behaves as (3.55), with asymmetric exponents ∆ + = ∆ − . This feature was already discussed in detail in Sec. 3.3.2, and is in sharp contrast with the IR behaviour of the two-point function in the U(n) 2 model. As pointed out in Sec. 3.1.5, the SL(2, R) subgroup of the low energy reparameterization symmetry (3.38)-(3.39), including the scale transformations τ → ατ for α > 0, is now spontaneously broken, even though the phase is gapless. It is thus a phenomenon of spontaneous scale symmetry breaking without the generation of a mass gap.
Plots of the energy, entropy and heat capacity divided by T are provided in Fig. 35. At zero temperature these quantities are all vanishing for m > m * , but become nonzero for m < m * , very much like in the U(n) 2 model (compare with Fig. 27). The low-temperature properties of the U(n) 2 and U(n) models are thus quite similar as far as the thermodynamical potentials are concerned. For example, they both have Re Σ 1/2 m Figure 34: Plot of Re Σ 1/2 Σ (0) as a function of the mass (black), for the U(n) model at T = 0.01 in natural units. We observe a clear crossover from a gapless phase, which must haveΣ(0) = m at zero temperature (red), as in (3.49), to a gapped phase, for whichΣ(0) 0 at very large masses. The transition occurs at m = m * = 0.44 1 , consistently with the appearance of a non-vanishing effective mass, see Fig. 33.
non-vanishing zero-temperature entropy and a heat capacity vanishing linearly with the temperature in the gapless phase m < m * . A basic difference with the U(n) 2 model is that the entropy is now continuous. Its derivative is discontinuous, however, and this translates into a finite jump for C/T at m = m * .
We conclude this Section with a few comments on some subtleties, both conceptual and numerical, associated with the new IR behaviour found in the U(n) model for 0 < m < m * .
We have argued that the deep IR behaviour changes discontinuously from the standard SYK behaviour at m = 0, with symmetric exponents ∆ + = ∆ − = 1/4, to the new asymmetric behaviour with ∆ + = ∆ − for any m > 0, even if arbitrarily small. In this sense, the mass operator is relevant at m = 0 (and marginal for all 0 < m < m * ), even though no mass gap is created: as soon as it is turned on, the IR fixed point of the model jumps from one gapless solution to another very different gapless solution. The fact that the jump must be discontinuous can be understood e.g. as a consequence of the formulas (3.54) and (3.57) for b 4 + , which show that this coefficient is discontinuous when one goes from ∆ + = ∆ − to ∆ + < ∆ − , even if the difference between ∆ + and ∆ − is very small. This discontinuity at m = 0 contrasts with the fact that all the thermodynamical quantities are found to be perfectly continuous when m → 0, including at T = 0, as shown in Fig. 35. At m = 0, they match with the results for the U(n) 2 model, which is essentially just SYK. But there is no contradiction here; the discontinuous change in the IR behaviour occurs only for a range of frequencies |ω| m that shrinks to zero size when m → 0. At very small m |ξ|, the model thus has an interesting RG flow: from the asymptotically free regime at frequencies ω |ξ|, it flows to an approximate conformal phaseà la SYK with ∆ = 1/4 in the regime m ω |ξ| before eventually flowing to the non-standard deep IR behaviour with ∆ + < ∆ − when ω m.
The fact that the correct IR behaviour can emerge only for ω m implies that the numerical characterization of the deep IR regime at small mass is particularly difficult. To do this, one should probe temperatures satisfying βm 1 as well as β|ξ| 1, but the former condition is much stronger than the latter in the regime m |ξ|.
To illustrate this difficulty, we plot in Fig. 36  values of m, evaluated numerically by following exactly the same procedure as in Fig.  15 of Sec. 3.3.2. We see that ∆ + stays reasonably close to the value ∆ + = 1/4, as required by (3.57), with small variations that may be attributed to numerical artifacts surviving the T → 0 extrapolation. This is not surprising, since the value ∆ + = 1/4 coincides for both the standard SYK conformal phase and the new non-standard phase with ∆ + < ∆ − . A precise determination of ∆ + is thus expected as soon as β|ξ| 1, even if m |ξ|. The situation is very different for ∆ − . The numerical values we find exhibit a marked mass dependence, starting from ∆ − = 0.250 1 at m = 0 and increasing up to ∆ − = 0.446 2 at m = 0.29.
The fact that the numerical result for ∆ − varies smoothly from its m = 0 value is an automatic consequence of the discussion above: at extremely small masses the lowest temperature we can study numerically will always be in the near-SYK regime 1/|ξ| β 1/m and the genuine deep IR regime β 1/m is not probed. It is therefore not possible to decide, based on the numerical analysis alone, whether the actual value of ∆ − in the phase ∆ + = ∆ − does depend continuously on m or not. Two possibilities thus remain open: i) The deep IR value of ∆ − depends smoothly on m and our numerical analysis can actually provide a reasonable approximation to the actual value, possibly even for small masses.
ii) The deep IR value of ∆ − > ∆ + is fixed and does not depend on m, but our numerical analysis cannot reliably determine it for small m. The mass dependence of ∆ − in Fig. 36 may then be interpreted as mimicking a RG flow, since the deep IR is probed better and better when m increases (as explained below, this is no longer true when the mass is too close to m * ).
Note that, in either case, the fact that ∆ + < ∆ − when m > 0 is unambiguously established. As will be explained elsewhere [70], the possibility ii) seems to be strongly favoured.
Finally, let us mention that obtaining numerical solutions at very low temperatures becomes increasingly harder as the mass approaches the critical value m * = 0.44 1 , which explains why the highest mass appearing in Fig. 36 is m = 0.29. Indeed, the solutions near the quantum critical point are quite subtle. They must be very near the low entropy, perturbative-like solutions above some effective mass scale that goes to zero when m → m * from below, but at the same time should reproduce the non-standard IR behaviour with ∆ + < ∆ − below this small effective mass scale.  (3.31). They can be solved numerically using a method outlined in Sec. 5.1.3. Note that our results generalize and are consistent with known formulas, presented for example in [9] for the SYK model; the only true novelty here may be the rigorous justification of the derivations, especially in Sec. 5.1.2, and the pedagogical presentation.
Real time two-point functions can be used to compute quasi-normal frequencies.
The results are presented in Sec. 5.1.4. We find a very interesting behaviour near the finite temperature critical point in the U(n) 2 model. Real time two-point functions will also serve as a crucial building block in the analysis of the four-point functions in Sec. 5.2, allowing the calculation of Lyapunov exponents in the following subsections.

General definitions and properties
For any operator O and complex "time" u we define O(u) = e Hu Oe −Hu , (5.1) denoting τ u = Re u , t u = Im u .
If u = τ u is real (5.1) is equivalent to the Euclidean time evolution (3.2). If u = it u is purely imaginary, on the other hand, (5.1) yields the standard real time Heisenberg evolution. We can then define a generalized two-point function for complex time arguments u and v satisfying 0 < |τ u − τ v | < β, where T indicates an ordering with respect to the real parts of u and v, i.e. with respect to the Euclidean time. This is an important point: it is the Euclidean time that always dictates the ordering of operators in a correlator. Clearly G(u, v) depends on u − v only, and we shall thus also note G(u, v) = G(u − v). The KMS condition can be used to extend G(u) to all values of τ u . It is sometimes convenient to decompose the generalized two-point function in terms of Wightman correlators, which are defined by 23 At tree level, generalizing (3.7). The Euclidean spectral decomposition (3.12) is generalized to 8) 23 We stick to the matrix model notation from now on.
where the spectral function ρ is defined in (3.9). We can also use the Fourier decomposition It is important to note that the correct analytic continuation from Euclidean time to complex time cannot be obtained by replacing τ → τ + it in the Euclidean Fourier expansion (3.4). The analytic continuation to complex time is simple in correlators for which the operator ordering is fixed, because the complex time operator evolution (5.1) is a simple analytic continuation of the Euclidean time evolution (3.2), but is non-trivial in time-ordered correlators.
Purely real time correlators are obtained from the complex time correlators by letting τ go to zero, the ordering of operators being fixed by the way the limit is taken, i.e. τ → 0 + or τ → 0 − . Eq. (5.8) implies that all the real time two-point functions, for any ordering of the operators, can be expressed in terms of the spectral function ρ. For instance, It is also convenient to introduce the retarded and advanced two-point functions for which From the definition of the resolvent in (3.22), we see that χ r (ω) = −R(ω + i0 + ) and χ a (ω) = −R(ω − i0 + ) . (5.17) The retarded and advanced Fourier transforms χ r and χ a can thus be analytically continued to holomorphic functions on the upper and lower half-planes respectively, an important and well-known property. Of course, they are not independent since the reality of the spectral function implies that Real time two-point functions are in this way completely encoded in the spectral function ρ, whereas the Euclidean two-point function is encoded in the Fourier-Matsubara coefficients G k . The real time and Euclidean time data thus look very different: frequencies in real time are always continuous, whereas their Euclidean counterparts are discrete. There is a powerful theorem due to Carlson [71] that shows that these two different-looking sets of data are nevertheless strictly equivalent: the real time data determines uniquely the Euclidean time data and vice versa. That the real time data determines the Euclidean time data is easy to understand: from the spectral function ρ one can get the G k 's from (3.22) and (3.23). The converse is much more subtle. The knowledge of the Matsubara-Fourier coefficients fixes the values (3.23) of the resolvent on a discrete set of points z = iν k on the imaginary axis. Moreover, we know that the resolvent is a holomorphic function for Im z > 0 and Im z < 0 and has the simple asymptotic behaviour (3.24). Carlson's theorem then implies that R(z) is determined uniquely for all complex values of z.
The theorem itself does not provide an explicit procedure to reconstruct ρ from the G k 's, but such a procedure was discovered in [72] and generalized in [57]. The explicit formulas in [57] actually allow to reconstruct the resolvent from the knowledge of the Euclidean data above an arbitrary cut-off k 0 in frequency, i.e. the set of Matsubara-Fourier coefficients {G k , k ≥ k 0 > 0} . A startling consequence of this fact is the existence in any quantum mechanical theory of a rigorous version of the Renormalization Group, as a consequence solely of analyticity properties [57]. 24 The explicit reconstruction formulas are thus very interesting conceptually, but in practice their use is limited by the fact that a reliable reconstruction requires the knowledge of the G k 's to an enormous precision. They are therefore impractical except in simple cases where ρ is known to be very smooth, without sharp features. Since we do not want to make any such a priori assumptions, we prefer to base our analysis on the real time version of the Schwinger-Dyson equations derived in the next subsection. Interestingly, the derivation of these equations that we now present is going to use Carlson's theorem in a central way.

Melonic real time two-point functions
The standard general approach to real time correlation functions is based on the Schwinger-Keldysh formalism. This formalism is a direct consequence of the path integral formulation, suitably adapted to cases where operators are inserted at arbitrary complex times, as in (5.1). One introduces an oriented closed Keldysh contour γ on the complex time cylinder where u and u + β are identified. The contour starts at u = 0, ends at u = β and goes through all the operator insertion points, say O 1 (u 1 ), . . . , O r (u r ). The correlator associated with such a contour is where the symbol T γ orders the operators according to their position along the contour. For instance, to compute the two-point function G(u 0 ), Eq. (5.3), we can use the contour depicted in Fig. 37, with an insertion of ψ † µ at u = 0 and an insertion of ψ µ at u = u 0 . The standard Euclidean formalism corresponds to the special case where all the insertions are on the Euclidean time axis Im u = 0, so that we can then use the contour Im u = 0. In practice, it is convenient to use contours along which Re u is monotonically increasing, so that the contour ordering coincides with the usual Euclidean time ordering. This is why the real time piece of γ in Fig. 37 was chosen to have a triangular shape. at u = u 0 are represented by bullets along the contour. The width of the triangleshaped piece of the contour may be taken to zero. Note that the contour may be chosen such that the insertion point at u = u 0 sits at the top of the triangle, but this is not required.
The advantage of this approach is that the Feynman rules take a familiar form. The only differences with the Euclidean case is that integrals over τ are replaced by contour integrals along γ, and propagators are associated with the complex time treelevel two-point function (5.7). The arguments leading to the Euclidean Schwinger-Dyson equations (3.30) and (3.31), in their path integral or diagrammatic versions for the disordered and matrix models respectively, can be repeated without change in the present, more general context. The self-energy for the contour of Fig. 37 is defined by the equation Traditionally, one proceeds by using the decomposition (5.5) to rewrite (5.20) in the form of the so-called Kadanov-Baym equations; see e.g. [73] for a pedagogical introduction. The Kadanov-Baym equations can then be integrated numerically using (5.21) or (5.22). 25 We are not going to use this rather tedious route here. Instead, we present a simple argument based on Carlson's theorem, leading to a very convenient expression of the solution which entirely avoids having to deal with the details of the Kadanov-Baym equations.
Let us start by defining Using (3.22) and denoting x = Re z and y = Im z, we get From the positivity of the spectral function ρ, which is manifest from its definition (3.9), we deduce that Im R(z) = 0 away from the real axis, and thus also R(z) = 0 there. Since R is holomorphic and has no zeros for Im z > 0 or Im z < 0, Eq. (5.25) implies that R Σ is also holomorphic for Im z > 0 and Im z < 0. Moreover, the expansion (3.24) shows that All the conditions for the applicability of Carlson's theorem to the function R Σ are thus met: if we are able to find some functionR Σ (z), holomorphic for Im z > 0 and Im z < 0, behaving as O(1/z) when z → ∞ and such thatR Σ and R Σ match for all z = iν k , then the theorem states thatR Σ = R Σ for all z: We can combine the equations (5.21) or (5.22) with the spectral decomposition (5.8) to guess a formula for R Σ (z). We shall then check that the above criteria are satisfied.
In the case of the U(n) 2 model, we consider The convergence of the integral is trivial to check. The functionR Σ so defined is manifestly holomorphic for Im z > 0 and Im z < 0 and is O(1/z) at large z. Moreover, if we evaluate Σ k by taking the Fourier transform of (5.21), in which we use the spectral decomposition formula (5.8) for G and the explicit form (5.7) for G (0) , we immediately find that Σ k = −R Σ (iν k ) is satisfied. From (5.24) we conclude that R Σ (iν k ) = R Σ (iν k ) and thus from (5.27)R Σ = R Σ .
Exactly the same reasoning can be applied to the U(n) model starting from and then showing thatR Σ (iν k ) = R Σ (iν k ) and thusR Σ = R Σ .
We should stress that it is crucial to understand that the propertyR Σ (iν k ) = R Σ (iν k ) = −Σ k is not sufficient in itself to concludeR Σ (z) = R Σ (z). For example, a naive way to do the analytic continuation from Euclidean time to real time in the U(n) 2 model is to start from the Euclidean equation (3.30) written in terms of the Matsubara-Fourier coefficients, This might induce us to propose Although this expression satisfies by construction R Σ (iν k ) = −Σ k , it does not have the required analytic properties; indeed, the right-hand side of (5.31) is holomorphic only within bands ni < Im z < (n + 1)i where n ∈ Z, but not on the whole half-planes Im z > 0 and Im z < 0 as required.
For the purpose of improving the efficiency of the numerical algorithm presented in the next subsection, it is useful to perform explicitly one of the integrals over the frequencies in (5.28) and (5.29). For example, in the U(n) 2 model we write We can then use (3.25) and the well-behaved properties of the integrand at infinity to get where C + and C − are semi-infinite rectangles on the upper and lower half complex w-plane respectively, see Fig. 38. Picking up the poles at w = iν k for all k ∈ Z + 1 2 , as well as at w = Ω − z, yields Overall, we thus obtain for the U(n) 2 model In the U(n) model, similar manipulations yield Equations ( Let us conclude this derivation with the following remark. The Kadanov-Baym equations obtained from (5.20) must of course be satisfied by our solution, and we have checked this explicitly. This is an interesting albeit somewhat tedious exercise in contour integration. One must consider several cases depending on the position of u and u along the Keldysh contour γ. It is convenient to use double Fourier decompositions of the form (5.9) for both G(u) and Σ(u). All the required identities follow from the analytic properties of R and R Σ discussed above.

Numerical solution of real time equations
As already mentioned, all the real time two-point functions can be straightforwardly derived from the spectral function ρ(ω). The basic formula relating ρ to the resolvent R is (3.25), and noting that for any real we can rewrite this equation as We will therefore focus on the computation of the resolvent through equations (5.23) and (5.36) or (5.37). In practice, we evaluate R(ω + i ) for a small positive , and assume throughout The real time Schwinger-Dyson equations can be solved numerically using an iterative algorithm that is very similar to the one applied in the Euclidean case, described in detail in Sec (ω + i ) by plugging these in the right-hand side of (5.36) or (5.37), according to the model under study. To complete the iteration, we define the resolvent at step i + 1 to be , (5.41) where α is an a priori arbitrary real weighting factor.
Because equations (5.36) and (5.37) involve integrals over ω, in practice we need to choose a finite number of discrete frequencies over which the resolvent is evaluated, This defines a unidimensional grid on a finite interval [ω min , ω max ] = [δω p i , δω p f ], outside which we use asymptotic freedom to identify the solution with the tree-level result (5.44), thus transforming all integrals into finite sums. The discrete step δω is chosen consistently with the parameter > 0, and sets an intrinsic upper limit on the resolution with which ρ can be reconstructed from (5.40).
Note that the Euclidean two-point function also appears on the right-hand side of equations (5.36) and (5.37) through its Matsubara-Fourier coefficients G k . We will therefore assume that a solution for these is already known, possibly by having run the Euclidean algorithm beforehand, so that they are simply constant coefficients in this context.
The iteration of the algorithm is repeated until the equations are solved up to a prescribed tolerance ε, i.e.
As in the Euclidean case, we will take ε = 10 −12 (in natural units λ = |ξ| = 1), ensuring that the solution is found with a very good accuracy while staying safely above machine-precision issues.
Before moving on, let us make a few further comments on the algorithm described above. First, the choice of initial approximation is very important for its convergence. We may use the tree-level result corresponding to R (0) 44) or any previously found solution deemed to be close to the actual solution we are seeking. As in the Euclidean case, this allows us to map the parameter space by taking small steps in either the m or T directions to "trace" solutions.
Another common element with the Euclidean algorithm is the presence of an arbitrary weighting factor α. This comes about for the same reasons discussed in Sec. 3.3, having a similar impact on convergence and performance as was seen there, and again we find that α = 0.75 works well for most cases.
A new ingredient in the real time algorithm is the arbitrary frequency-domain discretization (5.42), which was naturally given by the Matsubara frequencies in the Euclidean setting. While as already mentioned δω is related to the resolution of the resulting spectral function, its support is determined by p i and p f . Taking Eq. (3.21) into account, these parameters need to be chosen such that the interval [ω min , ω max ] is approximately centered around ω = m, and extends to either side at least up to the regime where asymptotic freedom takes over.
More generally, the choice (5.42) used here has the advantage that it is additive and symmetric about the origin, so that it satisfies the convenient properties ω p 1 + ω p 2 = ω p 1 +p 2 and − ω p = ω −p . (5.45) These properties simplify the practical implementation of the algorithm and improve its performance by completely avoiding costly interpolations. Furthermore, they allow for the efficient pre-calculation of the Matsubara-Fourier coefficient sums in (5.36) and (5.37), which after truncation of the Euclidean data can be done in O(P ) memory and O(P × k UV ) time, instead of the naive O(P 3 ) and O(P 3 × k UV ), respectively. It is for this reason that we chose to use the formulas (5.36) and (5.37), with two integrals over frequencies, instead of (5.28) and (5.29), with three integrals over frequencies. 26 There is no guarantee that the algorithm described here terminates, but if it does so in I α (ε) steps then the total running time is O(I α (ε) × P 3 + P × k UV ). Additional pre-processing in O(P 2 ) time and memory is useful in practice to reduce constant factors and the incidence of numerical errors, but this is not relevant asymptotically. For our C++ implementation running on a desktop computer one may then take P = O(10 3 ), which makes the pre-processing subleading for typical values of k UV = O(10 5 ). Depending on the model and values of the mass, this is enough to reach temperatures of order T = O(10 −2 ) in natural units, but quickly becomes insufficient below that point.
To illustrate the results obtained with this algorithm, we show in Fig. 39 examples of the spectral function computed in the U(n) 2 model, for various masses and temperatures. At m = 0 it is an even function as a consequence of the particle-hole symmetry, and coincides with the spectral function of the U(n) model. These properties are lost when m > 0. At small masses and temperatures, ρ(ω) approaches the conformal form (3.53), whereas for large masses it converges to the perturbative result δ(ω − m).
In order to check the precision of our numerical methods, we show in Fig. 40a the values of the first two moments of the spectral function, ρ 0 and ρ 1 defined in Eq. (3.17), which are to be compared with the exact results (3.20) and (3.21), i.e. ρ 0 = 1 and ρ 1 = m. Another check is provided in Fig. 40b, where we plot the energy computed from the real time data via Eq. (3.103), and compare with the results obtained from the Euclidean formula (3.102). Similar results are obtained for the U(n) 2 model at temperatures other than T = T c , used here for illustrative purposes, as well as for other thermodynamic quantities such as the fermion number (computed from Euclidean or real time data through (3.96) and (3.97), respectively).
While Fig. 39 and 40 correspond to the U(n) 2 model, very similar results are also obtained for the U(n) model. The main difference in this case resides in the fact that the spectral function has a more involved structure for m > 0, presenting multiple peaks instead of just one, see Fig. 41. While this is not intrinsically problematic, these peaks become sharper as the temperature is lowered, thus increasing our resolution requirements and therefore raising the lower bound on the temperatures we can probe numerically in a reasonable amount of time.  Fig. 12 for the solution at T = 0.01). Plot (c) corresponds to the HO-like phase, where the low-temperature limit ρ(ω) → δ(ω − m) is already evident for T = 0.03 (vertical axis is truncated for clarity). All solutions were obtained with P = 1001, ∆ω = δω(p f − p i ) = 10 and = 10 −5 , except for the lowest temperature in (c) for which ∆ω = 0.1 and = 10 −7 .

Application: quasi-normal frequencies
Quasi-normal frequencies ν (p) qn , p = 0, 1, 2, . . ., govern the way a system responds to a small perturbation. Linear response theory relates these frequencies to the long time behaviour of the retarded two-point function G r , with frequency ν (p) qn contributing a term proportional to exp(−iν (p) qn t). Equivalently, they correspond to poles of χ r on the lower-half complex frequency plane. 27 A non-zero negative imaginary part in ν (p) qn 27 Or to poles of χ a on the upper-half plane, see (5.18). indicates an exponential decay at long times of the associated mode, with time scale In the context of black holes, having Im ν (p) qn < 0 is a landmark feature of horizon physics, signalling the absorption of any small external perturbation by the black hole and its subsequent return to equilibrium.
For example, for the U(n) 2 model at low temperatures in the SYK-like phase, we can use (3.50) and (3.43) to write down the retarded function in the regime βλ 1, As expected, the large t asymptotic expansion of this function is a sum of decreasing exponentials of the form exp(−iν Amongst all the quasi-normal frequencies ν (q) qn , the one having the largest time scale (5.46) dominates at long times. It is called the leading quasi-normal frequency and will be denoted simply by ν qn , The associated time scale t qn = − 1 Im ν qn (5.50) is often referred to in the literature as the thermalization or dissipation time scale.
In Fig. 42, we plot − Im ν qn and Re ν qn as functions of the temperature, for various values of the mass in the U(n) 2 model. As explained in Sec. 4, for masses m < m * = 0.213 2 the IR physics is dominated by the SYK-like solution; our numerical results are then consistent with the asymptotic values (5.48) for p = 0. For masses m > m * , on the other hand, we find that the behaviour of − Im ν qn is o(T ).
The situation in the U(n) model is illustrated in Fig. 43. For low enough masses, in the gapless phase, the low-temperature behaviour appears to be linear for both the real and imaginary parts of ν qn , but the slope is seemingly dependent of the mass. Note that we have no analytical results to contrast with this observation. For larger masses, in the gapped phase, numerics do not allow us to reliably probe the very low temperature regime.
Just as with the thermodynamic quantities studied in Sec. 4.2.3, the behaviour of the leading quasi-normal frequency near the critical point of the U(n) 2 model turns out to be very interesting. In Fig. 44 we plot − Im ν qn as a function of the mass at T = T c and as a function of the temperature at m = m c . The results are consistent with a critical power-law behaviour,

Real time four-point functions
This subsection is devoted to the computation of real time four-point functions, in particular in the out-of-time-ordered (OTOC) case which is relevant to probe chaos [6,34,75,76]. We review basic properties in Sec. 5 In particular, we find that the chaos bound of [34] is saturated in the SYK-like phase of the U(n) 2 model, but find evidence that this is not so in the U(n) model except at strictly m = 0. We also find that the Lyapunov exponent has a singular behaviour near the critical point of the U(n) 2 model, governed by non-trivial critical exponents.

General definitions and properties
A four-point function for complex time arguments u i satisfying 0 < τ u i < β, generalizing (3.120), can be defined by where as usual T indicates the ordering with respect to the Euclidean times τ u i = Re u i . This definition can be extended to all values of τ u i using the KMS condition φ(u 1 + n 1 β, u 2 + n 2 β, u 3 + n 3 β, u 4 + n 4 β) = (−1) n 1 +n 2 +n 3 +n 4 φ(u 1 , u 2 , u 3 , u 4 ) , n i ∈ Z . (5.54) This four-point function satisfies the properties generalizing the similar conditions (3.122), (3.123) and (3.124) in the Euclidean setting. One may decompose where the sum is over all the elements of the group S 4 of permutations of four elements, orders the Euclidean times and φ σ is one of the 4! = 24 possible four-point correlators with a fixed given ordering of the operators. Note that, by construction, the Euclidean times appearing as the arguments of φ σ are always ordered according to the permutation σ, with τ u σ(1) > τ u σ(2) > τ u σ(3) > τ u σ(4) and 0 < τ u σ(1) − τ u σ(4) < β. The real times (imaginary parts of the u i ), on the other hand, are unconstrained. Any purely real time four-point function can be obtained by taking the Euclidean times to zero while keeping their ordering fixed. Spectral decomposition formulas read with spectral functions that can be expressed in terms of sums over the eigenstates |p and eigenvalues E p of the Hamiltonian, The coefficients α σ pqrs depend on the matrix elements of the operators and on the permutation σ in a natural way. If we use the explicit notation Correlators with the same cyclic ordering of the operators are related by the KMS conditions (5.54). These conditions follow, for instance, from the cyclicity property of the coefficients α σ pqrs , In terms of the spectral functions (5.61), they read The KMS conditions thus reduce from 4! = 24 to 3! = 6 the number of independent spectral functions or correlators. Moreover, in our particular case one may check that the additional symmetry (5.56) leaves only four independent spectral functions or correlators, The other two a priori independent correlators are given by (5.72) The above extends the discussion of Sec. 5.1.1 for the two-point functions, for which there is only one independent spectral function. It can be straightforwardly generalized to arbitrary p-point functions, with generically (p − 1)! independent spectral functions, see e.g. [77] for a thorough discussion.

Melonic real time four-point functions
By time translation invariance, we can restrict our attention to φ(u 1 , u 2 , u 3 , u 4 = 0). For the purpose of studying chaos, along the lines of [6,34], we can further limit ourselves to a real u 3 = τ 3 .

(5.73)
We also fix the real parts of u 1 and u 2 to be τ + and τ − (or τ − and τ + ) such that 0 < τ − < τ 3 < τ + < β , (5.74) and we note t 1 = Im u 1 , t 2 = Im u 2 . (5.75) The associated Keldysh contour Γ for these choices of insertion points has two triangleshaped pieces and is depicted in Fig. 45. The height of the triangles must be at least max(t 1 , t 2 ) but may conveniently be taken to be infinite. As already emphasized in Sec. 5.1.2, the advantage of the Keldysh contour formalism is that the Feynman rules and diagrammatics take the familiar form, with only simple and obvious modifications with respect to the standard Euclidean case. In particular, the arguments presented in Sec. 3.4 go through. The leading Feynman graphs are still the ladders depicted in Fig. 17, 18 and 19, but now the positions of the vertices have to be integrated over the Keldysh contour Γ instead of the Euclidean interval [0, β]. If we decompose with integrals taken along Γ, (5.77) is solved by resumming a geometric series or equivalently inverting the matrix I − K, The points u 1 and u 2 in (5.77) may be chosen to be anywhere on the contour Γ. The real time OTOCs correspond to insertions on the two triangles, i.e.
which yield (5.82) These correlators are decomposed, following (5.76), as The two connected OTOCs (5.85) and (5.86) can be inserted on the left-hand side of Eq. (5.77) by setting u 1 and u 2 to the appropriate values (5.80). However, it is important to realize that several other types of correlators then necessarily appear on the right-hand side of (5.77) when one performs the integrals over Γ. This includes, for example, correlators like 87) when u and u are both on the triangle at Re u = Re u = τ + , similar correlators when u and u are both on the triangle at Re u = Re u = τ − , and also correlators with only one real time turned on or purely Euclidean correlators when both u and u are on the Euclidean segments of Γ.
It is thus impossible to write an exact closed-form system of integral equations similar to (5.77) involving the OTOCs alone. 28 However, we shall explain in Sec. 5.2.4 that such an approximate system can be used if we are only concerned with evaluating the asymptotic, long time behaviour. It is this behaviour which is relevant to compute the chaos Lyapunov exponents, which we discuss in the following.

On Lyapunov exponents
Let us start by summarizing the standard lore on the late real time behaviour of fourpoint functions such as φ(u 1 , u 2 , τ 3 , 0) in models like the ones we are studying. Recall that t 1 and t 2 are defined in (5.75); we want to discuss the regime where t = t 1 = t 2 is very large.
We assume that the number N of degrees of freedom is large but finite. This is important, because time scales of the order of ln N play a role in the discussion. Much larger time scales, like the Poincaré recurrence time scale, which are exponentially large with N , are irrelevant for our purposes and can be considered to be infinite.
Depending on the ordering of the operators in the correlators, i.e. on the positions of u 1 and u 2 on the Keldysh contour of Fig. 45, two qualitatively distinct behaviours are expected: i) If u 1 and u 2 are on the same triangle-shaped piece of Γ, then the correlator will factorize and approach a constant value as soon as t = t 1 = t 2 is much larger than the dissipation time scale defined in (5.50), The function G + or G − appears on the right-hand side of (5.88) depending on the relative position of u 1 and u 2 along Γ. Note that the result (5.88) is consistent with the large N expansion, since the right-hand side of this equation matches with the first term on the right-hand side of (5.76). In other words, the long time and large N limits commute for this class of correlators. 29 ii) If u 1 and u 2 are on distinct triangle-shaped pieces of Γ, which is the case corresponding to the OTOCs (5.81) and (5.82), it is expected that the correlator will go to zero when t = t 1 = t 2 → ∞, Intuitively, this is because the insertions at u 1 and u 2 are infinitely far apart on the Keldysh contour Γ in the limit t → ∞. The result (5.89) can actually be used as a good criterion to detect chaos in a strongly coupled quantum system [34]. Indeed, (5.89) is inconsistent with the first term on the right-hand side of (5.76), which predicts (5.88) at N → ∞. In other words, the long time and large N limits do not commute for the OTOCs. The F/N term, which naively is a small correction at large N , must eventually become of order one to cancel out the constant GG term. More precisely, it is expected that as soon as t = t 1 = t 2 t qn , the exponential growth of the "small" correction F/N signalling the onset of chaos. The constant η L is called the Lyapunov exponent. Eventually, the OTOCs should vanish on a time scale, called the scrambling time scale, that can be estimated by setting F/N ∼ 1 in (5.90), A few clarifying remarks may be made at this stage. First, in the strict N → ∞ limit studied in the present paper the scrambling time t scr ∼ ln N is not seen and the vanishing of the OTOCs cannot be rigorously derived. This is a non-perturbative effect in the large N expansion.
On the other hand, the behaviour (5.90) can in principle be derived by solving the exact equation (5.77). In practice, this can be done only in the low-temperature limit and in the conformal phases of the models. The method, presented in [9] for SYK (see also [31,[78][79][80]), does not use (5.77) explicitly. Instead, one first works in the Euclidean setting and exploits the conformal symmetry to explicitly invert the operator I − K in (3.127), and then performs the appropriate analytic continuations to real time on the exact solution. The result is the famous formula for the Lyapunov exponent in the zero-temperature limit, which matches with the expectations from black hole physics and holography [76]. The universality of the result (5.92) can be understood as being the consequence of the small explicit breaking of the low energy reparameterization invariance (3.38)-(3.39) by UV effects, which yields a universal Schwarzian action for the reparameterization zero modes from which the leading behaviour (5.90) can be obtained.
At non-zero temperature one can in principle derive the behaviour (5.90) by solving (5.77) numerically. We first compute the kernel by using the solution for the two-point functions presented in Sec. 5.1, and then invert the operator I − K to get F from (5.79). This procedure is in principle straightforward and allows the computation of the OTOCs from the very short time scales 0 ≤ t 1 , t 2 t qn to the chaos regime t 1 , t 2 t qn . However, the size of the matrix I − K that one needs to invert to obtain reliable results when t 1 , t 2 t qn is very large. To the best of our knowledge, this has never been done and we let this interesting challenge for the future. Now, if the detailed time evolution of the OTOCs is not required, but one is only interested in the regime t t qn , the exact equation (5.77) can be simplified. This is sufficient to extract the Lyapunov exponents, and is explained in the next subsection. The drawback of this approach is that it relies on the standard expectations concerning the behaviour of the real time correlators reviewed above, discarding in a somewhat heuristic way many terms on the right-hand side of (5.77). It is therefore not a first principles derivation strictly proving that this standard lore is indeed correct. However, the results obtained in this way are fully consistent with the original hypothesis and we believe that they are very reliable, see Sec. This bound is supposed to be valid in a large class of models for which a large hierarchy exists between the dissipation time scale and the scrambling time scale, t scr t d [34]. It is saturated by black holes and, as Eq. (5.92) shows, by the SYK model and many generalizations thereof, including the U(n) 2 model in the standard conformal phase (as emphasized in Sec. 4, the standard conformal phase is not realized in the U(n) model; this will be further discussed in Sec. 5.2.5).
The OTOCs we are considering are labelled by an arbitrary triplet (τ − , τ 3 , τ + ) constrained by (5.74). The most natural choice, which yields physical, real time finite temperature correlation functions, corresponds to taking the limits We will refer to these as the "standard" correlators, which can be interpreted as particular expectation values in a model for which we take two copies of the original Hilbert space and in a state which is a deformed version of the standard thermofield double state [34].
It is not obvious whether the bound (5.93) is valid independently of the choice of triplet (τ − , τ 3 , τ + ), even though we believe it is very likely to be so in melonic models. We shall compute Lyapunov exponents for different choices of triplets (τ − , τ 3 , τ + ) in Sec. 5.2.5 and our results support this claim. However, it is important to stress that the proof presented in [34] does not apply to the general case of arbitrary (τ − , τ 3 , τ + ) but only to the special choice τ − = β/4, τ 3 = β/2, τ + = 3β/4. 30 Indeed, the reasoning crucially depends on the possibility to analytically continue the OTOC from the real time t = t 1 = t 2 to a strip −β/4 < Im t < β/4 in the complex t-plane. This is indeed possible when τ − = β/4, τ 3 = β/2, τ + = 3β/4, but, unfortunately, no such analytic continuation exists in the physical case. In general, an analytic continuation is possible in a strip − min(τ 3 − τ − , β − τ + ) < Im t < min(τ − , τ + − τ 3 ), whose width can never exceed, but is in general strictly less than, β/2; the argument of [34] may then be used to derive an upper bound which is less strict than (5.93).

Computing Lyapunov exponents
Lyapunov exponents are determined by analyzing the long time exponential growth of the connected OTOCs (5.85) and (5.86), and comparing with (5.90). We start by writing (5.77) for the relevant values (5.80) of u 1 and u 2 . Using the complex time version of (3.128), we get The idea is that we can simplify the right-hand sides of these equations, because we are interested in the regime in which the OTOCs are exponentially large. We thus keep only the terms that can develop such an exponential growth in the integrals and neglect all the others. From the discussion in Sec. 5.2.3, we expect that these terms correspond to the integration regions u = τ + + it and u = τ − + it or u = τ − + it and u = τ + +it , for which F(u, u , τ 3 , 0) coincides with the OTOCs (5.85), (5.86). We are left with integrals running along the two sides of the triangles, which are non-trivial because the kernel function K can undergo discontinuities when u or u goes through an insertion point at u 1 or u 2 . To take into account these discontinuities, it is useful to decompose the kernel along the lines of (5.58), We get in this way a closed system of equations involving C 1 and C 2 only, (5.101) The ≈ signs in (5.99) and (5.100) emphasize the fact that this is a truncated version of the exact equations, to be used only to evaluate the late time behaviour of the OTOCs. These equations can be presented in a convenient matrix form where and we define the -product between components in (5.102) as The above formulas can be applied in any melonic model with ladder kernel K. For the special cases of the U(n) 2 and U(n) models (3.129) and (3.130), the right-hand sides of (5.101) can be conveniently expressed in terms of the retarded and advanced functions defined in (5.14) and (5.15), as well as the particular Wightman two-point correlators (5.108)

Numerical implementation
The equations (5.99) and (5.100), or equivalently (5.102), can be used to compute C 1 (t 1 , t 2 ) and C 2 (t 1 , t 2 ) numerically by introducing a cutoff t max , discretizing the real time variables (t 1 , t 2 ) ∈ [0, t max ] 2 , and inverting the resulting finite size matrix. Note that the integrals effectively extend from zero to t 1 or t 2 due to the presence of the Θ functions in (5.101), so that this procedure allows us to obtain C 1 (t 1 , t 2 ) and C 2 (t 1 , t 2 ) for all t 1 , t 2 < t max without further truncations. The choice of the cutoff is also related to the support of the kernels (5.107) or (5.108), which can be estimated at low temperatures using the conformal forms derived from (5.8) with ρ(ω) given by (3.53), or alternatively determined beforehand from the explicit numerical solution for the real time two-point functions.
In order to see the exponential growth of the OTOCs we need to reliably evaluate C 1 and C 2 up to long times. To accurately approximate the integrals and achieve an appropriate resolution we may therefore need a significant number of points P in the [0, t max ] grid. The runtime of the algorithm is dominated by the P 2 × P 2 matrix inversion, which for the generic matrices we are dealing with here can be done in O(P 2δ ) for some 2 < δ < 3. This is the main reason why the truncated system (5.99)-(5.100) is much more convenient than the exact equation (5.77), since the latter involves more correlators and integrals thus enlarging the final size of the matrix to be inverted.
In practice, our Mathematica implementation allows for P ∼ 10 2 points in the interval [0, t max ], resulting in the inversion of ∼ 10 4 × 10 4 matrices performed in a few seconds with a desktop computer. 31 The typical time scale beyond which the exponential behaviour of C 1 and C 2 is seen is of the order t qn , and we have shown in Figs. 42a and 43a that t −1 qn = − Im ν qn behaves linearly with T at lowtemperatures in the SYK-like phase of the U(n) 2 and U(n) models. We are thus required to take t max ∼ β, so that the resolution requirements then impose a lower bound on the temperatures we can reach. We find that this bound is of the same order of that imposed by the real time algorithm used for the computation of the two-point functions themselves, so it is only occasionally a limitation.
Finally, we should note that while the physical correlators defined by (5.94) are better motivated, there is a practical advantage associated to the standard choice corresponding to (5.95). Indeed, in the latter case the system (5.102) decouples when taking the combinations C 1 ± C 2 , further reducing the size of the problem. Moreover, C 1 + C 2 is real and displays exponential growth, whereas we find that |C 1 − C 2 | is bounded and need not be computed if only Lyapunov exponents are sought. We show in Fig. 46 the combination C 1 + C 2 , obtained from (5.102) as explained above, for typical values of the parameters m and T . 32 The exponential behaviour is clearly seen at long times and yields a precise determination of the associated Lyapunov exponents through linear fitting in log-scale. Note that the values at short times are not reliable, since we use the truncated system (5.102) and not the exact equations (5.77).
We can also use an alternative method to obtain the Lyapunov exponents, due to Kitaev, that can be justified in the following way. We start from the truncated system (5.99)-(5.100) in which we plug the following ansatz, for some a priori unknown functions D 1 and D 2 . In the limit of very large t 0 , the first term on the right-hand sides of (5.99) and (5.100) drops. It is convenient to use 31 Much improvement can be expected from a GPU implementation, which might enable the numerical solution of the full set of equations (5.77). We found however that this was not necessary for our purposes here. 32 The truncated system (5.102) can be used only to compute the exponentially growing contributions to the correlators and thus does not yield a reliable evaluation of the combination C 1 − C 2 . All we can say in this case is that |C 1 − C 2 | is not exponentially growing. the new variables More explicitly, if used. In practice, we find P ∼ 10 3 and ε = 10 −5 to be feasible on a desktop computer within a few seconds.
We used both methods in all our evaluations of the Lyapunov exponents below, and the results turn out to be perfectly consistent.
Remark : Kitaev's method is sometimes justified by taking the large time t = t 1 = t 2 limit directly in (5.99) and (5.100), or (5.102), and using the quasi-normal behaviour of the two point functions to neglect the first term on the right-hand sides of these equations. This is not correct. In particular, the matrix K 1,1 K 1,2 K 2,1 K 2,2 does not have an eigenvalue one. If it had one, the system (5.102) would not be invertible! Note also that the integrals in (5.99), (5.100) or (5.105) run from zero to infinity, unlike the integrals in (5.111) or (5.113) for which one integrates from minus infinity to plus infinity.

Applications
Zero-temperature limit and chaos bound In Fig. 47 and 48, we plot η L /(2πT ) as a function of inverse temperature β, for various values of the mass parameter in the U(n) 2 and U(n) models. These results were computed for the standard case (5.95), but complete agreement is seen with the Lyapunov exponents of the physical one (5.94), as exemplified for m = 0.2 in Fig. 49. β ) computed for SYK in [9]. The bound is saturated in the SYK-like phase in the zero-temperature limit. β ) computed for SYK in [9]. Saturation of the bound is only clearly seen for m = 0 and seems to be violated for all m > 0, including in the gapless phase 0 < m < m * .
In the case of the U(n) 2 model, the results are as expected. We find that the bound η L ≤ 2πT of (5.93) is saturated for the SYK-like phase in the zero-temperature limit, including for m > m * where this solution is metastable (e.g. for m = 0.24 in Fig. 47).
The situation seems to be very different in the case of the U(n) model. As soon as m > 0, we find that the curve for η L /(2πT ) reaches a maximum and then starts to decrease as the temperature is lowered. Our numerics do not allow us to probe the deep IR regime, so we must remain cautious about the conclusions we may draw for the strict zero-temperature limit. However, the behaviour is clearly qualitatively different from the U(n) 2 model and suggests that the bound (5.93) is not saturated already in the phase where the non-standard ansatz discussed in Sec. 3.1.5 is valid. The fact that the SL(2, R) subgroup of the reparameterization symmetry (3.38)-(3.39) is broken, including the scale invariance, is certainly crucial to understand this property. This seemingly very different behaviour of the Lyapunov exponent is contrasting with other properties that we have discussed previously, like the existence of a nonvanishing zero-temperature entropy, that are similar to the SYK-like phase of the U(n) 2 model. Clearly, this raises interesting questions with respect to the gravitational interpretation of the U(n) model, in particular in relation with a possible description in terms of a suitable generalization of the Schwarzian action, and deserves further study.  .95), with values for the physical case (5.94) represented with crosses and circles, for the U(n) 2 and U(n) model respectively. All results were computed using Kitaev's method with P = 1001 and ε = 10 −5 . As discussed in Sec. 5.2.4 equations (5.102) can be decoupled in the standard case, leading to a faster and more precise solution allowing for slightly lower temperatures. An equally good agreement is observed between Lyapunov exponents for standard and physical correlators for other masses, both in SYK-like and HO-like phases.
Remark : The non-saturation of the chaos bound discussed above in the gapless phase of the U(n) model must not be confused with the usual and expected non-saturation of the bound found in the gapped phase, for m > m * . The gapped phases of both the U(n) 2 and U(n) models are similar to what is found in the SYK model with a random mass term or in the Maldacena-Qi model, where non-saturation of the chaos bound was already observed [39,55,66].
Critical properties of the Lyapunov exponent We close our analysis with a discussion of the very interesting behaviour of the Lyapunov exponent near the critical point (m c , T c ) of the U(n) 2 model. In Fig. 50a, η L /(2πT ) is plotted as a function of m for T = T c . We find a singular behaviour near m = m c , consistent with We have also studied η L /(2πT ) as a function of T at m = m c , see Fig. 50b. It is very plausible that critical exponents can be defined in this case as well, but our data is not precise enough to extract them in a reliable way.  Fig. 29. In (b) a singular behaviour is clearly seen when T → T c , but numerical data is not precise enough to extract critical exponents in a reliable way.

A q-body interacting Hamiltonians
One can consider the following extension of the U(n) 2 model, which has interesting features discussed in Sec. 4.2. We generalize the random Hamiltonian (2.17) to where the sum is over even integers q ∈ q = {q 1 , . . . , q p }, 4 ≤ q 1 < q 2 < · · · < q p , and the λ i 1 ···i q/2 j 1 ···j q/2 are random Gaussian couplings satisfying constraints similar to the λ kl ij in (2.14) and (2.15). One can also provide a matrix model version of the model, using the complete interaction vertex studied in [4]. The effective action (3.86) is generalized to  If only one coupling λ q is turned on, these formulas simplify to F N = m − T ln 1 + e βm + T

B Algebraic derivation of the kernel formula
The goal of this Appendix is to provide all the elementary details of the algebraic derivation of the formulas (3.127) and (3.133) given in Sec. 3.4 for the connected four-point function.
We start from the equations (3.139), (3.140) and (3.142). In order to simplify the notation we write functions or operators like G(t 1 , t 2 ) or O Σ as G ij or (O Σ ) ij , integrals over t are replaced by sums and functional derivatives by partial derivatives. In these notations, and using (3.142) to express the free energy in terms of the on-shell effective action, the four-point function (3.139) is given by For the models we are studying, the relevant s[G] functionals are (3.143) and (3.144), or in the notation of this Appendix

C On the numerical errors
Many of the results presented in this paper are numerical in nature, consisting either of specific values for various relevant quantities (critical temperatures, masses, exponents, etc.) or observations concerning the solutions of Schwinger-Dyson equations, such as e.g. the existence of monodromies in the space of parameters of the U(n) 2 model. It is therefore important to understand the limitations of the numerical methods employed, which we briefly discuss in this Appendix.
Lacking statistical fluctuations in our results given the exclusive use of deterministic algorithms, there are two main classes of numerical errors or uncertainties to consider: i) Systematic errors are incurred in whenever imposing a cutoff to truncate an infinite set of equations or realize a limiting procedure, such as in (3.113) for the Euclidean setting and (5.40) for the real time one; or introducing a discretization to approximate continuous quantities, as in (5.42) for the real time case. This is essentially unavoidable, and completely trumps the numerical errors inherent to the machine-precision floating-point arithmetic used during all computations. 33 ii) Methodological uncertainties are present whenever fitting numerical results to analytical models, either for parameter extraction (scaling dimensions, critical exponents, etc.) or for zero-temperature extrapolations. They also naturally appear whenever using binary search to determine quantities such as critical temperatures and masses, or Lyapunov exponents as explained in Sec. 5.2.5.
As already mentioned in Sec. 3.2.3 and 3.3, systematic errors translate into any and all quantities computed using the numerical solutions we obtain. While we have not rigorously bounded them in terms of the various parameters involved, they can be roughly estimated by comparing with exact analytical results whenever possible. For example, to estimate the precision of the thermodynamic functions computed for the U(n) 2 model at T = 10 −3 in natural units λ = 1, we may choose a large enough mass to make sure that the only solution to the Schwinger-Dyson equations is the perturbative one, very close to (3.7) at this temperature. In this case taking m = 1 is certainly enough, so that using a cutoff of k UV = 2 16 we find F/N 6 × 10 −5 instead of F/N = 0 as expected for (3.7) at T = 0. We may then assume that a systematic error of ∼ 10 −4 is present when computing free energies at T = 10 −3 with a cutoff k UV = 2 16 , also for other large masses in the U(n) 2 model. An analogous procedure can be carried out for the entropy at small masses comparing with the exact formula (3.111), yielding an estimate of ∼ 0.1 for the systematic error in the numerical results when using a cutoff k UV = 2 16 at T = 10 −3 . This is consistent with the β = 10 3 factor relating the entropy and free energy in (3.104), so that the systematic errors are observed to be approximately mass-independent in this case, at least for large enough cutoffs and low enough temperatures.
The methodological uncertainties, on the other hand, originate in most cases both in the fitting models and techniques we use, as well as in the actual data we choose to analyze. The effect of the latter is generally the most relevant one, with the statistical error of the fitting procedures being considerably smaller than the typical fluctuations observed when modifying the underlying data. While we have not attempted here a rigorous treatment of these issues, we can again roughly estimate their effects. For example, when computing scaling dimensions at small temperatures by fitting the zero-temperature ansätze (3.42) or (3.55), as shown in Fig. 10 and 14, we need to choose a fitting interval satisfying 1 τ β in natural units λ = |ξ| = 1. One simple choice is to take a 10% section of the full Euclidean time interval, say τ /β ∈ [0.2, 0.12] for some large β, and keep track of the fluctuations of the fitted scaling dimensions when considering 10% fluctuations of its endpoints, i.e. when fitting the ranges τ /β ∈ [0.1, 0.11] and τ /β ∈ [0.3, 0.13]. Similar considerations apply to the calculation of critical exponents, where there is an additional uncertainty in the critical temperatures and masses to be used, so that fluctuations in these parameters need to be accounted for in an analogous manner.
Other error estimation methods are of course possible, but the important observation here is that results show certain universality: for all reasonable choices, even if arbitrary, the numerical fluctuations we encounter are consistent. The methodological error bar obtained in this way is then fairly robust. Moreover, when the method above is applied systematically it serves the additional purpose of providing an indicator for relative certainty in these determinations. It is in this sense that we may begin to quantify the naked-eye impression that the determination of ∆ − in Fig. 14 is less precise than that of ∆ + .
In practice, both types of numerical errors and uncertainties described above are intertwined. For example, the systematic error in numerically computed thermodynamic functions prevents us from distinguishing branches of solutions to the Schwinger-Dyson equations with arbitrary precision, thereby forcing us to prematurely interrupt binary searches for critical temperatures (thus prompting a methodological uncertainty). Alternatively, we may use extrapolation procedures as shown in Sec. 4.2.2 to account for sytematic errors, obtaining numerical results that have better precision than one could naively expect, see e.g. Fig. 25.
The error bars we report in the various quantities discussed in this work account for the methodological uncertainties, estimated as explained above, and for systematic errors only as far as they are encompassed by the former. They are therefore not exhaustive, and subject to vary at least slightly for different data and methodologies. It is for these reasons we do not present them in figures, but provide them only for the most important values explicitly referenced in the main text. In any case, the general features and conclusions we derive from these values never depend on the fine-grained details. Indeed, wherever this could be an issue, such as in the detailed characterization of the low-temperature conformal phase of the U(n) model, we opt for the more conservative route and only provide definite assertions if the numerical results are conclusive beyond reasonable doubt.