{'fraction_non_alphanumeric': 0.039109754603277154, 'fraction_numerical': 0.011050093012375658, 'mean_word_length': 4.227624759797356, 'pattern_counts': {'":': 1, '<': 28, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The effect of quantum fluctuations on the coloring of random graphs

We present a study of the coloring problem (antiferromagnetic Potts model) of random regular graphs, submitted to quantum fluctuations induced by a transverse field, using the quantum cavity method and quantum Monte-Carlo simulations. We determine the order of the quantum phase transition encountered at low temperature as a function of the transverse field and discuss the structure of the quantum spin glass phase. In particular, we conclude that the quantum adiabatic algorithm would fail to solve efficiently typical instances of these problems because of avoided level crossings within the quantum spin glass phase, caused by a competition between energetic and entropic effects.

We present a study of the coloring problem (antiferromagnetic Potts model) of random regular graphs, submitted to quantum fluctuations induced by a transverse field, using the quantum cavity method and quantum Monte-Carlo simulations. We determine the order of the quantum phase transition encountered at low temperature as a function of the transverse field and discuss the structure of the quantum spin glass phase. In particular, we conclude that the quantum adiabatic algorithm would fail to solve efficiently typical instances of these problems because of avoided level crossings within the quantum spin glass phase, caused by a competition between energetic and entropic effects.

I. INTRODUCTION
Among the fascinating potential applications of quantum mechanics, the promise of quantum computing is to make use of the laws of quantum mechanics to enhance computers' calculation power. Besides the great effort of research towards the physical realization of quantum computers, a lot of activity has been devoted to develop quantum algorithms and understand their efficiency. A central problem encountered in almost all branches of science is to optimize irregularly shaped cost functions: the Quantum Adiabatic Algorithm (QAA) [1][2][3][4][5] is a generic and universal procedure to tackle such problems. Suppose one wishes to find the ground state of a Hamiltonian H P acting on N qubits. The idea of the QAA is to implement an interpolating Hamiltonian H(Γ) = H P + Γ H Q such that the quantum computer can easily be initialized in H Q 's ground state, and to decrease the interpolation parameter Γ from a very large value down to zero. If this is done slowly enough the adiabatic theorem [6] ensures that with high probability the system remains at all times in the ground state of the interpolating Hamiltonian. In particular, at the end of the evolution, it is in the ground state of H P and the original problem is solved. The crucial question is how slow the evolution should be in the thermodynamic limit N → ∞, and whether the total evolution time T has to grow polynomially or exponentially fast with N . This is reminiscent of the classical complexity classes [7,8], which indeed have quantum counterparts [9][10][11]. Quite generally, the adiabaticity condition requires T to be larger than the inverse of the squared gap between the ground state and the first excited state of H(Γ), for all Γ. Hence, the efficiency of the QAA mainly depends on the rate of closing of the minimal gap along the interpolation, even though more subtle issues such as determining the residual energy after the interpolation [12,13] would require a further understanding of the annealing dynamics.
The possibility that the QAA could outperform classical algorithms has triggered a lot of work on optimization problems in quantum fields (see [14][15][16][17] for recent reviews), trying in particular to pinpoint the rate of closing of the minimal gap along the interpolation path. As it has now become usual, these studies used random ensembles of constraint satisfaction problems [18,19], focusing on results valid with high probability in the thermodynamic limit. Early results generated excitement by reporting polynomial scaling of the minimal gap for some classically hard problems [5,20,21]; however, these results where hampered by strong finite-size effects or by the fact that the instances considered were not typical of the underlying problem. On the other hand, [22][23][24][25] obtained negative results for the success of the QAA on a certain class of optimization problems: namely, they found a first order quantum phase transition during the interpolation between the quantum and the classical Hamiltonians. Such a phase transition is not surprising in the context of fully-connected quantum spin glasses [26][27][28][29], and is known to generically lead to a gap vanishing exponentially fast with the system size N [22][23][24] and thus to a blow-up of the time needed by the QAA. Failure of the QAA because of another mechanism was also discovered for different models in [30][31][32]: because the classical energies have a non-trivial perturbative expansion in Γ, the energy of an excited state may decrease much faster than the one of the ground state, leading to avoided "perturbative" crossings near the classical end of the algorithm. However, these results were obtained either on a toy model, or using well-chosen Hamiltonians H P with a peculiar structure of the low-energy spectrum of H P , that were not typical of the optimization problem considered.
Henceforth, one would like to understand what happens more generically for models with multiple and not necessarily isolated ground states. This is expected to happen for typical optimization problems, which are known to possess complex and intricate configuration spaces that were unveiled by the use of statistical mechanics tools [33][34][35][36]. In these more common cases, do the avoided crossings remain finite and isolated in Γ (hence leading to singularities in the ground state energy for N → ∞) or do they proliferate and accumulate, leading to a continuum of level crossings and a gapless phase? Hints for an answer have been obtained on a toy model in [17,32]; this paper aims to highlight their conclusions on a more realistic optimization problem. This question is also important from a practical point of view as it has been argued that a finite number of level crossings could be eliminated by a suitable redefinition of the quantum Hamiltonian [37][38][39][40]. Also note that the Xorsat problem studied in [24] already possessed in some cases an exponential degeneracy of its ground states, but they had a particular structure that smoothened the effect of quantum fluctuations.
The model that we shall study is the coloring one, a famous problem in combinatorics, which is known to be classically very hard to solve -more precisely, NP-complete [41]. Given a graph and q col colors, it consists in coloring the vertices such that no edge connects two vertices of the same color. In physical terms, this corresponds to a q col -states antiferromagnetic Potts model. This model can be studied on any kind of graph, but we shall focus in the following on random regular graphs, with four colors (q col = 4) and connectivity c = 9 and c = 13. Thanks to the quantum cavity method [42,43], we will compute the phase diagram of this model in both cases and unveil the nature of its quantum spin glass phase. Our main results are (i) the nature of the quantum phase transition that occurs along the interpolation from the quantum phase to the classical one at low temperature, that we find to be continuous in the spin glass language and thus of third order thermodynamically, and (ii) the presence of a continuum of level crossings within the spin glass phase, that are induced by entropic effects and in particular by the clustered structure of the spin glass phase. This last feature should lead the QAA to fail to solve the coloring problem efficiently. Our results will also be corroborated by Monte-Carlo simulations. Note that result (ii) confirms the predictions made in [32] on a much simpler toy model.
The plan of the paper is as follows. We first briefly recall in Section II the classical features of the model and the definition of its quantum version. Section III presents the phase diagrams that we obtain for the quantum coloring problems, while Section IV presents in greater details the structure of its quantum spin glass phase. Technical details are deferred to a series of Appendices. A brief part of our results (mostly Sec. IV A) already appeared on the review paper [17].

II. THE COLORING PROBLEM: CLASSICAL PICTURE AND DEFINITION OF ITS QUANTUM VERSION
A. The classical coloring problem Let us consider a graph G = (V, L) with V a set of N vertices and L a set of M edges between pairs of vertices. We introduce a Potts variable σ i ∈ {1, . . . , q col } on each vertex i of the graph, and denote σ = (σ 1 , . . . , σ N ) the global configuration of the variables. To each of these configurations we associate an energy, or cost function, where the sum runs over the edges of the graph, and δ denotes the Kronecker symbol. Interpreting the Potts variable σ i as a color given to vertex i, this cost function counts the number of monochromatic edges in the configuration σ.
In physical terms there are antiferromagnetic interactions between pairs of variables linked by an edge of the graph. The Gibbs-Boltzmann probability measure at inverse temperature β is then defined as with the partition function Z(β) = σ e −βE(σ) ensuring its normalization. As mentioned in the introduction the coloring problem is NP-complete (for q col ≥ 3) in its decision version, namely there does not exist any algorithm able to decide the existence of a proper coloring (a configuration of zero energy) in a number of operations bounded by a polynomial in N, M for every graph G. Computer scientists [19] thus turned to random ensembles of graphs in an attempt to study the typical difficulty of the coloring problem, typical meaning here "with a probability going to one in the large size (thermodynamic) limit N → ∞". The interesting regime is the one of sparse random graphs, with the number of edges M of the same order than N . The absence of an underlying Euclidean structure in these random graphs makes the corresponding statistical model of the mean-field type, hence the methods devised for the study of mean-field spin glasses models [44] could be applied to characterize the typical behavior of the partition function Z and of the measure µ. For the coloring problem, see [45] for the use of the replica method and [35,[46][47][48][49] for the application of the cavity method. The main outcome of these studies is the unveiling of phase transitions in the behavior of the free energy density f = − 1 N β ln Z (that concentrates around its average in the thermodynamic limit) and in the properties of µ, as a function of the inverse temperature β and of the parameters of the random graph ensemble. In the following we shall sketch these phase transitions and the cavity methodology used to derive them; for more details and justifications of the equations the reader is refered to the above quoted references for the coloring problem, and to [36,50] for more generic presentations of the cavity method. We will focus for technical simplicity on the case of random c-regular graphs: the graphs G are chosen uniformly at random among all the graphs on N vertices in which all vertices have the same degree (number of neighbors, or connectivity) c. In the thermodynamic limit (N → ∞), such graphs are locally tree-like, meaning that if one selects one of its vertex at random, the shortest loop around it is larger than any fixed length with a probability which goes to one when N → ∞ [51]. The cavity method exploits this tree-like character, building on the exact solution for finite trees that is easily obtained by recursion, and dealing with the boundary conditions induced by the long loops of the random graphs in a self-consistent way. Consider an arbitrary variable i in a large random c-regular graph, and its marginal probability η(σ i ) that would be obtained from the Gibbs-Boltzmann measure (2) if one of the c edges around it were removed. Forgetting the possible effects of the long loops of the graphs, the translational invariance and the recursive structure of the tree implies that η is solution of the following self-consistent equation, called replica-symmetric (RS) cavity equation: Here z({η i }) is a normalization factor ensuring that σ g(η)(σ) = 1. The cavity method then proposes an expression for the free energy density f in terms of the solution η of this self-consistent equation. The equation (3) and the associated RS prediction for the free energy density are correct when the spatial decay of correlations between variables is fast enough, which is equivalent to the existence of a single pure state in the Gibbs-Boltzmann measure (or of a finite number of them, related by simple symmetries). At low enough temperature and in presence of many constraints between variables (i.e. when c is large enough) this assumption breaks down. The Gibbs-Boltzmann measure is then split in an exponentially large number of pure states (or clusters), the correlation decay assumption underlying the RS computation is then valid only inside one pure state, not in the complete Gibbs-Boltzmann measure. The 1-step of replica-symmetry breaking (1RSB) cavity method deals with this situation by making a further ansatz on the structure of the pure states. One assumes that the typical number of pure states of internal free energy density f is exponentially large, at the leading order e N Σ(f ) , with a rate function Σ(f ) called complexity or configurational entropy. This function is computed via its Legendre transform [52] φ(m) = min f {f − T Σ(f )/m}, the parameter m conjugated to f being called the Parisi replica-symmetry breaking parameter. The thermodynamic potential φ(m) is obtained from the solution of the 1RSB cavity equation that generalizes (3), the order parameter solution of this self-consistent equation becoming P (η), a distribution of the RS probabilities η with respect to the choice of the pure states weighted proportionally to e −N mβf : where Z is a normalization constant and the functions g and z are defined in Eq. (3). One can check that the replica-symmetric equation (3) is recovered if P is a delta-peaked distribution. The structure of the Gibbs measure and its decomposition into pure states can be probed by the inter-and intrastate overlaps. The latter are defined, respectively, by: where P is defined by the right hand side of (4), but with a product over c terms instead of c − 1. The energy function (1) is symmetric under the exchange of two colors; hence for all σ, d P (η)η(σ) = 1 q col , and q 0 is always equal to 1/q col . Moreover, on the RS solution, this symmetry enforces η(σ) = 1/q col for all σ. The appearance of a non-trivial 1RSB solution to the equation (4) can then be detected by the fact that the intra-state overlap q 1 takes a value strictly larger than 1/q col .
Depending on the values of the temperature and of the degree c the pure state decomposition of the Gibbs-Boltzmann measure exhibits qualitatively different properties, that can be read on the type of solutions of the RS (3) and 1RSB (4) cavity equations: • if the 1RSB equation with m = 1 admits only a RS solution, then almost all configurations are part of a single pure state, the RS hypothesis is correct, and its prediction for the free energy density is valid. We shall call this situation a RS phase, or paramagnetic phase (P) here.
FIG. 1: Classical phase diagram of the coloring problem for q col = 4 in the connectivity-temperature plane, see [48,49] for the actual numerical values. The solid orange line separates the paramagnetic phase (P) from the dynamical paramagnet (dP) while the dashed blue lines marks the boundary of the genuine spin glass (SG) phase. The model is defined only for integer connectivites (indicated by squares and circles); lines serve as a guide to the eye. The connectivites c d = 9 and cK, cs = 10 are defined in the text. Pink arrows indicates the connectivities that will be studied in the quantum case.
• if the 1RSB equation with m = 1 admits a non-trivial solution, and if the min in the definition of φ(m = 1) is reached at f int with Σ(f int ) > 0, then the system is in a dynamic 1RSB (d1RSB) phase, or dynamic paramagnetic phase (dP) here. Almost all configurations belong to pure states of internal free energy density f int , and there are e N Σ(fint) such pure states. It turns out that the total free energy density f int − T Σ(f int ) coincides with the RS prediction, yet the splitting of the configurations into pure states has drastic consequences on the dynamics of such a model, that becomes non-ergodic (on sub-exponential timescales), hence the name of dynamic 1RSB phase.
• if the 1RSB equation with m = 1 admits a non-trivial solution, but predicts a negative complexity Σ(f int (m = 1)) < 0, then the Parisi breaking parameter m has to be set to the value m s ∈ [0, 1] such that Σ(f int (m s )) = 0. In such a 1RSB phase almost all configurations belong to a sub-exponential number of pure states of internal free energy density f int (m s ); as the complexity vanishes this value is also the 1RSB prediction for the total free energy density of the model. Thise phase corresponds to a genuine spin glass (SG) phase.
In general these three phases are encountered in this order upon lowering the temperature (for large enough degrees) or increasing the degree (for low enough temperature). The transition between the RS and d1RSB phase occurs at the dynamic temperature T d (c) (or at c d (T )), and is not accompanied by any singularity in the free energy. The transition between the d1RSB phase and the 1RSB phase, at which the Gibbs-Boltzmann measure condenses on the few largest existing clusters, occurs at T K (c). This transition corresponds indeed to the Kauzmann transition in the Random First Order Theory (RFOT) of structural glasses [53,54], which is a second order transition from the thermodynamic point of view, even if the order parameter jumps discontinuously ("first order") when the transition is crossed. In this context the dynamic transition T d corresponds to the divergence of the relaxation time of the Mode Coupling Theory (MCT) of supercooled liquids [55]. This pattern of transitions has indeed been found for the coloring of random regular graphs with q col ≥ 4 colors [48,49], and is illustrated on Fig. 1. The case q col = 3 is singular, there is then no intermediate d1RSB phase, the transition between the RS and the 1RSB phase occurs continuously, via a local instability of the RS solution towards the 1RSB one at a temperature T i (c). In the context of optimization problems the zero temperature limit plays a particularly important role; the transition lines defined above end in this limit at the critical connectivities c d and c K . A further transition can be defined at zero temperature: the satisfiability threshold c s is such that the ground state energy vanishes for c ≤ c s (in terms of the coloring problem the typical answer to the decision question is yes, there are proper colorings of such graphs) and is strictly positive for c > c s .
In the following we shall study the quantum version of this model for q col = 4 and c = 9 or c = 13; these two cases correspond to the situations illustrated by arrows on Fig. 1. Before turning to this quantum framework let us make a few remarks. The coloring problem can be studied for other ensembles of random graphs, notably for Erdös-Rényi random graphs in which the degree of vertices are Poisson random variables of average c. The transitions described above are also encountered in this case [48,49], with the difference that c is now a continuous parameter. We restricted ourselves to the regular case for technical reasons: to deal with the fluctuating degrees in the Erdös-Rényi case one has to introduce a more complicated order parameter, i.e. in the 1RSB cavity treatment a distribution over the distributions P to deal with the fluctuations of the probabilities η with respect both to the choice of the pure states and of the local connectivity. Although this complication is affordable in the classical case, it would become impossible to treat in the quantum setting. We also neglected in this brief presentation the phenomenon of full replica-symmetry breaking [44]: the 1RSB description of the Gibbs-Boltzmann measure is only the first one in a hierarchical construction in which the pure states are themselves grouped in clusters of pure states, and so on and so forth, that become relevant at a so-called Gardner transition [56], computed in [48,49] for the classical antiferromagnetic Potts model. This whole hierarchy can be dealt with in fully-connected models, most notably the Sherrington-Kirkpatrick one, yet it cannot be handled with present techniques in finite connectivity models (defined on sparse random graphs), even for classical variables, not to speak about quantum models. The full RSB structure allows to study some important physical properties of these models, most notably the marginal stability of glass states; however, the 1RSB approximation is expected to give quantitatively good predictions for thermodynamic properties even when it is unstable towards a full RSB one, and therefore it will be sufficient for the purposes of this study.

B. Definition of the quantum version of the coloring problem
In order to define the quantum version of the problem we introduce the q N col -dimensional Hilbert space spanned by {|σ , σ ∈ {1, . . . , q col } N }. We then define the "problem" Hamiltonian H P corresponding to the cost function (1) as the operator diagonal in this (so-called computational) basis with diagonal elements equal to the classical energies: We introduce quantum fluctuations, i.e. off-diagonal matrix elements, by defining for each site i an operator T i that flips the color σ i to any other different color: This is the simplest generalization from the Ising spin case (that corresponds to q col = 2) of the Pauli matrix σ x i (taking as the computational basis the eigenvectors of the σ z i matrices). The intensity of the quantum fluctuations is controlled by the "transverse field" Γ, the total Hamiltonian reading The partition function at inverse temperature β is Z(β, Γ) = Tr[e −β H(Γ) ], thermodynamic averages being denoted • = Tr[• e −β H(Γ) ]/Z(β, Γ). In particular we will call m x = T i ∈ [0, q col − 1] the "transverse magnetization", by analogy with the Ising spin case. This also explains the name of "transverse field" for Γ.

C. A sketch of the quantum cavity method
A standard way to compute the partition function of such a quantum model consists in using the Lie-Suzuki-Trotter formula to disentangle the two non-commuting parts of the Hamiltonian H(Γ), introducing copies of the original degrees of freedom σ i along an "imaginary time" axis of length β. In the limit where the number of such copies goes to infinity one thus obtains an exact representation of the quantum partition function as a path integral of a classical model, the price to be paid being the replacement of the spins σ i ∈ {1, . . . , q col } by more complicated degrees of freedom σ i (we emphasize the use of a bold font here), which are piecewise constant periodic functions σ i (t) from [0, β] to {1, . . . , q col }.
Apart from this replacement the classical model thus obtained from the quantum one has the same topology of interactions as the classical one (1); in particular if the latter can be treated with the classical cavity method (i.e. if it is defined on a locally tree-like random graph), then its quantum version can be handled with the quantum extension of the cavity method where the spins σ i are replaced by the imaginary time trajectories σ i . This observation was first put to work in [43] with a finite number of Suzuki-Trotter slices, the continuous imaginary time limit being taken in [42]. These two works developed the quantum cavity method at the RS level, and the quantum 1RSB framework was then introduced in [24,57]. Detailed explanations and derivations, covering the case of the coloring in presence of a transverse field, can be found in the review [17], we shall thus content ourselves here with a few remarks, some technical details being presented in App. B, in particular the parallelization of the numerical code we used to accelerate the resolution of the cavity equations.
As explained above the quantum cavity method has exactly the same structure as the classical one, sketched in Sec. II A. The main modification is the fact that the basic object η appearing in the RS (3) and 1RSB (4) cavity equations is now a probability distributions over the space of piecewise constant periodic functions from [0, β] to {1, . . . , q col }. This space is obviously much larger than {1, . . . , q col }, and in consequence a single η, that could be represented by q col − 1 real numbers in the classical case, has now to be represented in an approximate way for the numerical resolution of Eqs. (3), (4). A convenient representation is provided by a finite sample of N traj random trajectories, η being approximated by the empirical distribution of this sample.
Finally, the quantum overlaps are defined in analogy with the classical case, taking an extra average over the imaginary time: As in the classical case, because of the symmetry between colors, q 0 = 1/q col , and q 1 ≥ 1/q col , with equality on the RS solution.
D. Expected features for the quantum model

Quantum phase transitions
The ground state of the classical coloring problem is diagonal in the computational basis, while on the other side of the interpolation the ground state of the quantum Hamiltonian H Q is diagonal in the tensorial product of the eigenvectors of the flipping operators T i . This defines two phases with very different physical properties. Slowly interpolating from high Γ down to Γ = 0, the spins have to rotate to go from the quantum paramagnetic phase to the classical spin glass one. If the classical Hamiltonian was simple or resembled the quantum one, for example if it acted as a product of single spin operators, this rotation would be easy and nothing special would happen along the way; but in our case, because of the complicated nature of the interaction (6) and of its associated classical phase, one can expect collective effects to become strongly relevant and a quantum phase transition [58] to occur along the way.
It is well established that the gap vanishes at least polynomially fast with N at a quantum second order critical point [25,58] (it vanishes exponentially fast in some cases in presence of disorder [59]), while it vanishes exponentially in N at a first order phase transition [22][23][24]. As the scaling of the minimal gap is, according to the adiabatic theorem [6], critical for the success or failure of the QAA, we have to determine in our case the order of this quantum phase transition. This will be the subject of Sec. III.

Level crossings and the role of entropy
Apart from a quantum phase transition at relatively large values of the transverse field, we also expect the addition of the transverse field to have very important effects also for much smaller Γ. A first example of such a phenomenon was put forward in [30,31]; we sketch it very briefly here. Consider an instance of an optimization problem which has an isolated solution, and another local minimum of its cost function with only one violated clause, and such that these two configurations are far apart in Hamming distance. Computing the continuations of these energies using standard perturbation theories, one can find that, if the model is well chosen, they cross for some value of Γ. This value shrinks to zero as N → ∞, hence the name of perturbative crossings in the literature. This mechanism is for the moment only understood in perturbation theory, and therefore it holds whenever perturbation theory holds, that is, if at small enough Γ the full eigenstates of the quantum problem remain close enough to the classical eigenstates. The other important ingredient is a very particular construction of the instances of the problem, that admit only one solution. But typical instances of generic random optimization problems, even close to the satisfiability threshold, have an exponentially large number of solutions, and so we expect that non-degenerate perturbation theory should not hold [60] and that the spectrum should be much more complex. Therefore we would like to understand what happens generically in problems that have multiple and not necessarily isolated solutions.
The quantum random subcubes model [17,32] is built as a quantum extension of the classical random subcubes model [61] (note that this model makes use of Ising spins, but these could be changed into Potts spins with any number of states q col , with only irrelevant changes in numerical prefactors). In this model, the classical spin glass phase is represented by a set of random subcubes of the Hilbert space (representing clusters). These subcubes are supposed to be disjoint, of different sizes, and to be associated with random (classical) energies. Adding quantum fluctuations to this model gives rise, under reasonable hypotheses on the distribution of the sizes of the subcubes, to a series of level crossings induced by a combined energetic-entropic effect. More precisely, if one introduces the internal intensive entropy s of a cluster (such that a cluster contains e N s configurations), the complexity Σ(e, s) such that there exist e N Σ(e,s) clusters of intensive energy e with such an entropy, and s max (e) = sup{s, Σ(e, s) ≥ 0}, then it can be shown that the ground state energy of the model is given by: Henceforth, as soon as Γ > 1/((ln 2)s ′ max (e GS (Γ = 0))), the minimum is obtained for a different value of e for each value of Γ. In this region the ground state changes abruptly from one cluster to another upon changing Γ by an infinitesimal amount, similarly to what is called temperature chaos in classical spin glasses [62,63]. Because the clusters are at Hamming distances proportional to N , we expect these crossings to be avoided at finite N by producing exponentially small gaps.
The analysis of the quantum random subcubes model can also be done at finite temperature, which reveals the existence of a condensation transition at a temperature T K (Γ) (similar to the one discussed in Sec. II A). The effect of the quantum fluctuations is particularly drastic when the classical model, at Γ = 0, remains un-condensed down to zero temperature (i.e. when T K (Γ = 0) = 0). Indeed in this case, for infinitesimal values of Γ > 0 the model suddenly condenses on the largest clusters of classical groundstates, because the latter undergo the largest decrease in energy when the transverse field is turned on. One can study more precisely the limit Γ, T → 0, in which the partition function is equivalent to The slope of the condensation line T K (Γ) in the small Γ limit can then be deduced from a saddle-point evaluation of this expression: it corresponds to the critical value of βΓ such that the maximum of the argument of the exponential is reached at the upper limit of integration s max (e GS (Γ = 0)), i.e.
As discussed in Sec. II A, we expect the coloring problem for q col = 4 and c ≥ 9 to possess a clustered classical phase and thus to exhibit a phenomenology very close to that of the quantum random subcubes model; in particular, crossings within the spin glass phase (see Sec. IV) and a linear condensation temperature T K (Γ) at small field for c = 9 (see Sec. III B). The fact that we shall observe these phenomena in the coloring problem confirms the relevance of this scenario, first proposed in [17,32], for generic random optimization problems. Let us finally briefly explain why the same properties did not appear when studying the random regular Xorsat problem [24], even when it has a degenerate ground state (for example 4-Xorsat on a regular random graph of connectivity 3): because such a formula is locked [64], each solution is isolated and no entropic effects appear when adding quantum fluctuations. Hence the quantum coloring case is the first realistic optimization problem displaying the phenomena outlined above.

III. THE PHASE DIAGRAMS OF THE COLORING MODEL IN A TRANSVERSE FIELD
In this section we present our results concerning the phase diagrams of the quantum model as a function of the temperature T and of the transverse field Γ, for two different connectivities c, emphasizing the modifications of the classical transitions (at Γ = 0) once the quantum fluctuations are turned on. Most technical and numerical details of the analysis are deferred to the Appendices. A further investigation of the structure of the spin glass phase is presented in Sec. IV. Let us start with the case c = 13. We recall that in the classical case (Γ = 0) the model exhibits a clustering (dynamical) transition at T d (Γ = 0) = 0.455, a condensation transition at T K (Γ = 0) = 0.450, and a RS instability transition at T i (Γ = 0) = 0.441 [49] (we will not discuss the Gardner transition towards full RSB at T G = 0.185 [49]). From a thermodynamic point of view the only relevant transition in this classical limit is T K , the free energy has a discontinuity in its second derivatives at this temperature, while it is non-singular at T d and T i . The latter marks the limit of local stability of the RS solution, but it is here irrelevant because of the discontinuous transition towards the 1RSB phase that occurs at higher temperature.
We report on Fig. 2 the phase diagram that we obtain with the quantum 1RSB cavity method. We do find, for Γ large enough, only one transition line, T i (Γ), that corresponds to a local instability of the RS solution towards a non-trivial solution of the 1RSB equation. This transition is continuous, in the sense that the overlap q 1 defined in Eq. (9) grows continuously from its replica-symmetric value 1/q col when one enters the spin glass 1RSB phase. In order to locate it as precisely as possible we resorted to a finite population size scaling analysis, as explained in more details in Appendix C 1. On this line the free energy of the model is singular, and this corresponds to a thermodynamically third order transition. This pattern of a second order RFOT like condensation transition (at T K for Γ = 0) that becomes a thermodynamically third order transition upon changing one parameter was actually met in earlier studies of mean-field disordered models, namely in the fully-connected spherical p-spin model in presence of a magnetic field [65,66]. For the convenience of the reader we gathered in Appendix A the main aspects of the analysis of this (classical) model; the third order character of the thermodynamic transition derives, in both cases, from the linear growth of the overlap with respect to the distance to the continuous transition. A further evidence for the order of the transition is presented on Fig. 3: on the left panel one observes a rather good collapse of the complexity curves Σ(m), for various values of the transverse field Γ at a fixed temperature T , rescaled by (Γ i (T ) − Γ) 3 . These curves vanish at the static Parisi parameter m s (Γ, T ), that admits a finite limit (smaller than 1) when Γ → Γ i (T ). As shown on the right panel of Fig. 3 this value is proportional to T when T is reduced for a fixed value of Γ.
In the classical limit Γ → 0, the line T i (Γ) falls below the classical transitions T d and T K . This means that these two classical transitions should extend into two lines at finite but small Γ, and at some finite Γ the transition must change nature from discontinuous to continuous. However, in this model all this happens in a range of T and Γ too small to be resolved within the accuracy of the numerical resolution of the quantum 1RSB equation.
We should emphasize here that the incorporation of the replica-symmetry breaking effects was crucial to unveil the phase diagram of this model. As a matter of fact the RS computation predicts, incorrectly, a first order phase transition as a function of the transverse field at low enough temperature. We present in Appendix C 3 the detailed numerical results and arguments we used to rule out this spurious prediction of the RS computation. A further confirmation on the continuous nature of the RSB transition will be provided by the Monte-Carlo simulations presented in Sec. IV. We now turn to the case c = 9. Before discussing the quantum case, we recall [49] that for q col = 4 and c = 9, the model is classically satisfiable, meaning that the ground states' energy is zero and that graphs are typically colorable. These ground states are exponentially numerous and are arranged in an exponentially large number of clusters; this corresponds to the region c d ≤ c < c K , meaning that T K (Γ = 0) = 0, while T d (Γ = 0) = 0.153 is finite. The model is thus described, for Γ = 0 and T ≤ T d , by the 1RSB equation with Parisi parameter m = 1.

Solutions of the cavity equations and spinodal lines
As discussed in details in [17,48], the central quantity that is computed by the cavity method is the free energy function φ(m), which is the Legendre transform of the complexity function, φ(m) = min f {f − T Σ(f )/m}. The thermodynamic free energy of the system is obtained by maximizing φ(m) in the interval m ∈ [0, 1], therefore in the following we restrict ourselves to these values of m.
When solving the cavity equation, we found that in addition to the trivial RS solution, in some regions of (T, Γ) there exist two different non-trivial 1RSB solutions. The first one develops continuously from the RS solution through a linear instability: hence the overlap q 1 grows continuously from its minimal value 1/q col , and we refer to this solution as the "low overlap" solution 1RSB lq . The second solution appears discontinuously, as in the classical coloring problem: the overlap q 1 is typically larger, hence we refer to this solution as the "high overlap" solution 1RSB hq . The coexistence of two different 1RSB solutions is quite unusual, but had been observed before: a good pedagogical example is discussed in details in Appendix A.
We find that the 1RSB lq solution exists in an interval m ∈ [0, m sp lq ], while the 1RSB hq solution exists in an interval m ∈ [m sp hq , 1], as illustrated in Fig. 4. Both m sp lq and m sp hq depend on T and Γ, to lighten the notations we keep implicit these dependencies. To delimit the region of existence of these solutions in the (T, Γ) phase diagram, we define the following transition lines: • The clustering (or dynamic) transition line T d (Γ) is defined by the condition m sp hq = 1. It is the continuation of the classical clustering transition T d (Γ = 0) to the quantum regime. Above this line, the 1RSB hq solution only exists outside the interval m ∈ [0, 1] and is therefore irrelevant for the thermodynamics of the model.
. In this sense the line T ⋆ (Γ) corresponds to a "critical point": below this line, there is a first order transition in m between the two non-trivial solutions; on the line, the first order transition disappears, a transition still exists and is of second order; above the line, there is a unique analytic 1RSB solution that we still call 1RSB lq because it connects continuously to the RS one.
As illustrated in Fig. 4, these three lines divide the (T, Γ) plane in four different regions. At high enough temperature, T > max(T i , T d ), only the RS solution exists. In the region below T d and above T i , the RS solution remains stable but it coexists with the 1RSB hq one (right, lower panel). In the region below min(T d , T i , T ⋆ ), the two non-trivial 1RSB solution coexist (right, middle panel). In the region below T i and above max(T d , T ⋆ ), only the 1RSB lq exists, either because the 1RSB hq has moved outside the interval m ∈ [0, 1] at T d , or because the two solutions have merged at T ⋆ (right, top panel). Note that in the numerical computation each solution, for the same point of the phase diagram, can be reached by following different routes: starting from high Γ and decreasing it to reach the 1RSB lq solution, or starting from Γ = 0 and increasing it to reach the 1RSB hq solution.

Thermodynamic phase diagram
The thermodynamic free energy of the problem corresponds to the maximum of φ(m), which is reached in m = m s (T, Γ). We have seen that in some regions of the phase diagram several solutions for φ(m) can coexist, leading to different branches of this function. This leads to several thermodynamic phase transitions that we now describe. For clarity, in Fig. 5 we report in the left panel the thermodynamic lines only, while in the right panel we report the complete phase diagram including the spinodal lines that were already discussed in Fig. 4.
Let us recall the following properties of the 1RSB solution. First of all, φ(m = 1) is always equal to the replicasymmetric free energy. Second, the complexity Σ(m) = βm 2 φ ′ (m) is related to the first derivative of φ(m). We find the following phase transitions in the model: • At high temperature the system is in a RS paramagnetic (P) phase. Upon lowering the temperature at low enough Γ, the dynamical transition T d (Γ) is met and below it the 1RSB hq appears. As in the classical limit Γ = 0, for Γ small enough the new solution is such that Σ(m = 1) ∝ φ ′ (m = 1) > 0. By convexity, for any m < 1 the free energy of the 1RSB hq solution is smaller than the one of the RS solution and the system remains in the paramagnetic phase. However, as in the classical case the paramagnetic phase is not ergodic because it is a superposition of many metastable states, and we refer to it as "dynamical paramagnet" or dP. • Upon lowering further the temperature, the complexity Σ(m = 1) decreases and reaches zero at the condensation temperature T K (Γ). At this point the slope of φ hq (m) in 1 becomes negative, in such a way that φ hq (m) has a maximum in the interval m ∈ [0, 1], and its value at the maximum is necessarily larger than the RS solution ( Fig. 4, right, lower panel). Hence at T K a thermodynamic transition happens between the RS and 1RSB hq phase. This transition is a standard RFOT transition: it is of second order thermodynamically, because the maximum of φ hq (m) moves smoothly enough away from m = 1, but it is of first order from the point of view of the order parameter because the system is jumping discontinuously from the RS to the 1RSB hq phase. We call the spin glass phase below T K the SG hq phase. Note that we find numerically that the line T K (Γ) grows linearly with Γ at small Γ: this is a consequence of the clustered structure of the classical spin glass phase, as will be discussed in Sec. IV A.
• If the temperature is lowered below T i (Γ) at high enough Γ, the RS phase becomes linearly unstable before the line T K is met. In this case, the possible presence of the 1RSB hq solution is irrelevant as we already discussed.
On the other hand, the solution 1RSB lq appears at T i (Γ) and we always find that it appears with a negative complexity at m = 1, Σ(m = 1) < 0, which implies that the slope of φ lq (m) in m = 1 is negative. The maximum of φ(m) is thus assumed on the 1RSB lq in this case (Fig. 4, right, upper panel). The transition at T i (Γ) is a thermodynamical third order transition (see Appendix A), and the overlap grows continuously from its replica-symmetric value. We call the phase below T i the SG lq phase.
• Finally we should discuss the behavior in the spin glass region, for T < min(T K , T i ). The line T K (Γ) can be continued in this region, and above it φ ′ hq (m = 1) > 0, hence the 1RSB hq is surely metastable with respect to the 1RSB lq solution. However, below this line, φ ′ hq (m = 1) < 0, hence both the 1RSB lq and 1RSB hq solutions have a maximum for m ∈ [0, 1]. This can lead to a thermodynamical first order transition between these two solutions if the values of φ(m) at the two maxima cross. This must happen at some temperature T fo (Γ) < T K (Γ), because we know that the 1RSB lq solution is stable at high temperature while the 1RSB hq is stable at low temperatures. Unfortunately, the free energy differences are so small in this region that our numerical accuracy does not allow us to determine the first order transition line. The line drawn in Fig. 5 is therefore schematic. We can only say that the variation of the complexity at m = 1 is found to be much higher in the 1RSB hq solution than in the 1RSB lq one. This is actually related to the fact that T K is a second order transition while T i is a third order transition. It implies in particular that T fo must be tangent to T K when they separate, and it also suggests that T fo (Γ) must be always very close to T K (Γ). Note also that the line T fo (Γ) must necessarily end on the line T ⋆ (Γ) where the distinction between the two 1RSB phases disappears.
As a final remark, we note that within our numerical accuracy it seems that the three lines T d (Γ), T K (Γ) and T ⋆ (Γ) cross at a single point. This finding is perfectly compatible with all of our data. To better understand this point let us call (T 0 , Γ 0 ) the point where the lines T d (Γ) and T ⋆ (Γ) cross. If T K (Γ) does not cross the lines in the same point, there are two possibilities: (i) T K crosses T ⋆ at some T < T 0 . However, this can be excluded because at T ⋆ the two 1RSB solutions must merge and convexity implies that the slope of φ hq in m = 1 must be negative. (ii) T K crosses T d at Γ < Γ 0 . This scenario seems to us logically possible. However one must stress that in most cases the lines T K and T d do not cross, but merge into some kind of critical point. Therefore the scenario of a triple crossing at (T 0 , Γ 0 ) seems likely. Unfortunately, we could not devise a more solid argument.
In summary, the thermodynamic phase diagram contains a paramagnetic (P) phase, a dynamical paramagnet (dP), a "high overlap" spin glass (SG hq ) and a "low overlap" spin glass (SG lq ). The transition at T d from the P to the dP phase is a standard dynamical (clustering) transition, and has no thermodynamic consequences. The transition T K from dP to SG hq is also a standard RFOT transition, i.e. a second order transition. The transition T i from P to SG lq is a third order thermodynamic transition. Finally, a first order thermodynamical transition between SG lq and SG hq must exist, even if our numerical accuracy is not sufficient to determine it precisely.
In particular we find that in the limit of zero temperature, the system becomes a spin glass under the action of an infinitesimal transverse field. At very low temperatures the two spin glass phases transform smoothly into each other, and the only phase transition is the third order transition at Γ + i (T → 0) ≃ 0.47. The artifact of the spurious first order transition predicted by the RS computation, discussed for c = 13 at the end of Sec. III A and in Appendix C 3, occurs also in the case c = 9, even if at even lower temperatures with respect to c = 13.

IV. STRUCTURE OF THE SPIN GLASS PHASE
A. Clusters crossings in the spin glass phase As sketched in Sec. II D 2, apart from the various transition lines that appear for this model when the quantum field is turned on, we also expect the spin glass phase to exhibit an interesting quantum behaviour due to its rich classical structure. The intuition based on the study of the quantum random subcubes model [17,32] recalled in Sec. II D 2 is that a continuum of (avoided) level crossings occurs in the spin-glass phase, because the classical model admits a pure state decomposition with clusters of various energies and entropies. As clusters with a larger classical energy have also a larger entropy, their total (quantum) energy is reduced faster when Γ increases, hence the crossings which lead to exponentially small gaps because of the extensive Hamming distance between clusters. The linear behavior of the line T K (Γ) reported in Sec. III B was a first confirmation of this mechanism, we report here a further evidence in its favour (the following discussion was already partially undertaken in [17]).
The coloring problem for q col = 4, c = 9 is in its clustered phase at zero temperature and zero field.
The key observation is that, as in the quantum random subcubes model [32], we expect m x (A, Γ) to be finite when Therefore, the energy of the ground state of a cluster is linear at small Γ, e(Γ) ∼ −Γm 0 x (A). Largest clusters yield the greatest decrease in energy when quantum fluctuations are switched on, and they dominate at zero temperature as soon as Γ > 0. Because these are the states with maximal entropy, they correspond to Σ(s max ) = 0. Hence as soon as Γ > 0, the zero temperature complexity abruptly drops to zero. In other words, we expect the system to condense into the largest clusters under an infinitesimal amount of quantum fluctuations. This in particular implies that, as in the quantum random subcubes model, a non-zero T K (Γ) should emerge, and that T K (Γ) ∝ Γ for small Γ, the effect of the transverse field entering the computation of the free energy under the combination βΓ. This feature is confirmed by the phase diagram presented on Let us now probe the hypothesis m 0 x (A) > 0, first numerically. Our path integral formalism does not allow to work directly at zero temperature, hence we had to run several Monte-Carlo simulations at low temperature and extrapolate the results to the limit T → 0. In this case, because the system is in a dynamical 1RSB phase at Γ = 0, it is possible to use "quiet planting" [67] to construct configurations equilibrated at T = 0. We started all the simulations in the same quietly planted configuration at Γ = 0, for the same instance of the coloring problem; in this way we assumed that the simulations all follow the evolution of the same cluster. Extrapolating the results to T = 0 gives the ground state properties of the cluster as m 0 x (A) = lim Γ→0 lim T →0 m x (A, T, Γ), as shown on Fig. 6. Note that at the very low temperatures we investigated, and given the size of the graph, the planted solution at T = 0 is actually an equilibrium configuration of the model at Γ = 0 for all temperatures of Fig. 6.
It is also possible to get a theoretical understanding of the existence of this finite limit for lim Γ→0 lim T →0 m x (A, T, Γ). First of all in the quantum random subcubes model, one can derive the exact expression m x (A, T, Γ) = s(A) ln 2 tanh(Γ/T ). Hence m 0 x (A) = s(A) ln 2 , and a perturbation theory in Γ to compute m x (A, T, Γ) will be divergent when T → 0, the coefficient of the term linear in Γ being proportional to m x (A)/T . For the coloring problem, no closed expression can be derived for m x (T, Γ), but it is still possible to compute it perturbatively in Γ within the path-integral representation of the partition function. The details of the computation are given in Appendix D. The important conclusion is that in the limit of small Γ and T (with the limit Γ → 0 taken before T → 0) one gets the asymptotic where σ is an uniformly chosen groundstate of the cluster A, and q auth (σ, i) denotes the number of colors the site i can take without creating a monochromatic edge if the rest of σ is kept fixed (in other words the number of colors that do not appear in the neighbors of i in σ). The expansion with Γ → 0 before T → 0 is thus the same as for the quantum random subcubes model. Altough we cannot prove that m 0 x (i.e. with the order of the limits reversed) is again proportional to q auth − 1 A in this case, we expect that it will still be positively correlated with it, and with the internal entropy of the cluster. This could be checked numerically by repeating the above computation for many different clusters of different entropies. In any case the demonstration given above of a non-zero m 0 x (A) is an additional piece of evidence in favour of the mechanism of crossings induced by the competition between the energy and the entropy of the classical pure states.

B. Quantum Monte Carlo annealing
Let us conclude this section by pointing out the consequences of this clustered structure of the spin glass phase for a quantum annealing of the problem. One expects that, because of this complicated structure, a quantum annealing of the problem starting from high Γ down to Γ = 0 will remain stuck in large clusters of higher energy when Γ gets small, and will not be able to find the smaller clusters that contain the ground states. However, it is not possible to simulate the real-time Schrödinger quantum dynamics of such a diluted model for reasonable values of N (the size of the Hilbert space is equal to 4 N ). Still, the effects of clustering can be seen on a classical alternative to quantum annealing, namely the annealing of a path integral Monte-Carlo computation (PIMC), whose results are shown on Fig. 7. These simulations were performed in two ways: in the first run, we prepared a typical graph together with one of its typical configurations at a very small temperature T = 0.06 via the quiet planting technique [67]; we initialized the PIMC at Γ = 0 in this configuration (which, given the size of the graph, is actually at zero energy), and then slowly increased Γ. Note that the results of this run perfectly match the results of the 1RSB quantum cavity computation, confirming the validity of our analysis. In the second run, we initialized and equilibrated the PIMC at Γ = 2 (in the paramagnetic phase), and then decreased Γ down to Γ = 0. On the scale of the figure, no difference between the two PIMC runs is observed, however a closer look (inset of left panel in Fig. 7) reveals that decreasing Γ one obtains a positive residual energy at Γ = 0. The latter is found to be larger than the residual energy after an infinitely slow thermal annealing, which is obtained by preparing a quietly planted configuration at T d and performing a slow classical annealing down to T = 0 [67,68]. Although a direct comparison is not possible because the PIMC annealing has not been extrapolated to the infinitely slow limit (we used steps of transverse field dΓ = 0.02 with 1000 sweeps per step), this result suggests that an annealing of a PIMC simulation is not more efficient than a thermal annealing for this model. We believe that this is once again related to entropic level crossings inside the spin glass phase, as in the quantum random subcubes model. Note that the residual energy of the PIMC annealing is a priori not directly related to the one of a true quantum annealing.
Similar results are obtained for connectivity c = 13, however in this case the quiet planting technique cannot be used at T = 0 because T K > 0. Therefore we used as a starting point for the increasing Γ run a classical configuration obtained through a slow thermal annealing from T d . These results are reported on Fig. 8.

V. CONCLUSIONS
In this paper we presented a study of the classical coloring problem, to which a quantum transverse field is added to represent the action of a quantum computer performing a quantum annealing with the aim of finding solutions to the classical problem. The study of this particular problem was motivated by the fact that among classical random optimization problems, it is the simplest one that shows an exponential degeneracy of solutions inside clusters (in the jargon of [64], it is a non-locked problem). Therefore, this is the simplest non-trivial model where the predictions of [32], that were obtained on a toy model, could be tested. We believe that similar conclusions would be reached on similar problems, such as random K-SAT.
Our main results are the characterization of the phase diagram in the (T, Γ) plane of parameters, and a further description of the spin-glass phase at low temperature and transverse field. Concerning the latter issue, we found evidences that the mechanism described in [32] is at work also in the case of the coloring, more precisely that (i) at low Γ, the line T K (Γ) is linear in Γ for c = 9, with T K (Γ = 0) = 0 and T d (Γ = 0) > 0 (see Fig. 5); and that (ii) the transverse magnetization of the ground state of a cluster goes to a positive constant in the limit Γ → 0 (see Fig. 6). Both these results are direct consequences of the exponential degeneracy of solutions in the clusters [32], and they show that the general scenario proposed in [32] using a toy model is at work in realistic random optimization problems.
A new feature of the present study, which could not be expected based on the analogy with the random subcubes model, is that the quantum phase transition between the spin glass and quantum paramagnetic phases at low temperature is of third order thermodynamically. Although we do not know the scaling of the gap at such a transition, the analogy with the results of [25] leads us to believe that the gap will be polynomial right at the quantum phase transition. On the other hand, based on the results of [32], we expect that the spin glass phase is characterized by a continuum of level crossings with an everywhere exponentially small gap. It would be extremely interesting to study numerically the gap in given instances of this model; exact diagonalization is impossible due to the rapid growth of the Hilbert space with N (as 4 N ), however this might be doable through Quantum Monte Carlo following [25,69].

Acknowledgments
We warmly thank Laura Foini and Florent Krzakala for useful discussions related to this work. Numerical computations were performed in part at the MesoPSL computing center, with support from Région Ile de France and ANR, in part using HPC resources from GENCI-CCRT/TGCC (Grant 2012056924), and in part using local computational resources that were bought thanks to the PIR grant of ENS "Optimization in complex system".
Appendix A: A reminder on the fully-connected p-spin spherical model in a longitudinal field

Definition and phase diagram
During the study of the quantum version of the coloring problem we have encountered several phenomena (in particular a third order continuous transition towards a spin glass phase, and multiple RSB solutions) that also appear in a much simpler, classical disordered model, namely the fully-connected p-spin spherical model in a field [65]. Analytical computations, and in particular expansions around the transition line, can be performed explicitly in this model and have constituted a very useful guideline for the analysis of the quantum coloring model. For these reasons we briefly collect in this Appendix the main results on this well-known model that are enlightening with this application in mind and refer the reader to [65,66,[70][71][72][73] for more details on this model.

Its Hamiltonian is
where the quenched couplings J i1,...,ip are independent Gaussian random variables of zero mean and variance p!/2N p−1 , and the degrees of freedom σ i are continuous real variables subject to the spherical constraint N i=1 σ 2 i = N . The partition function is defined as the integration of the Gibbs-Boltzmann factor e −βH over the N -dimensional sphere of radius √ N , and the associated free energy density concentrates in the thermodynamic limit around its quenched average. The latter can be computed with the help of the replica method and reads where the variational function is In this expression m is the Parisi parameter discussed in the main text, while q 0 (resp. q 1 ) is the typical overlap between two configurations in different (resp. in the same) pure states. Note that the 1RSB potential φ(m) discussed in the main text corresponds to the maximization of this function with respect to q 1 and q 0 , with m fixed. The phase diagram of this model, shown on Fig. 9, is obtained by solving the maximization problem defined in (A2,A3); depending on the values of β, h the supremum of φ 1rsb is either reached in the subspace of parameters q 0 = q 1 (this corresponds to a replica symmetric situation), or on a non-trivial 1RSB solution with q 0 < q 1 . The phase transition that separates these two regimes (and that reveals itself as a singularity in φ(β, h)) as a function of the temperature, changes qualitatively depending on the value of h. For h smaller than h c the transition is of the RFOT type, as described in Sec. II A: there exists a dynamical transition temperature T d (h) and a condensation (Kauzmann) transition temperature T K (h) < T d (h), the free energy undergoing a thermodynamic phase transition at the latter. On this line T K (h) the overlap order parameter is discontinuous, yet the value of the static Parisi parameter m s , i.e. the one maximizing (A3), is m s (h, T K (h)) = 1, which makes this discontinuous transition second order from a thermodynamic point of view. On the contrary for h > h c (but of course not too large) there is a single transition temperature T i (h), at which the order parameter of the RSB phase grows continuously with m s (h, T i (h)) < 1, and which is thermodynamically of third order. The three lines T d , T K and T i meet at h = h c . In the following we give some details on the derivation of these properties of the phase diagram, and also on some features that are not directly relevant for its thermodynamic properties (a spurious first order RS transition and the coexistence of different 1RSB solutions at m = m s ).

RS spinodals
Ti  Tc). The orange dash-dotted line is the limit of coexistence of the two RS solutions, whose free energies cross at a spurious first order phase transition not shown on this figure.

The RS solution(s) and its stability limit
Let us define the replica-symmetric variational free energy, function of a single overlap q, by substituting q 0 = q 1 = q in (A3); one then sees that this expression is independent of m, and is equal to The stationary points of this function of q are solutions of the following equation, Depending on the values of the parameters β, h there exists either one or three solutions to the RS equation (A5). In the latter case the intermediate one corresponds to a minimum of φ rs and can be discarded, while the two extreme ones compete to maximize φ rs and thus cross at a first order transition line. The boundary of the domain of existence of multiple RS solutions in the (h, T ) plane is more easily expressed parametrically, as a function of q; on this boundary one has the relation (A5) and in addition The orange dash-dotted line in Fig. 9 represents this spinodal limit of existence of multiple RS solutions. The first order transition between these two solutions occurs on a line, not shown on the figure, that starts at some T > 0, h = 0 and joins the spinodals at their cusp. We shall see that this transition is irrelevant thermodynamically, as it occurs inside the 1RSB phase.
It is easy to check that if q is a stationay point of φ rs , then the first derivatives of φ 1rsb with respect to q 0 , q 1 and m vanish in (q 0 , q 1 , m) = (q, q, m), for all values of m, i.e. that (q, q, m) is a stationary point of φ 1rsb . Let us discuss more precisely its nature, assuming q is a local maximum of φ rs . All second derivatives of φ 1rsb with respect to q 0 , q 1 , m which involve at least one derivative with respect to m vanish, hence one can concentrate on the Hessian of φ 1rsb with respect to q 0 and q 1 , and in particular on its determinant. A short computation leads to det Hess(m, q) = As m ∈ [0, 1] and ∂ 2 φ rs /∂q 2 ≤ 0 on a maximum of the RS potential, the sign of the determinant of the Hessian can only change when the last factor of this equation vanishes, i.e. when The black line T i (h) (or, alternatively, the two branches h ± i (T )) in Fig. 9 corresponds to the parametric representation, as a function of q, of the conjoint solution of (A8) and of the RS equation (A5). It starts from (T = 0, h − i (T = 0) = 0) (corresponding to q = 0), reaches a maximum in T c = 2 (p−1)(p−2) p−2 hitting the zero temperature line again at h + i (T = 0) = p(p−2) 2 (for q = 1). Outside the region enclosed by T i (h), i.e. at high temperature or field, the determinant of the Hessian is positive, corresponding to two negative eigenvalues. A local maximum of φ rs then corresponds to a local maximum of φ 1rsb in the larger RSB subspace of parameters, and such a RS solution is at least locally correct. On the contrary once T i (h) is crossed one of the eigenvalues of the Hessian becomes positive, there is one direction that increases φ 1rsb starting from the RS local maximum, which is thus locally unstable (for all values of the Parisi parameter in [0, 1]). We will come back in more details in Appendix A 4 on the behavior of the model around the instability line T i (h). Let us also mention as a technical detail the issue of the stability of the RS solutions in their domain of coexistence. The latter is traversed by the small field branch of the instability line, h − i (T ); in the intersection of the domain of coexistence of two RS solutions with the interior of the instability domain T < T i (h) both RS solutions are locally unstable. When two RS solutions coexist with T > T i (h) the one continuously connected to the high temperature regime is locally stable, while the other one is locally unstable. In any case even the locally stable solution is irrelevant as the non-trivial 1RSB solution will have a larger free energy, as discussed below.

The RFOT-like transitions at small h
We have seen above one mechanism for the appearance of replica-symmetry breaking, namely the transformation of a local maximum q of φ rs from a local maximum of φ 1rsb with q 0 = q 1 = q to a saddle. Another possibility is the appearance of a local maximum of φ 1rsb with q 0 = q 1 , i.e. far away from the replica-symmetric subspace (hence the "discontinuous" character of such a transition). At variance with the local instability, which is independent of the value of m, the occurence of the discontinuous transition is m-dependent. We shall start by considering the most important case m = 1, and come back on this m-dependence in Appendix A 5.
At m = 1 the stationary conditions of φ 1rsb , ∂φ 1rsb ∂q0 = ∂φ 1rsb ∂q1 = 0 imply that q 0 is solution of the RS equation (A5), while q 1 ≥ q 0 verifies: Only for some range of parameters the equation f (x) = f (q 0 ) admits a non-trivial solution q 1 > q 0 . On its limit of existence one has the additional condition f ′ (q 1 ) = 0, i.e.: The line T d (h), corresponding to this limit of existence and drawn in dashed orange on Fig. 9, can be obtained parametrically: from (A10) one obtains β as a function of q 1 , then q 0 is expressed in terms of q 1 as a solution of f (q 0 ) = f (q 1 ), and finally the field is obtained from (A5), q 0 being solution of this RS equation. The relevant interval for this parametrization is q 1 ∈ [ p−2 p , p−2 p−1 ]. At the lower limit the line T d merges with the local instability T i at h c ; in this limit the two extrema of f merge, q 1 → q 0 and (A10) indeed reduces to (A8). The upper limit q 1 = p−2 p−1 corresponds to h = 0, q 0 = 0, and one can find the explicit expression T d (h = 0) = p(p−2) p−2 2(p−1) p−1 . The non-trivial 1RSB solution just found at m = 1 might be the global maximum of φ 1rsb , or not. To test this one has to compute the derivative of φ 1rsb with respect to m, in this point. As discussed in Sec. II A this is proportional to the complexity; when the latter is positive m = 1 is indeed a maximum of φ 1rsb , otherwise there will be a maximum at a value m s < 1. The Kauzmann, or condensation, transition separates these two regimes, and corresponds to the line drawn in solid blue on Fig. 9.

The third order transition at large h
We now come back in more details on the behavior of the free energy around the RS local instability. We want to justify that the perturbations in the overlaps q 0 and q 1 with respect to their common RS value are linear in the distance from the instability line, and that the first discontinuities appear in the third derivatives of the free energy φ(β, h) (hence the transition is thermodynamically of third order). In particular, the complexity takes a scaling form in (δT + δh) 3 when one approaches the instability line (see Fig. 10). These last two facts have been used in the analysis of the quantum coloring problem.
As explained in Appendix A 2, when crossing the instability line T i (h) the RS extremum goes from a maximum to a saddle (in the larger 1RSB space of parameters). By definition the determinant of the Hessian matrix written in Eq. (A7) vanishes in (h, T i (h)). More precisely, one can diagonalize the matrix in this point and find that it has one eigenvalue 0 with eigenvector −(1 − m) m and one strictly negative eigenvalue associated to the eigenvector Consider now some parameters (β, h) inside the unstable region, denote q rs (β, h) the associated solution of the RS equation, and expand the 1RSB potential around this point in powers of δ + , δ − as As q rs is a RS solution, the expansion begins with second order terms. In the following we shall only need the coefficients of the terms of order δ 2 + and δ 3 + , denoted a/2 and b/3 respectively, with Consider now that (β, h) is at a small distance, call it ǫ, of the instability line; a short reasoning reveals that the coefficients of the terms in δ 2 + and δ + δ − are of order ǫ, while all others are of order 1 in this limit (in fact, note that the perturbation δ + is in the direction of the vanishing eigenvalue of the Hessian). One then sees that the maximization of δφ over δ + and δ − will lead to δ + = O(ǫ), δ − = O(ǫ 2 ), and that at the lowest order δφ is simply given by the maximization of aδ 2 + /2 + bδ 3 + /3. We thus obtain in this limit where we made some simplifications using Eq. (A8) verified by q, the value of the overlap at the point of the instability line approached in this limit. The term in square brackets is positive in the RS unstable region, and vanishes linearly on the instability line. As φ rs (β, h) is regular across this line this demonstrates the third order character of the transition. The dependence on m of the 1RSB potential can then be easily studied, in particular its maximum is reached in m s (q) = (p−2)(1−q) 2q . This behaviour is checked numerically by the scaling plots of the complexity, see Fig. 10. We emphasize that as δ + = O(ǫ) on the maximum of φ, the overlap order parameter q 1 − q 0 grows linearly away from the instability.
It could seem at this point that the above study is valid along the two branches h ± i (T ); this is however not the case, as further considerations reveal. The maximization over the 1RSB overlaps must enforce the condition q 1 ≥ q 0 , which translates here into δ + ≥ 0. For the local maximum of aδ 2 + /2 + bδ 3 + /3 to happen on this side of the origin one must have b < 0 (we always have a > 0 in the unstable regime). Assuming m ∈ [0, 1], and after some simplifications, one sees this condition to be equivalent to m s (q) + m(1 − 2 m s (q)) > 0. Two cases can now be distinguished: • The high-field branch h + i (T ) corresponds, in the q-parametrized representation, to q ∈ p−2 p , 1 . Then m s (q) ∈ [0, 1], thus the condition q 1 ≥ q 0 at the local maximum is respected for all m ∈ [0, 1], and the maximum of φ(m) is reached in an acceptable value m s (q) ∈ [0, 1] of the Parisi parameter.
• On the contrary for the low-field part of the instability line h − i (T ), m s (q) > 1. The condition q 1 ≥ q 0 is only fulfilled for m ∈ [0, m sp lq ], where m sp lq = ms(q) 2 ms(q)−1 . For these values of m there exists a non-trivial solution of the 1RSB equations that grows continuously from the RS one. But for m ∈ [m sp lq , 1] one sees that the maximum of the expansion would correspond to δ + → ∞, in other words the maximum of the 1RSB potential is far from the RS solution, and corresponds to the discontinuous solution.

Coexistence of RSB solutions
We have just seen that close to the low-field branch h  Fig. 11). However, once one also extremizes over m, no sign of the continuous solution to the 1RSB equation remains, and the physical solution is always the discontinuous one. To justify this point it is convenient to consider the reduced potential: where the dependencies on T and h are kept understood, and the normalization has been chosen such that V vanishes on the replica-symmetric solution. A plot of V as a function of q 1 for h slightly larger than h − i (T ) and several values of m is shown on the right panel of Fig. 11. It can be seen that, as expected from the previous discussion, V has for m small enough (m < m sp lq (T, h)) a local maximum near the replica-symmetric point q 1 = q 0 (indicated by a tick on the figure), corresponding to the continuous solution. On the other hand, there exists for m sufficiently large (m ≥ m sp hq (T, h)) another local maximum of V at a value of q 1 which is further away from q 0 . This second maximum corresponds to the discontinuous solution. When m → m sp lq (T, h), the continuous solution merges onto the discontinuous one. Increasing h, m sp lq (T, h) shrinks and the limit of coexistence of both solutions (at fixed m) is given by the critical field h ⋆ (T ) such that m sp lq (T, h ⋆ (T )) = m sp hq (T, h ⋆ (T )). In the region h ∈ [h − i (T ), h ⋆ (T )] thermodynamical observables acquire two branches when computed at fixed m, as shown for the complexity on the left panel of Fig. 11. Finally, the relevant solution from a thermodynamical point of view is the one which maximizes V (m, q 1 ) over both m and q 1 ; because the continuous maximum is always very small and merges with the discontinuous one upon increasing m, one finds that the relevant solution is always the discontinuous one. by this change is the replacement of the functions g, z by (see [17] for more details) z(η 1 , . . . , η c−1 ) being defined by normalization, and n(σ) counting the number of discontinuities in the trajectory σ.
Note that for Γ = 0, only constant trajectories remain and one recovers the definition given in (3).
The population dynamics method [50,74] is a convenient way to solve these equations numerically. Probability distributions are approximated by weighted samples of representative elements: We therefore have two deal with two level of populations: one of N pop messages η, each of these messages being represented by N traj trajectories σ. . This induces both systematic deviations of the empirical mean from the exact value and noise in its estimation; extrapolations to N pop , N traj → ∞ via finite size analysis can in principle be performed to reduce these effects, we show one example of such treatement in App. C 1. Further difficulties arise when the weights a i (or b i,j ) become very heterogeneous: then the population is effectively supported only on the representants with the largest weights, hence the effective size of the population can be much smaller than the number of samples. These population representations of probability distributions are known in statistics as particle approximations [75]; many resampling techniques are known to fight this impoverishment of the sample representativity, but there does not seem to be an universal way to avoid it.
To increase the speed of our numerical code we made use of the parallelization opportunities offered by multi-core computers. Let us sketch how the 1RSB equation resolution can be distributed on several processing units. The update procedure that we alluded to above can simply be thought of as a way to generate a sample at step τ + 1, {(η i (τ + 1), a i (τ + 1))}, from a sample at step τ , {(η i (τ ), a i (τ ))}. The important point is that the new representants (η i (τ +1), a i (τ +1)) are generated independently one from the other (apart from the normalization condition i a i = 1 that can easily be enforced). The update procedure can then be parallelized as follows: each core is first sent the whole sample of messages and weights {η i (τ ), a i (τ )} 1≤i≤Npop . Then each core c, with 1 ≤ c ≤ N core , generates independently a new sample of N pop /N core messages: {(η c,i (τ +1), a c,i (τ +1))} 1≤i≤Npop/Ncore . These new messages are then gathered together to form the full sample at step τ +1: {(η i (τ ), a i (τ ))} 1≤i≤Npop . This method does not change the memory limit of the procedure, because each core is sent the whole population in the first step, but allows for a gain in time roughly proportional to N core (the communication between the processors usually takes much less time than the generation of the samples). Typically we used N core between 12 and 64; values of N pop and N traj are limited by the amount of memory available on a single core, leading in our case to the constraint N pop × N traj 3.10 7 /(βΓ). A moment of thought reveals that it is also possible to avoid the step of population gathering, and to keep the population of messages {η i (τ ), a i (τ )} 1≤i≤Npop spread over the N core cores. This allows for larger populations, but strongly increases the time needed to update the population, as much more information has to be exchanged between cores; therefore, we dit not use this second procedure in this work. the Parisi replica-symmetry breaking parameter m. It is thus possible to take m = 0, which is very convenient in practice: the weights z({η i }) m are all equal to 1 (at non-zero temperature the factors z never vanish), the external population representation has thus homogeneous weights and N pop can be reduced to a rather small value without too much negative effects on the numerical accuracy. We can then concentrate on the finite population size effect of N traj and build a scaling theory around the transition. To simplify the discussion we focus on the case in which the temperature T is fixed and Γ is lowered from the quantum paramagnetic phase at large Γ through the first continuous instability Γ + i (T ). First of all, the RS solution has a very simple overlap structure: because of the symmetry between colors, one has q 0 = q 1 = 1 q col on the RS solution. Introducing q 1 = q 1 − 1/q col , one finds that within the RS phase (Γ > Γ + i (T )), q 1 has a finite population size behaviour as 1/ N traj . On the other hand, we expect from the discussion of Appendix A that for Γ < Γ + i (T ), q 1 = −a Γ − Γ + i (T ) . These two behaviours can be matched by introducing a scaling function F such that: Moreover, F (x) admits the following asymptotic behaviours: lim x→−∞ F (x)/x = −a, lim x→∞ F (x) = O(1). The value of Γ − Γ + i (T ) is then determined in order to obtain the best collapse of the numerical data for different values of N traj ; as shown in Fig. 12 this scaling form gives a very good collapse of the numerical curves. This procedure therefore allows one to determine the line T i (Γ) with very good precision. Having discussed the determination of T i (Γ), let us now focus on the more complex determination of T K (Γ) and T ⋆ (Γ) that were defined in Sec. III B and appeared in Fig. 4 and 5. We assume here that the reader is already familiar with the discussion of Sec. III B but we will repeat some parts of the discussion for clarity.
For simplicity, we will discuss the procedure at a fixed low temperature T = 0.09, and explain how Γ K and Γ ⋆ are determined at this temperature. At Γ = 0, there is a single 1RSB solution of the classical cavity equation, the 1RSB hq one, that exists in an interval m ∈ [m sp hq , 1]. This solution can easily be found by solving the classical 1RSB equations, it is the relevant one and has positive complexity at m = 1, because T d (Γ = 0) > 0.09 and T K (Γ = 0) = 0. The static Parisi parameter is therefore m s (Γ = 0) = 1. The evolution upon increasing Γ from Γ = 0 has already been sketched in Sec. III B and in Fig. 4, but we give some additional details here (see a summary on Fig. 13): • For Γ < Γ K (T ), in addition to the RS solution, there exists only the 1RSB hq solution, which is defined for m ∈ [m sp hq (Γ), 1] and has a positive complexity at m = 1. It can be easily constructed by initializing the 1RSB population dynamics in this solution at Γ = 0 and following it at positive Γ. The value of m sp hq (Γ) can be determined numerically by following the evolution of the 1RSB hq solution upon decreasing m from 1 to 0 at constant Γ, as the value at which the solution disappears and the complexity jumps to zero (see Fig. 14  The region of coexistence is therefore delimited by Γ − i < Γ < Γ⋆ and m sp hq (Γ) < m < m sp lq (Γ). However, only the 1RSB hq solution is thermodynamically relevant: the value of ms(Γ) is equal to 1 for Γ < ΓK, and it is smaller than 1 for Γ > ΓK, but it always correspond to the 1RSB hq branch. • Upon increasing Γ, the coexistence region shrinks; the field Γ ⋆ (T ) at which m sp hq = m sp lq marks the limit of coexistence of both solutions. However, when the coexistence region is too small it is hard to determine the spinodals, so the point Γ ⋆ is obtained by extrapolating the two lines m sp hq (Γ) and m sp lq (Γ) from smaller Γ to determine their intersection (see Fig. 13).
• Finally above Γ ⋆ there is a single 1RSB solution, and again m s (Γ) is the point where its complexity vanishes.
The evolution of the complexity in Fig. 14 is the best way to provide numerical support to the schematic picture for the free energy φ(m) outlined in the right panel of Fig. 4. Unfortunately, a direct computation of φ(m) is not possible because the data are too noisy. The fact that the relevant solution for thermodynamics is the 1RSB hq one is therefore not evident from Fig. 14. However, this can be deduced from the consistency arguments presented in Sec. III B and from the fact that the complexity of the 1RSB hq is much bigger, thus suggesting a more rapid rise of its free energy from the RS value.
The dynamical transition temperature T d (Γ) can be determined quite easily as the point where the 1RSB hq solution at m = 1 disappears on increasing temperature at fixed Γ, like in the classical case. An example of this procedure is shown on Fig. 15. It is important to stress that when the two lines T d (Γ) and T K (Γ) merge, the overlap q 1 does not tend to its RS value, as shown in Fig. 16. This is an important difference with respect to the spherical p-spin model, where the two lines merge on a critical point where the transition becomes continuous. In the coloring problem, the overlap q 1 remains non trivial, indicating that the lines T d and T K do not merge on the instability line T i , consistently with the phase diagram reported in Fig. 5.  Fig. 11. For this temperature, ΓK ≃ 0.064, Γ − i ≃ 0.08, Γ⋆ ≃ 0.25. The left plot is for ΓK < Γ < Γ − i : in this region there exists only the 1RSB hq solution with Σ(m = 1) < 0. For Γ − i < Γ < Γ⋆ (middle and right panel) there exist two branches 1RSB hq and 1RSB lq , which can be seen to merge when Γ approaches Γ⋆. To find the thermodynamic solution, one has to maximize the free energy of the system, which amounts to enforce the condition Σ(m) = 0 (or Σ(m = 1) > 0): this condition always selects the 1RSB hq solution. 3. The spurious RS first order transition for c = 13: a study in the (Γ, m) plane As mentioned above, at the replica-symmetric level and for T 0.12, the q col = 4, c = 13 model exhibits a first order phase transition when Γ is varied, while the transition is only third order when the 1RSB computation at the static value m s (T, Γ) of the Parisi parameter is done. Altough the 1RSB equation (4) yields thermodynamic observables only for m = m s (T, Γ), the static value of the Parisi parameter, one can also solve it for any value of m. This gives a way to interpolate between the first order phase transition at m = 1 and the continuous transition at m = m s (T, Γ). The results of cavity computations performed for various m are shown on the left panel of Fig. 17 for T = 0.06. The value m = 1 corresponds to the replica-symmetric calculation, in which the hysteresis between the solution coming from Γ = 0 and from large Γ is very well marked. Decreasing m, the hysteresis gets smaller and smaller, and cannot be seen anymore for m = 0.5. In this case m s (0.06, Γ i (0.06)) ≃ 0.13 and therefore the physical transition is continuous. As for the c = 9 model at small field, it is possible to build a scenario to explain these numerical results and to understand in greater details the nature of the transition; we briefly sketch it herafter.
These results can be summarized by drawing a phase diagram in the (Γ, m) plane, for T fixed (right panel of Fig. 17). At small Γ there exists only one 1RSB solution, defined for all values of m. However for larger values of Γ there exists a range of m (between the blue and red lines of Fig. 17, right panel) where two solutions of the 1RSB equation coexist, similarly to what we explained in Sec. III B 1 for q = 4, c = 9 model when Γ was decreased below Γ ⋆ . In this case, one finds that the thermodynamic solution at m s always lies on the "low q" branch, defined for m ∈ [0, m sp lq (T, Γ)]. Upon The high overlap solution exists above the red line with squares, and the low overlap one below the blue line with circles. These two lines meet at the value of m below which no hysteresis remain. The relevant solution is always the low overlap one; it becomes RSB for Γ below the continuous instability shown as a vertical dashed black line. Finally, we also report in this region the static value of the Parisi parameter ms(T, Γ).
increasing the field, m sp lq reaches 1, and above the value of the linear instabilty Γ i (T ), this low overlap branch becomes m-independent: this is the (quantum paramagnetic) replica-symmetric solution. On the other hand, the other "high q" branch still exists for Γ > Γ i not too large, and m ∈ [m sp hq (T, Γ), 1]. It completely disappears only when m sp hq (T, Γ) reaches 1, which corresponds to the limit of existence of the corresponding replica-symmetric solution. This branch is never relevant, because maximizing over m leads to a solution of negative complexity, that has to be dismissed. However, it explains why a replica-symmetric computation gave the illusion of a first order phase transition.
This expansion can of course be obtained from the usual formulas for the perturbation of eigenvalues, taking into account their possible degeneracies. We now assume that H Q takes the form of an uniform "transverse field" − i T i , with the flipping operators T i defined in Eq. (7). For simplicity of notation we shall work with the free energy density instead of the partition function. The expansion for this quantity reads where σ (i,σ ′ ) denotes the configuration obtained from σ by replacing its i-th variable by σ ′ , and we introduced the classical Gibbs-Boltzmann probability law µ(σ) = e −βE(σ) /Z(β, Γ = 0). Deriving with respect to Γ we obtain a similar equation for the transverse magnetization, We would like to emphasize here that these first quantum corrections are expressed as an average of simple, local, observables with respect to the classical Gibbs-Boltzmann measure -unfortunately, for higher order terms in Γ, the corresponding classical quantities involve correlation functions between variables at arbitrary distances. In particular, if the classical energy E(σ) is such that µ is correctly described by the classical cavity method (be it in its RS or RSB version), then the first quantum correction can be obtained by solving the classical cavity equations, which are much simpler than the quantum ones. As a matter of fact (D2) can be obtained from the quantum cavity method, either by expanding the quantum messages η(σ) in powers of the number of discontinuities they contain (0 for the classical trajectories, 2 at the lowest order in Γ) and noting that their first quantum correction can be expressed from the classical component of the incoming messages (see Eq. B1) or by taking the derivative with respect to Γ 2 in the variational expressions of the free energy provided by the cavity method in terms of its order parameter (η at the RS level, P (η) at the 1RSB one). In the 1RSB case the quantum corrections are thus the sum over the postulated pure states of the correction inside each pure state; this tells indirectly something about the classical cavity method, namely that most pairs of configurations belonging to two different pure states are at Hamming distance strictly larger than 1, as it should. As usual in perturbation theory, the presence or absence of close-by configurations with degenerate energies leads to qualitatively different behaviors; this is particularly relevant at low temperatures, as the two corresponding terms in Eqs. (D2,D3) have different scalings with β. In particular, for the coloring problem, the transverse magnetization given in Eq. (D3) satisfies where in the r.h.s. the average is over the ground states σ of the model, and q auth (σ, i) is the number of colors that i is authorized to take by its neighbors (in other words q col minus the number of distinct colors the neighbors of i take in σ): the more "floppy" the spins are, the stronger they respond to the quantum field. As mentioned above these small Γ expansions can be obtained from the quantum cavity formalism, at any level of replica-symmetry breaking. In particular at the 1RSB level the formula (D4) can be interpreted as an expansion for the transverse magnetization of a single cluster, provided the average in the r.h.s. is restricted to this particular classical pure state. This provides the