Disorder-free spin glass transitions and jamming in exactly solvable mean-field models

We construct and analyze a family of $M$-component vectorial spin systems which exhibit glass transitions and jamming within supercooled paramagnetic states without quenched disorder. Our system is defined on lattices with connectivity $c=\alpha M$ and becomes exactly solvable in the limit of large number of components $M \to \infty$. We consider generic $p$-body interactions between the vectorial Ising/continuous spins with linear/non-linear potentials. The existence of self-generated randomness is demonstrated by showing that the random energy model is recovered from a $M$-component ferromagnetic $p$-spin Ising model in $M \to \infty$ and $p \to \infty$ limit. In our systems the quenched disorder, if present, and the self-generated disorder act additively. Our theory provides a unified mean-field theoretical framework for glass transitions of rotational degree of freedoms such as orientation of molecules in glass forming liquids, color angles in continuous coloring of graphs and vector spins of geometrically frustrated magnets. The rotational glass transitions accompany various types of replica symmetry breaking. In the case of repulsive hardcore interactions in the spin space, continuous the criticality of the jamming or SAT/UNSTAT transition becomes the same as that of hardspheres.

We construct and analyze a family of M -component vectorial spin systems which exhibit glass transitions and jamming within supercooled paramagnetic states without quenched randomness. Our system is defined on tree-like lattices with connectivity c = αM and becomes exactly solvable in the limit of large number of components M → ∞. We consider generic p-body interactions between the vectorial Ising/continuous spins with linear/non-linear potentials. The existence of self-generated randomness is demonstrated by showing that the random energy model is recovered from a M -component ferromagnetic p-spin Ising model in M → ∞ and p → ∞ limit. In our systems the quenched randomness, if present, and the self-generated randomness act additively. Our theory provides a unified mean-field theoretical framework for glass transitions of rotational degree of freedoms such as orientation of molecules in glass forming liquids, color angles in continuous coloring of graphs and vector spins of geometrically frustrated magnets. The rotational glass transitions accompany various types of replica symmetry breaking. In the case of repulsive hardcore interactions in the spin space, continuous the criticality of the jamming or SAT/UNSTAT transition becomes the same as that of hardspheres.  Simple spin models often provide useful grounds to develop statistical mechanical approaches for various kinds of phase transitions. For the glass transition [1,2], which is one of the most important open problem in physics, a family of mean-field spinglass models called as p-spin spinglass models have played important roles [3][4][5][6]. The concepts and techniques used in the spinglass theory have promoted substantial progress of the first principle theory for the glass transitions of supercooled liquids [7,8]. Most notably exact mean-field theory in large dimensional limit was constructed recently for the hardspheres [9][10][11][12][13][14][15]. The color angle 0 < θ < 2π, as in the standard HSV color map, can be represented by a XY spin, i.e. a vector with M = 2 component (green arrow). The example shown here is a solution to the requirement that color angle on adjacent vertexes must be greater than or equal to 2π/3. b) Geometrically frustrated magnets: vectorial spins (green arrows) with antiferromagnetic couplings on adjacent vertexes on corner sharing triangles (e.g. kagome lattice), tetrahedra (e.g. pylochlore lattice). The ground states are highly degenerate due to the loose connectivity of the lattices. c) Glass forming liquids of molecules or colloidal particles with 'spins': a simple molecule or a colloidal particle, like the Janaus particle, is symmetric under rotation around an axis whose direction can be specified by a spin.
There remains, however, a conceptual problem regarding the origin of randomness. The spinglass models [16], which have been developed originally by Edwards and Anderson to model a class of disordered and frustrated magnetic materials [17], have quenched randomness which is apparently absent in glass forming liquids. It is often emphasized in the studies of spinglass materials that both the quenched randomness and frustration are important. However it is believed that somehow the randomness is self-generated in structural glasses which are born out of supercooled liquid and thus the quenched randomness is not necessarily. Early seminal works [18][19][20] have suggested that self-generated randomness are actually realized in some spin models without quenched randomness. However a comprehensive understanding of the mechanism of the putative self-generated randomness and its possible relation to the quenched randomness in spinglass models is still lacking.
In order to shed a light on this issue, we explicitly develop and analyze a family of mean-field vectorial spin models. We show that they exhibit glass transitions within their supercooled paramagnetic phases without quenched randomness. Our model consists of M -component vectorial spins, which can take either the Ising ±1 or continuous values, put on tree-like lattices with connectivity c = αM , which becomes exactly solvable in the limit of large number of components M → ∞. We perform a unified study of the crystalline phase (e.g. ferromagnetic phase), supercooled paramagnetic phases and glassy phases of the same model. We clarify the condition needed to ensure local stability of supercooled liquids and glasses against crystallization. We demonstrate in particular that the theoretical results of the random energy model [21] and the p-spin spinglass models can be fully recovered from a M -component p-spin models with purely ferromagnetic interactions within their supercooled paramagnetic phases. This proves the existence of the self-generated randomness in our models. In a sense this observation strengthen the view that the p-spin spinglass models are good caricature spin models for glass transitions [6] because the quenched randomness is actually not needed. We show that the quenched randomness, if present, add on top of the self-generated randomness.
Glass transition of the rotational or spin degrees of freedom is an important problem by itself and can be found not only in the spinglasses but also in many other real systems. It should be noted first that most of the molecules and colloidal particles in glass-forming liquids are not simply spherical but have rotational degrees of freedom because of their shape or patches on their surfaces (see Fig. 1c)) and the rotational degree of freedom can exhibit glass transitions simultaneously or separately from that of the translational degrees of freedom. Sometimes the rotational degrees of freedoms alone exhibit glassiness on top of crystalline long-ranged order of the translational degrees of freedom. This happens for instance in the so called plastic crystals where the rotations of molecules slow down and eventually exhibit glass transitions [22]. Another important problem is the spinglass transition found in frustrated magnets but without quenched randomness (Fig. 1b)) [23,24]. Possibilities of disorder-free spinglass transitions have been a matter of long debate in the field of frustrated magnets. We expect our results provide a useful basis to tackle these problems theoretically.
Within our formalism we consider p-body interactions through generic non-linear potentials. In particular we apply the scheme to the case of a M -component continuous spins interacting with each other through a hardcore potential which enables jamming transition of the vectorial spins. This is relevant in the continuous constrained satisfaction problems such as the coloring of graphs or periodic scheduling [25] (Fig. 1 a)): the problem is to put continuous colors parametrized by "color angle" 0 < θ < 2π on the vertexes of a given graph such that angles on adjacent vertexes are sufficiently separated from each other. This is exactly a continuous version of the usual coloring problem where one is allowed to use only discrete colors like red, green and blue [26,27]. Increasing the coordination number c of the graph, the solution space exhibit clustering transition (glass transition) and eventually SAT/UNSAT transition (jamming) above which one cannot find a solution which satisfies the constraints. Given the continuous variables, an interesting question is the universality class of the SAT/UNSAT transition. Closely following the analysis done on hardspheres in d → ∞ limit [11], we will show that the jamming criticality of our model belong indeed to the same universality of the hardspheres. Our result extends the result on the perceptron problem [28,29] which can be regarded as a special case p = 1 of our models.
The organization of this paper is as follows. In sec. II we introduce a family of large M -component vectorial Ising/continuous spin models with a generalized p-body interaction described by linear/non-linear potentials. We introduce a disorder-free model that has no quenched randomness and also a model which interpolates between the disorder-free model and a fully disordered spinglass model. In sec. III we discuss possible crystalline orderings in our disorder-free models and possibility to realize supercooled paramagnetic states, which are crucial as the basis for glass transitions to take place without the quenched randomness. In sec. IV we show that the random energy model can be recovered from a M -component p-spin Ising ferromagnetic model with a linear potential in the limit M → ∞ and p → ∞. This demonstrates the presence of self-generated randomness in our models. In sec. V we derive the replicated free-energy functional in terms of the glass and crystalline order parameters. We also discuss stability of the supercooled paramagnetic state and the glassy states against crystallization. In sec. VI we establish the connection between our model with linear potential and the standard p-spin spinglass models. Then in the subsequent sections, we turn to study glassy phases of our model with non-linear potentials limiting our selves to the case of continuous spins. In sec. VII we discuss some general results within the replica symmetric (RS) ansatz. In sec. VIII we discuss some general results within 1 step and continuous replica symmetry breaking (RSB) ansatz. In sec. IX we analyze the model with a quadratic potential as the simplest case of non-linear potential. In sec. X we analyze in detail the model with a hardcore potential which exhibit jamming. Finally in sec. XI we conclude this paper with some summary and remarks. Some technical details are reported in the appendices.
• Continuous spin M -component continuous spin with length |S| = √ M which can continuously rotate in the M -dimensional space. It is known in some models that this case is closely related to the 'spherical model' which has just M = 1 component spins S i normalized by a global constraint 30,31]. The spins are put on the vertexes of lattices (graphs) which are locally tree-like with no closed loops as shown in Fig. 2. Spins are involved in p-body interactions represented by factor nodes (interaction node) in the figure. Each spin is involved in c = αM p-tuples. Thus the number of the p-tuples is given by, In the present paper, we will limit ourselves with M → ∞ limit with fixed α, where we will ind that exact analysis of glass transitions becomes possible. We note that the thermodynamic limit N → ∞ limit must be took before M → ∞ limit.
The interaction between the spins is given by a generalized p-body interaction, where Here 1( ), 2( ), . . . , p( ) represent the spins involved in a given p-tuple . The function V (r) represents a generic interaction potential. We will call the argument variable r as 'gap', whose meaning will become clear later, with δ ∈ R being a control parameter.
In the present paper we mainly study models without quenched randomness (disorder-free model) but we also discuss models with quenched randomness (disordered model).
• disorder-free model Here ξ µ s are mutually independent, quenched random variables which obey the Gaussian distribution with zero mean and unit variance. The parameter λ represents the strength of the 'disorder-free' part in the disordered model. Note that the disorder-free model is recovered by choosing λ/ √ M = 1. In the other limit λ/ √ M = 0 we have completely disordered, spinglass model. Thus we have a smooth interpolation between the two limits with this parametrization.
The free-energy of the system can be written as, where β is the inverse temperature. The partition function Z is defined as, here Tr S represents a trace over the spin space of the spin S, where dS i is an integration over the surface of the M -dimensional sphere with diameter √ M . We have also introduced a Fourier transform of the Boltzmann's factor,

B. Linear and non-linear potentials
The most simple potential is the linear potential, This is a p-spin ferromagnetic model. We will use this potential in order to establish connections to the random energy model (sec. IV) and the p-spin spinglass models (sec. VI). As a simplest non-linear potential, we will consider briefly in sec. IX the quadratic potential, For the hardcore potential Eq. (219) the spin S0 in panel c) is excluded from the cones around each of the neighboring spins S1,S2,S3. (Note that, for instance, S2 and S4 can overlap if they are not directly connected by a link). The size of the cones grows with decreasing the parameter δ. Thus the excluding volume effect becomes larger by decreasing δ or increasing the connectivity c.
We will study in detail in sec. X the case of more strongly non-linear potential, The hardcore potential is obtained in → ∞ limit. This amount to bring in an excluded volume effect in the spin space similarly to the interaction between the hardspheres (See Fig. 2 c)). With p = 2 body interaction it can be used for the continuous coloring problem shown in Fig. 1 a)): spins representing the color angles on adjacent vertexes are forced to be separated in angle larger than cos −1 (δ/ √ M ) for the hardcore potential (See Fig. 2 c)). In the case of p = 1, and in the presence of quenched randomness ξ µ 's Eq. (88), the problem becomes the perceptron problem [32] [28]. The case p = 2 was also studied in part by a seminal work [33]. In the present paper we study the cases of p ≥ 2.

C. Pressure, distribution of gaps and isostaticity
With the soft/hardcore potential Eq. (13), the system becomes more constrained as we decrease the parameter δ much as an assembly of hardspheres becomes more constrained as the diameter of the spheres increase so that the volume fraction increases. This motivates us to introduce 'pressure' as an analogue of that in particulate systems, The normalization factor N is simply the number of interaction links in the system which is given by Eq. (1). Then it is also useful to introduce the distribution function of the gap, In the 1st equation . . . is the thermal average. In the 2nd equation δ/δ ln e −βV (r) is a functional derivative. Apparently the distribution function of the gap g(r) is analogous to the radial distribution function in the particulate systems. The pressure Eq. (14) can be rewritten using ∂(−βF )/∂δ = dr δ(−βF ) δ ln e −βV (r) (ln e −βV (r) ) and g(r) defined above as, Π = drg(r)(ln e −βV (r) ) = drg(r)(−βV (r)). (17) This is the analogue of the virial equation for the pressure in the liquid theory [34].
Given N spins S i (i = 1, 2, . . . , N ) with M components, which are normalized such that |S i | 2 = M , the total number of the degrees of freedom is N (M − 1). Each spin is involved in c = αM sets of p-body interactions (See Fig.  2). We say the gap associated with such an interaction is closed if r(S i1 , . . . , S ip ) < 0. The fraction of the interactions or contacts whose gaps are closed can be written as where g(r) is the distribution function of the gap defined in Eq. (16). This means there are N f closed constrains. Then isostaticity implies, or in M → ∞ limit.

III. SUPERCOOLED SPIN LIQUID STATES, CRYSTALLINE STATES AND THEIR STABILITY
In this section we focus on the crystallization and possibility of super-cooling, realization of supercooled paramagnetic state which is at least locally stable against crystallization. This is an important step toward realization of glasses without quenched randomness. In the present section we consider the disorder-free model Eq. (4). The effect of quenched randomness will be discussed in sec. V B.
A. Crystalline order parameter and the free-energy functional Our disorder-free models given by the Hamiltonian Eq. (2), Eq. (3) and Eq. (4) have the following global symmetries. In sec. II A we introduced two types of spins: Ising and continuous spins. In the cases of Ising spins S µ i = ±1, and for even p, the system has a global symmetry with respect to S µ → −S µ for each component µ. Such symmetry is absent for the cases of odd p. In the cases of continuous spins S µ i ∈ R, and for p = 2, the system has a global continuous symmetry with respect to rotations of spins in the M -dimensional spin space. The continuous symmetry is lost for p > 2 and the residual global symmetries become just the same as those in the Ising cases.
To be more specific, suppose that the system has a ferromagnetic ground state S i = (1, 1, . . . , 1) for ∀i. This is achieved for example by choosing the linear potential V (x) = Jx with J > 0 Eq. (11). Because of the global symmetries mentioned above, there can be other equivalent ground states, e. g. S i = (−1, −1, . . . , −1) for ∀i. In order to study the possibility of spontaneous symmetry breaking which select one ground state out of the equivalent ones, we may apply an external field of strength h > 0 parallel to the ground state (1, 1, . . . , 1), (21) and examine the behavior of an order parameter, where · · · h represents a thermal average in the presence of the symmetry breaking field. The standard procedure to analyze the problem is as follows. 1) One first construct a free-energy −βG(h) in the presence of the field h and then perform a Legendre transform to obtain −βF (m) = −βG(m) + N mh and then 2) seek for a solution which satisfy ∂ m (−βF (m)) = h = 0. In addition, since we are considering to take the limit of large number of components M → ∞, we may also define a local order parameter,

Spin trace
The above discussion motivates us to introduce an identity, The integration over h and m corresponds to the steps 1) and 2) mentioned above. Using the identity spin traces can be expressed formally in M → ∞ limit as, Here the integration over h can be done (formally) by the saddle point method in the limit M → ∞. The saddle point h * (m) is given by the saddle point equation, and we find, where h * = h * (m) is given by Eq. (26). More specifically, by taking the spin traces explicitly we obtain the following expressions for the Ising and continuous spin systems, • Continuous spin We find using Eq. (9), Here we and performed the integrations dS µ assuming λ > 0. Then we performed integrations over and λ and h by the saddle point method.
Using Eq. (25) we find, for example, In Eq. (25) we notice that different spin components µ are decoupled in the average µ . . . µ . Then we obtain the following cumulant expansion which will become useful in the following, Now the partition function Eq. (7) can be rewritten formally in M → ∞ limit as, where we defined where ∂ i represents the set of interactions which involve S i . Now we are left with the integrations over m i s in Eq. (33). We can take advantage of the fact that the lattice is a tree and that each spin is connected to c = αM neighbors which is a large number. Let us rewrite the above expression as, where we introduced a distribution function Here ∂ i \ i represents the set of vertexes attached to i excluding the vertex i itself. Then we find from Eq. (33), Eq. (35) and Eq. (36) that the integration over each m i can be done by the saddle point method in M → ∞ limit. It is natural to expect that the saddle point equations, admit a uniform solution m * i = m for ∀i. As the result we obtain the free-energy associated with such a saddle point as, with where m must satisfy the saddle point equation It is also required to satisfy the stability condition, B. Possibilities of the crystalline states We have just considered a ferromagnetic phase but we can also consider other crystalline states. For example, suppose that there is a crystalline ground state in which the spin configuration can be represented by some configuration (σ i ) 0 ∈ (−1, 1). Then it is useful to perform a gauge transformation S µ i →S µ u ≡ σ i S µ i . The crystalline order parameter m can be defined again as Eq. (22) but replacing the spins S by the gauge transformed onesS.
By the same gauge transformation the gap Eq. (3) (with X µ = 1) is transformed to, where we defined The variable η takes ±1 values. For simplicity we limit ourselves to the ground states such that it is a constant η = η for all the interactions . Then the results in the previous section Eq. (38)-Eq. (41) holds just by changing the argument of the potential as: The simplest example is p = 2 model with the linear potential V (x) = Jx but with J < 0. Obviously the ground state is the antiferromagnetic one : σ alternates the sign across each of the interactions (note that we are considering tree-like lattices with no loops). In this case η = σ 1( ) σ 2( ) = −1 so that it becomes essentially the same as a ferromagnetic model with J > 0 after the gauge transformation.
The stability condition becomes, • Continuous spin The saddple point equation Eq. (29) becomes, The stability condition becomes, It can be seen that the paramagnetic solution m = 0 always verify the saddle point equations. We are especially interested with the possibility that the paramagnetic state with m = 0 remains as a metastable state after the crystalline transitions take place so that glass transitions within the paramagnetic phase become possible.
• p = 2 case: -if |V (δ)| > 0, a 2nd order ferromagnetic transition takes place at a critical temperature below which the paramagnetic solution m = 0 becomes unstable and the ferromagnetic or antiferromagnetic order with |m| > 0 emerges continuously. If V (δ) is positive (negative) the ordering is ferromagnetic (antiferromagnetic) and we should choose η = 1 (η = −1). Since the paramagnetic state m = 0 is unstable below T c , supper-cooled paramagnet state is absent and thus glass transitions is not possible without quenched randomness.
-If V (δ) = 0, there will be no ferromagnetic nor antiferromagnetic phase transitions at finite temperatures.
The m = 0 solution remains stable at all finite temperatures, This is a very interesting situation where the crystallization is totally suppressed opening possibilities of glass transitions without quenched randomness.
• p > 2 case: The paramagnetic solution m = 0 remains locally stable at all temperatures in the sense of Eq. (50). Thus in this case supercooled paramagnetic state exist opening possibilities of glass transitions without quenched randomness.
1. Linear potential: p-spin ferromagnetic model As a simplest example let us consider the case of the linear potential V (x) = Jx where J > 0, which means V (δ) > 0. It is a ferromagnetic model so we choose η = 1. The saddle point equation becomes for the Ising spins, and for the continuous spins, We see m = 0 always verifies the saddple point equations as it should, • For p = 2 case, a 2nd order ferromagnetic transition takes place at a critical temperature k B T c /J = α √ M . Supper-cooled paramagnet and thus glass transitions without quenched randomness are not possible as discussed above.
• For p > 2, a 1st order ferromagnetic transition take place at k B T c /J = O(α √ M ). On the other hand the paramagnetic state m = 0 remains locally stable at all temperatures as discussed above. Quite interestingly in [20] a p = 3 Ising ferromagnet with M = 1 component was studied via cavity method and Monte Carlo simulations and the supercooled paramagnetic state and the glass transition were discovered. Our result is consistent with this observation.
Let us consider p → ∞ limit with γ fixed. Then we easily see that a 1st order ferromagnetic phase transition takes place at In the next section IV we will find that the system becomes equivalent to the random energy model (REM) [21] by excluding the ferromagnetic state.

Non-linear potentials with flatness
If the potential V (x) has a flat part where V (x) = 0, it tends to suppress crystallization and thus enhances the possibility to realize glass transitions inside the paramagnetic phase.
The simplest example may be the quadratic potential, Thus for p = 2 and δ = 0, the system should remain paramagnetic at all finite temperatures. More interesting case is the soft/hard core potential Eq. (13) which is completely flat for x > 0, Let us consider again p = 2 case. For δ > 0, V (δ) = 0 so that the system is paramagnetic at all finite temperatures. On the other hand for δ < 0, V (x) < 0 antiferromagnetic phase emerges via 2nd order transition. Then we choose η = −1. The transition temperature Eq. (49) is found as, In the hardcore limit → ∞, the antiferromagnetic transition takes place as δ → 0 + . In the following we we will find that the system becomes essentially identical to the random energy model (REM) [21] in the p → ∞ limit as far as the supercooled states are concerned. This proves the existence of the self-generated randomness.
The hamiltonian is given by with J > 0. Here we are especially interested with the p → ∞ limit with γ = α/p introduced in Eq. (53) fixed. As we discussed in sec III C it exhibits a ferromagnetic phase transition at k B T c = γJ √ M / ln 2. We examine the distribution of the energies of the disorder-free model performing a similar analysis done for the p-spin Ising spinglass model with quenched randomness in [21]. To this end let us first introduce a flat average over the 2 N M spin configurations, Then the distribution of energy among all configurations is obtained as, Here (1)), since we take the p → ∞ limit with fixed γ as defined in Eq. (53). Next let us examine simultaneous distribution of energy E 1 associated with an arbitrary chosen spin configuration (S 1 , S 2 , . . . , S N ) and the energy E 2 of another configuration (S 1 , S 2 , . . . , S N ). Here the latter is created from the former by flipping, say according to a deterministic rule, a fraction (1 − q)/2 with 0 < q < 1 of the elements of the former. In other words the overlap between the two configurations is Here we evaluated the expectation value . . . S by performing expansion Here we realize that in the p → ∞ limit, the distribution function decouples, because 0 < q < 1. The above observations imply that the present ferromagnetic system without any quenched randomness behaves essentially as a REM in p → ∞ limit: over the majority of the 2 N M spin configurations , excluding the negligible fraction of the spin configurations close to the ferromagnetic ground state, it is as if each microscopic states is assigned a random energy drawn from a Gaussian distribution with 0 mean and variance √ γJ. The we see that all the exact results of the standard version of the REM [21], which corresponds to γ = 1/2, can be used in the present system by replacing J of the standard model by √ 2γJ. Then we readily find that the system exhibit the Kauzmann transition, i. e. ideal static glass transition at, within the supercooled paramagnetic states. At temperatures below T K , the internal energy becomes stack at E/N M = − √ 2γ ln 2J among the disordered states, while the ferromagnetic ground state energy is given by In this section we setup a formalism to study the glass transitions using the replica method. We first develop a free-energy functional of the disorder-free model Eq. (4). Then we also consider the model with quenched randomness Eq. (5).
A. Disorder free model We consider a system of replicas a = 1, 2, . . . , n of the disorder-free model Eq. (4), where The free-energy of the replicated system can be expressed as, with the replicated partition function where Tr a i represents a trace over the spin space in replica a. In order to detect spontaneous glass transition, we can follow steps analogous to the one we took in sec. III A for the crystalline (ferromagnetic) transition. Namely we can explicitly break the replica symmetry as [35], and study the behavior of the glass order parameter matrixQ Here . . . represents the thermal average in the presence of the symmetry breaking field ab . Note also that Q aa = 1 due to the spin normalization. Just as the case of ferromagnetic transition discussed in sec. III A, we can consider the following steps to analyze the problem: 1) One first construct a free-energy −βG( ) in the presence of the field and then perform a Legendre transform to obtain −βF (Q) = −βG(Q) + N Q and then 2) seek for a solution which satisfy ∂ Q (−βF (Q)) = = 0.
Since we are considering M → ∞ limit we may also define a local glass order parameter, The above discussion motivate us to introduce an identity, (71) The integration over ab and Q ab corresponds to the steps 1) and 2) mentioned above. Using the latter spin traces can be expressed formally in M → ∞ limit as, Here the integration over ab can be done by the saddle point method in M → ∞. The saddle point equation which determines the saddle point * ab (Q) is given by, and we find, where * ab = * ab (Q) determined by Eq. (73). More precisely, by taking the spin traces we obtain the following expressions for the Ising and continuous spin systems, • Continuous spin: We find using Eq. (9) and introducing aa = λ a and Q aa = 1 (spin normalization) for convenience, Here we have performed integration over ab by the saddle point method which yields a saddle point * =Q −1 (76) • Ising spin: We find using Eq. (8), Here we performed the spin trance formally as Tr Sc e − a<b ab S a S b = Tr Sc e − a<b ab Using Eq. (72) we find, for example,

Evaluation of the free-energy
In Eq. (72) we notice again that different spin components µ are decoupled in the average µ . . . µ . Then we can evaluate the spin trace in the replicated partition function Eq. (67) in the M → ∞ limit using the cumulant expansion Eq. (31) and Eq. (79) as, The above expression is a crucial result because it reveals the self-generated randomness in our 'disorder-free' model. Collecting the above results, the partition function Eq. (81) can be rewritten formally in M → ∞ limit as, where we defined To derive the last expression we used Eq. (10) and performed integrations by parts (see Eq. (A36) for the same calculation.). Now repeating the same argument as in sec. III A (see Eq. (35)), we realize that the integrations over each (Q i ) ab can be performed in the M → ∞ limit by the saddle point method. We can assume that admit a uniform solutionQ * i =Q for ∀i. As the result we obtain the free-energy associated with such a saddle point as, where −F int represents the interaction part of the free-energy (free-entropy). Importantly Q ab must satisfy the saddple point equations, It is also required to satisfy the stability condition, i.e. the eigen values of the the Hessian matrix, in n → 0 limit, must be non negative. The exact free-energy functional Eq. (84) with Eq. (85) can be derived also using a density functional approach as we show in appendix A for the case of continuous spins. It is done closely following the strategy used in the recent replicated liquid theory for hardsphere glass in large dimensional limit [36][37][38].

B. Interpolation between disorder-free and completely disordered model
In the present paper we are most concerned with systems where the randomness is self-generated. However it is very instructive to consider also the model with quenched randomness Eq. (5), Here ξ µ is a random variable with Gaussian distribution with zero mean and unit variance. With this parametrization, we have a continuous interpolation between the disorder-free model λ/ √ M = 1 and completely disordered, spinglass model λ/ √ M = 0. Analysis of this model is useful for the following reasons.
• We can show that the self-generated randomness and the quenched randomness act additively.  64). Actually the fact tat T K is lower than T c is natural by itself but they become too much separated in the disorder-free model in the M → ∞ limit. With the choice of the disordered model Eq. (5) we can bring the energy scales of two transitions much closer. This is achieved in two steps: (i) choose a much weaker interaction energy scale of order O(1/ √ M ) for the original 'disorder-free' model (ii) then add additional quenched randomness such that the energy scale of the glass transition is kept constant with respect to the variation of λ. Bringing the crystallization and glass transition temperatures closer, it becomes easier to discuss competitions between liquid, glass and crystalline phases.
Averaging over the quenched randomness of the replicated partition function Eq. (67) we obtain, where the overline denotes the average over the quenched-randomness and we introduced, For the order parameters we may consider both the glass order parameters Eq. (69) and the crystalline one Eq. (22). Since we are considering M → ∞ limit we naturally define the following set of local order parameters, and the corresponding identities, Using the identities shown above spin traces can be expressed in M → ∞ limit as, Here we introduced a short hand notationm = (m 1 , m 2 , . . . , m n ). The last equations are the saddle point equations for the integrations over ab and h a which fix the saddle points * ab = * ab [Q,m] and h * a = h * a [Q,m]. More precisely, by taking the spin traces we obtain the following expressions for the Ising and continuous spin systems, • Continuous spin: We find similarly to Eq. (75), and Here we have performed integration over ab and h a by the saddle point method which yield a saddle point • Ising spin: We find similarly to Eq. (74), where we performed the spin trace formally as Tr Sc e − a<b ab S a S b − a haS a = e − a<b ab Using Eq. (95) we find, for example,

Evaluation of the free-energy
In Eq. (95) we notice again that different spin components µ are decoupled in the average µ . . . µ . Then we can evaluate the spin trace in the replicated partition function Eq. (67) in the M → ∞ limit using the cumulant expansion Eq. (31) and Eq. (102) as, Here we point out that the last term in the exponent of the last equation is the summation of the contributions of two different different kinds of disorder: (1) quenched randomness of amplitude 1 − (λ/ √ M ) 2 (see the 2nd term in Eq. (90)) (2) self-generated randomnesss of amplitude (λ/ √ M ) 2 . Now it is now clear that parametrization Eq. (5) is chosen such that the energy scale of the glass transition does not change between the disorder-free limit (λ/ √ M = 1) and completely disordered limit λ/ √ M = 0. Collecting the above results, the disorder averaged replicated partition function Eq. (89) can be rewritten formally in M → ∞ limit as, As before in the M → ∞ limit we can assume that, admit uniform solutions :Q * i =Q andm * i =m for ∀i. As the result we obtain the free-energy associated with such a saddle point as, where Q ab and m a must satisfy the saddple point equations It is also required to satisfy the stability condition, i.e. the eigen values of the the Hessian matrix, in n → 0 limit, must be non negative. Finally let us note again that the disorder-free model can be recovered by choosing λ/ √ M = 1 in the above expressions. For the disorder-free model discussed in previous sections, we gave free-energy functional in terms of the crystalline order parameter m Eq. (39) and and that in terms of the glass order parameter Q ab Eq. (85) separately just to simplify the presentations. In any case here we now have complete free-energy functional where both the crystalline and glass order parameters are present.

C. Stability against crystallization
Given the complete free-energy functional in term of both the crystalline and glass order parameters, it is very interesting for us to investigate the stability of glassy phases against crystallization. Here assume limit ourselves with a glassy phase with vanishing crystalline order parameter m = 0 and do not consider possible 'glassy crystals' with m > 0. First we note that, holds. This can be checker by taking the derivatives explicitly. For the entropic part we find, The last equation follows from Eq. (101) which implies that m a = 0 requires h * a = 0. For the interaction part of the free-energy we find, The last equation holds for p > 1. Thus Eq. (111) must hold. Then the local stability of the glassy phase with m = 0 against crystallization is solely determined by the matrix, where we used ab .) Here we see that the 2nd term on the r.h.s of Eq. (114) vanishes in two cases (i) p > 2 (ii) p = 2 with non-linear potential with the flatness V (δ) = 0. Remarkably in these cases the matrix becomes independent of λ and its the eigen values are simply the inverse of the eigen values of the matrix Q ab . We expect the latters are positive for physical solutions (for instance note that in the case of the continuous spins the entropic part of the free-energy has a term ln det(Q) with m a = 0 (see Eq. (97)). Thus the eigenvalues of the matrix Q ab are needed to be positive.) The above observation is very interesting because including the regime of large enough λ where we naturally expect crystalline order as the true equilibrium phase, supercooled paramagnetic phase m = 0 (for which Q ab = δ ab ) and also the glassy phase with m = 0 remains locally stable against crystallization for the two cases: (i) p > 2 (ii) p = 2 with non-linear potential with the flatness V (δ) = 0. Note that this is a generalization of our analysis in sec. III C which was limited to the supercooled liquid (paramagnet) regime (before the glass transition).
Finally we note that the situation can change in systems with finite connectivity. The supercooled paramagnetic phase can disappear for sufficiently large λ. In the context of statistical inference problems this is an important issue because one has to find the hidden crystalline state (ground truth) in the immense sea of wrong solutions (glasses) [39].

VI. LINEAR POTENTIAL: CONNECTION TO THE STANDARD p-SPIN ISING/SPHERICAL SPINGLASS MODELS AND THE RANDOM ENERGY MODEL
Let us discuss here the simplest case, the linear potential Eq. (11) which reads, The interaction part of the free-energy Eq. (107) becomes, where s ent (Q,m) is given in Eq. (97) for the continuous spin case and Eq. (100) for the Ising case. This is exactly the same as those of the standard p-spin Ising/spherical spinglass models with M = 1 but with global couplings [21] [40,41] by choosing α/p = 1/2. Note that such a correspondence has been known for the case of a p = 2 continuous spin model with global coupling [31]. Let us summarize below some important known results of the p-spin spinglass models. The case p = 2 with the Ising spin corresponds to the SK (Sherrington-Kirkpatrick) model [42]. It exhibits a continuous phase transition from the paramagnetic to the spinglass phase accompanying the full replica symmetry breaking (RSB) [43] while the spherical version of it exhibits a continuous phase transition but without RSB [44]. The SK model is the standard mean-field model for spinglasses. On the other hand p > 2 system exhibit 1 step RSB [40] with a discontinuous transition from the paramagnetic to spinglass phase. These models show the essence of the glass phenomenology such as the dynamical and static glass transitions so that they are regarded as prototypical theoretical model to capture the physics of structural glasses [6]. Among the latter models with the Ising spin exhibit yet another glass transition to enter the full RSB phase at lower temperatures [45]. In p → ∞ limit the random energy model is recovered [21].
We emphasize that the result Eq. (118) is valid also in the disorder-free limit λ/ √ M = 1. Indeed in sec. IV we have shown that the random energy model is recovered in the disorder-free limit. Thus the disorder-free model have sufficient amount of self-generated randomness to realize glass transitions. The supercooled paramagnetic state and the glass phase which emerge there are stable against crystallization for p > 2 as discussed in sec. V C. In the case p = 2, however, we have to use the model with quenched randomness to realize glassy phases. (This amount to yield nothing but the SK model for the Ising spins and spherical SK model for the continuous spins mentioned above.)

VII. REPLICA SYMMETRIC (RS) ANSATZ
For the rest of the present paper we study glass transitions of our model, which emerge within the supercooled paramagnetic phase with no crystalline order m = 0. And we limit our selves with the continuous spin models for the rest of the present paper. In the present section and in the next section we derive some generic results within the replica symmetric (RS) and replica symmetry breaking (RSB) ansatz. We apply these schemes to systems with non-linear potentials in later sections.
Our starting point is the free-energy functional Eq. (84) which reads, with Eq. (85) which reads, Here F int represents the interaction part of the free-energy. For the entropic part in Eq. (85) we used the expression given by Eq. (75) and we omitted irrelevant constants n 2 + n 2 ln(2π) for simplicity. The pressure Eq. (14) can be computed as, and similarly the distribution function of the gap Eq. (16) as, Before passing let us recall the discussion in sec III C 2 and sec. V C that the supercooled paramagnetic m = 0 states and glassy states of the model is locally stable against crystallization if p > 2. But for the case p = 2 we must have the flatness V (δ) = 0. Although we may not mention these points often in the following but we must keep these in our minds.

A. Formulation
In the replica symmetric (RS) ansatz we assume the following form of the overlap matrix parametrized by a single parameter q, so that the entropic part of the free-energy is obtained as The interaction part of the free-energy is obtained as, where we used the formula exp a 2 and the following short hand notations: γ a (x) is a Gaussian with zero mean and variance a, by which we write a convolution of a function A(x) with the Gaussian as, where Collecting the above results we obtain the variational free-energy Eq. (119) within the RS ansatz as − βf RS (q) = ∂ n s RS (q)| n=0 = 1 2

The saddle point equation
The saddle point equation for the order parameter q is obtained as, where we introduced

Pressure and distribution of gap
Using Eq. (131) we obtain the pressure Eq. (14) as, and similarly the distribution of the gap Eq. (122) as, which is properly normalized such that drg(r) = 1. One can also check easily that the 'virial equation' Eq. (17) for the pressure is satisfied, as it should be.

p = 2 case
From the results in appendix B 1 we find for the q = 0 solution for p = 2, from which we find the eigenvalues of the Hessian matrix as (see Eq. (B26)), which vanishes at, For α < α c (δ), the eigenvalues are positive so that the q = 0 solution is stable but it becomes unstable for α > α c (δ).
Interestingly we see that at the critical point α = α c (δ), q = 0 solves also G(q) = 0 (see Eq. (133)). Since q = 0 solution must solve G(q) = 0, this suggests a possibility of continuous phase transition at the critical point such that q = 0 solution emerges continuously. The situation is similar to the Sherrington-Kirkpatrick model which exhibits a continuous spinglass transition at the d'Almeida-Thouless (AT) instability [46] line.

p > 2 case
For p > 2, using the results reported in appendix B 1, we find λ R = λ L = λ A = 2 > 0 (see Eq. (B26)) so that q = 0 solution is always stable. Thus contrary to the p = 2 model, we find the liquid phase described by the q = 0 RS solution is always (meta)stable. The situation is very similar to the usual p-spin spherical spinglass models [40]. Then we are naturally lead to consider the possibility of a discontinuous glass phase represented by 1 step replica symmetry breaking (1RSB) much as the usual p-spin SG models (see sec VI) including the random energy model (see sec. IV)). The latter is the standard random first order (RFOT) scenario [5,6] which is established theoretically for hardspheres in d → ∞ limit [9].
Here we continue the previous section and obtain some generic results based on the replica symmetry breaking (RSB) ansatz for the M -component continuous vector spin system with generic potential V (x).

A. Parisi's ansatz
We assume the following structure of the glass order parameter in the glass phase, which is the Parisi' ansatz [43], where I m ab is a kind of generalized ('fat') identity matrix of size n × n composed of blocks of size m × m. (see Fig. 3) The matrix elements in the diagonal blocks, are all 1 while those in the off-diagonal blocks are all 0. The Parisi's matrix has a hierarchical structure such that which becomes in the n → 0 limit. The expression Eq. (140) becomes valid also for the diagonal part by introducing which reflects the normalization of the spins. Let us note that we may sometimes extend the labels m's Eq. (142) introducing an additional label m k+2 just for conveniences when we deal with some recursion formulas.

B. Free-energy
Let us evaluate the free-energy Eq. (119) using the above ansatz. To compute the entropic part of the free-energy one needs to evaluate ln detQ k−RSB . Given the hierarchical structure of the Parisi's matrix, one can obtain them in a recursive fashion and one finds [47], Remembering that m 0 = n we find, with which implies The interaction part of free-energy can also be evaluated in a recursive fashion [48]. One finds, In the 2nd equation of Eq. (149) we used I m k+1 =1 ab = δ ab and introduced where we used the formula Eq. (127). The expression Eq. (149) naturally motivates a family of functions g's which obey a recursion relation, for l = 0, 1, . . . , k. Then it is easy to see that Equivalently we can introduce a related family of functions f 's, which follows a recursion relation, for i = 0, 1, . . . , k + 1 with the boundary condition, where m k+1 = 1. We may also express the boundary condition as, by introducing m k+2 just as an additional label for convenience. Remembering that m 0 = n we find the interaction part of the free-energy becomes, Finally collecting the above results we obtain the free-energy within the k-RSB ansatz as, with where we used G 0 = 1 − (m 1 − m 0 )q 0 and that m 0 = 0 and m 1 = 1. In the 2nd equation we used Eq. (158). The result agrees with Eq. (131) as it should.

k = 1 case: 1 RSB
For the k = 1 RSB case we find, where q 0 and q 1 must be determined through the saddle point equations which we discuss in sec. VIII E 2.
An important quantity is the complexity or the configurations entropy Σ(f ), which describes the exponentially large number of states ∝ e N Ξ(f ) df with free-energy density between f and f + df in the glass phase. Using Monasson's prescription [49], one can construct the complexity function Ξ(f ) using m = m 1 as a parameter; Thus extremization of the free-energy with respect to m, 0 = ∂ m βf 1−RSB (q 0 , q 1 , m) amounts to force the complexity to vanish Σ * (m) [49]. We readily find the following explicit expressions, ln Dz 2 e −βV (Ξ(δ,q0,q1)) and ln Dz 2 e −βV (Ξ(δ,q0,q1)) 3. k = ∞ case: continuous RSB In the limit k → ∞, the overlap matrixQ is parametrized by a function q(x) with n < x < 1. Then Eq. (145) becomes, From the above expression we find with Then the free-energy Eq. (161) can be written as, where the function f (x, h) obeys,ḟ Here and in the following we denote a partial derivative with respect to the 1st argument by a dot, e. g. ∂ x f (x, h) = f (x, h) and that with respect to the 2nd argument by a dash e. g. ∂ h f (x, h) = f (x, h). The partial differential equation Eq. (174) is the continuous limit of recursion formula Eq. (157). The boundary condition Eq. (158) becomes, C. Variation of the interaction part of the free-energy Later we will often meet the needs to consider variation of the free-energy. This happens when we wish to compute the pressure Π Eq. (121), the distribution of the gap g(r) Eq. (122) and also to obtain the saddle point equation for the order parameters q l 's (sec. VIII E). Here we consider a strategy to deal with the variation of the interaction part of the free-energy Eq. (149).
As we discussed in sec. VIII B, the interaction part of the free-energy is is constructed in a recursive way. This fact naturally motivates us to introduce for 0 ≤ i ≤ j ≤ k + 1, Using the chain rule we can write, where as one can easily find from the recursion relation Eq. (157). Then we find the recursion relation, with the 'boundary condition' Here ⊗ h stands for the convolution with respect to the variable h. A useful property to note is that the recursion relation Eq. (180) preserves the 'normalization', which can be easily proved using Eq. (157). Now let us discuss how to analyze variations of the interaction part of the free-energy Eq. (149). The interaction part of thee free-energy is given by Eq. (160) which reads as, . It is convenient to introduce a short handed notation, where we used the chain rule, Eq. (179) and set m 0 = n → 0. Clearly it follows the same recursion formula as Eq. (180), with the 'boundary condition' which follows from Eq. (183) and Eq. (181). The functions P (m i , h) is also normalized such that reflecting Eq. (182).
As an application, let us show that P (m k+1 , r) is nothing but the distribution function of the gap g(r) defined in Eq. (16). The distribution of the gap g(r) Eq. (122) can be obtained using Eq. (161), Eq. (159) which reads βV (r) = f (m k+2 , r), Eq. (177), Eq. (183) and Eq. (184) as, We used m k+1 = 1 in the last equation. It can be seen that g(r) is properly normalized such that drg(r) = 1 reflecting Eq. In the k → ∞ limit, the function P (x, h) can be obtained by solving a differential equation, which is the continuous limit of Eq. (184). The boundary condition Eq. (185) becomes,

D. Pressure
The pressure Eq. (121) for the generic k-RSB ansatz is obtained using and Eq. (161) as, where we introduced, Here and in the following the prime stands for taking a partial derivative with respect to the 2nd argument h. The function π(m, h) also follows a recursion formula which can be obtained using Eq. (157) and Eq. (158) as, for i = 1, 2, . . . , k with the 'boundary condition', The previous result in the RS case (k = 0) Eq. (134) can be recovered using Eq. (193) in Eq. (190). It is also easy to verify that the pressure Eq. (190) can be recovered through the virial equation for the pressure Eq. (17). In fact the pressure can be expressed as, with any i = 1, 2, . . . , k + 2 ( so that π = dhP (x, h)π(x, h) for any 0 < x < 1 in the k → ∞ limit). Using the recursion formulas Eq. (184) and Eq. (192) one can check that dhP (m i−1 , h)π(m i , h) = dh P (m i , h )π(m i+1 , h ) so that the r. h. s of the above equation does not depend on the level i of the hierarchy. The case of the virial equation for the pressure Π = drg(r)(−βV (r)) Eq. (17) corresponds to the case i = k + 2 which can be seen by noting π(m k+2 , h) = −f (m k+2 , h) = −βV (h) (see Eq. (158)) and Eq. (187). On the other hand, the case i = 1 corresponds to the expression Eq. (190) which can be seen using Eq. (185). In the k → ∞ limit, the function π(x, h) ≡ −f (x, h) obeys a differential equation, which is the continuous limit of Eq. (192) and can be obtained from Eq. (174). The boundary condition for the latter is given by Eq. (193) which reads,

E. Saddle point equations for the order parameters
Here we derive variational equations to determine q i for i = 0, 1, 2, . . . , k. Since q i 's are related linearly to G i 's through Eq. (148), the saddle point equations can be written as, In the last equation we used Eq. (151) and Eq. (150). As we show in appendix C we find, where P i,j (y, h) is defined in Eq. (177) and π(m, h) is defined in Eq. (191) Collecting the above results we obtain the variational equations as where we introduced where P (m j , h) and π(m i , h) can be obtained by solving the recursion formulas Eq. (184) and Eq. (192) respectively together with their boundary conditions. We note that for p > 1, q 0 = 0 always solves the 1st equation of Eq. (199).  For the k = 1 case (1RSB) we have, After solving the above equations for G 0 and G 1 , the order parameters q 0 and q 1 can be obtained as (See Eq. (148)), To evaluate κ 0 and κ 1 in Eq. (203) we need more information. Suppose that we are given some initial guess for the values of q 0 and q 1 . Then we can recursively obtain functions f (m 2 , h) and f (m 1 , h) (see Eq. (158)) and Eq. (157)) as, where m 2 = 1. Similarly we can recursively obtain functions π(m 1 , h) and π(m 2 , h) (See Eq. (192) and Eq. (193)) as, Next we can recursively obtain functions P (m 0 , h) and P (m 1 , h) (see Eq. (184) and Eq. (185)) as, With these we are now readily to evaluate κ 0 and κ 1 using Eq. (203).
We note that the parameter m 1 remains. In order to study the equilibrium state m 1 is fixed by the condition of vanishing complexity Σ(m 1 ) = 0. (See sec. VIII B 2) The saddle point equations for a generic finite k-RSB ansatz with some fixed values of 0 < m 1 < m 2 < . . . < m k < 1 can be solved numerically generalizing the procedure explained above. 0. Make some quesss for the initial values of q i (i = 0, 1, 2, . . . , k). Then compute G i for i = 0, 1, . . . , k using Eq. (147). 4. Return to 1.
The above procedure 1-4 must be repeated until the solution converges.

k = ∞ case: continuous RSB
In the limit k → ∞, the variational equations Eq. (199) become, From the above equations we can derive some exact identities which become useful later. Taking a derivative with respect to x on both sides of Eq. (209) and using Eq. (210), Eq. (172), Eq. (188), Eq. (195), we find after some integrations by parts, Then taking another derivative on both sides of the above equation we find after some integrations by parts, F. Stability of the kRSB solution Stability of the k(≥ 1)-RSB ansatz must be examined by studying the eigenvalues of the Hessian matrix. As we note in sec B 2, there is a residual replica symmetry within each of the inner-core part of the replica groups. Here we do not study the complete spectrum of the eigen-modes of the Hessian matrix but focus on the replicon eigenvalue responsible for the replica symmetry breaking of the residual replica symmetry.

IX. QUADRATIC POTENTIAL
Now we are ready to study specific problems with non-linear potential V (x). Here we briefly study the simplest one, As we see below it is already a non-trivial problem. First we note that the interaction part of the free-energy Eq. (120) can be expanded in power series of the order parameter as, In contrast to the case of the linear potential discussed in sec VI, higher order terms of Q p ab appears. Thus even with p = 2 for which system remains RS for the linear potential case (spherical SK model [44]), it is natural to expect that the quadratic potential allow RSB because it is appears similar to the 2 + p spherical model [50] which exhibits various types of RSB.
In the present paper we do not explore the whole phase diagram but let us show the existence of full RSB for the p = 2 case. From Eq. (139) we find, which implies We solved numerically the continuous RSB equation approximated by k-step RSB (With k = 200) recursion formulas with appropriate boundary conditions as explained in sec. VIII E 3. In Fig. 4 we show the continuous RSB solution obtained numerically for the case of p = 2, α = 2 and δ = 1.0 at various temperatures below T c / = √ 2 − 1 = 0.414... As expected the glass transition takes place continuously. Moreover it accompanies full RSB.
Finally let us recall the discussion in sec. V C that for the case p = 2, we need the flatness V (δ) = 0 in order to ensure the locally stability against crystallization. However the above results imply α c (δ = 0) = ∞. Thus for the p = 2 case with the quadratic potential, we cannot avoid using the disordered model Eq. (5) in order to realize the glass transition.

X. HARDCORE POTENTIAL
We will now focus on a more strongly non-linear potential, namely soft/hardcore potential, which becomes a hardcore potential in the limit → ∞. Note that in the special limit p = 1 and with fully disordered choice λ/ √ M = 0 in Eq. (5), the model becomes identical to the perceptron problem [32]. Recent works [28,29] have established that the universality class of the jamming in the perceptron model with δ < 0 is the same as that of hardspheres [11]. We wish to clarify if the same universality holds for p ≥ 2 or not.
As we discussed in sec III C 2 and sec. V C, the supercooled paramagnetic states and glassy states of the model is locally stable against crystallization if p > 2. For the case p = 2 we must have V (δ) = 0 which means δ > 0 for the soft/hardcore potential (see Eq. (56)). In the following we study both δ > 0 and δ < 0 but we should keep this point in our mind.
Below we closely follow the analysis done on hardspheres in d → ∞ limit [11] and find indeed that many aspects are quite similar to those found there, especially at jamming.

A. Replica symmetric solution
Let us first study the RS solution discussed in sec. VII in the case of the hardcore model.

Free-energy
For the hardcore potential we find the RS free-energy Eq. (131) as, where we introduced a function Θ(x) with erf(x) being the error function, which behaves for x → ∞ as, This implies, The function G(q) defined in Eq. (133) becomes, where we introduced, which behaves asymptotically as, x → ∞ (227)

q = 0 RS solution and its stability
Within the liquid state q = 0, the p-dependence disappears. The pressure is obtained as, As shown in the left panel of Fig. 5, the pressure monotonically increases by decreasing δ as expected. We display in the right panel of Fig. 5 b) the behavior of the distribution of gap, We see that the peak around r = 0 develops by decreasing δ as expected.
As we found in sec. VII B q = 0 solution is always stable for p > 2 body interactions. For the p = 2 body hardcore interaction, we find the q = 0 solution becomes unstable for α > α c (δ) with which is obtained from Eq. (139). The result is shown in Fig. 6. . Note that these are independent of p.

Jamming within the RS ansatz
It is possible to look for glass transition within the RS ansatz by looking for q = 0 solution of the RS saddle point equation Eq. (132), which must solve G(q) = 0. In sec X D 1 we will examine the phase diagram within the RS ansatz for p = 2 case. Here we instead focus on the jamming limit where the order parameter saturates q → 1 within the RS ansatz.
The location of the jamming point can be analyzed as follows. We find G(q) given in Eq. (225) becomes in the limit q → 1, In the last equation we used the asymptotic behavior of the function r(x) given in Eq. (227). Thus we find the jamming line α = α j (δ), which is also displayed in Fig. 6. The pressure Eq. (134) becomes for the hardcore model, where we used the asymptotic expansion Eq. (227). Thus as expected the pressure diverges by jamming (see Fig. 7 b)). Next let us examine the distribution of the gap Eq. (135) which becomes for the hardcore model, The behavior in the jamming limit q → 1 can be viewed in the following two ways (see Fig. 7 c)).
2. In the vanishing region around r = 0 parametrized as r = (1 − q p )λ we find a different behavior as follows.
Assuming q ∼ 1 we find for r > 0, In the 1st equation we dropped contribution from δ −∞ dz 0 .. which can be neglected compared with the contribution from ∞ δ dz 0 .. and used the asymptotic behavior of the error function Eq. (223) which implies Θ(−X) ∼ e −X 2 2 √ πX for X ≫ 1. Thus in the jamming limit q → 1 we find diverging peak in the "contact region" around r = 0 whose height diverging as 1/(1 − q) and the width vanishing as 1 − q.

B. k-step RSB solution
We study now the k-RSB solution discussed in sec. VIII in the case of the hardcore model.

Inputs
Here we present some necessary inputs to study the glass phase and jamming of the hardcore model within generic k-RSB ansatz. With the following inputs, 1RSB solution can be obtained following sec. VIII E 2 and generic k-RSB solution can be obtained following sec. VIII E 3.
For the hardcore potential Eq. (219) we find where Θ(x) is defined in Eq. (221). Then The functions f (m i , h) and π(m i , h) are determined via recursion formulas Eq. (157) and Eq. (192) using the boundary values obtained above.
It is useful to study the asymptotic behavior of the functions f (m i , h) and π(m i , h) in the limit h → −∞ both for numerical and analytical purposes. Using Eq. (221) and Eq. (223) and the recursion formula Eq. (157) one finds for i = 1, 2, . . . , k + 1, Note thatΛ k+1 = m k+1 Λ k+1 = Λ k+1 = 1 − q p k . In the continuous limit k → ∞ this implies, with λ(x) = q p (x) defined in Eq. (175). The above observation suggests us to introduce a function j(m i , h) defined as, From Eq. (157) we find that the function j(m i , h) follows a recursion relation, with and the boundary condition, Correspondingly one finds that Eq. (192) becomes,

Rescaled quantities useful close to jamming
In the hardcore model jamming means disappearance of the thermal fluctuation q 1 → 1. Let us show below how to modify the numerical algorithm in sec VIII E 3 to solve the full RSB equations close to jamming. To this end let us first introduce several rescaled quantities.
The above steps 1-6 must be repeated until the solution converges.

Algorithm to look for the jamming point
We can also look for the k-RSB solution for a given, fixed ∆. This can be seen as the following. In the step 5 of the procedure explained above we obtain γ k using Eq. (261) but γ k = 1. Thus Eq. (261) for i = k can be considered as an equation to determine α = α(∆). In particular, the jamming point α J (δ) via k-RSB ansatz can be determined by choosing ∆ = 0.

C. Jamming criticality
Let us discuss properties of the system approaching the jamming. We expect the q(x) function of the full RSB solution has a continuous part for some range x m < x < x 1 and a plateau q(x) = q 1 for x 1 < x < 1. In the hardcore model jamming means disappearance of the thermal fluctuation q 1 → 1 as mentioned before. For a convenience we define, Then jamming implies ∆ 1 = 1 − q 1 → 0. We discuss below properties of the system encoded in the full RSB solution in the vicinity of the core x → x 1 which encodes physical properties of the system in the deepest part of the energy landscape. In the following we will find results very similar to those found in the hardspheres in d → ∞ [11] where it was shown that continuous RSB solution gives a qualitatively different result from finite k-RSB ansatz concerning the scaling behavior approaching jamming.
Following [11] and [28] we consider the following scaling ansatz at the core x → x 1 in the jamming limit ∆ 1 → 0, with an exponent κ.
From Eq. (188) and Eq. (195) we have, Based on the asymptotic behavior of the function π(x, h) given in Eq. (240) we expect, For x → x 1 and ∆ 1 → 0 we can assume, This implies the following scaling form for x → x 1 , with some scaling function P 0 (x).
2. Matching between (1) and (2): assuming we find the following relation is needed to eliminate the dependence on ∆, (266) we find the following two ordinary differential equations, which are subjected to the boundary condition One can check that the differential equations Eq. (285) with the boundary condition Eq. (286) is consistent with the scaling relations for θ and α given by Eq. (278) and Eq. (282).
Here we notice that the apparent dependence on p in Eq. (285) can be formally eliminated by the following replacement, This means that if we find a solution for the p = 1 case, the solutions for other values of p can be obtained using Eq. (287) in the reversed manner. Importantly such a solution satisfies the same desired asymptotic behaviors Eq. (286). This implies the universality does not change with p.
However as pointed out in [11] the above equations do not completely solve the problem. We are left with the exponent a undetermined while other quantities P 1 (z), π 1 (z) and the exponent c can be obtained in a form parametrized by a. (All other exponents are fixed given a and c.) The final task to fix the value of the exponent a can be done using the exact identity Eq. (212) which reads in the limit x → x 1 and q 1 = q(x 1 ) → 1 where we defined, We notice that the contribution of dhT 1 (h) vanishes for the p = 1 case. Thus we must carefully examine whether dhT 1 (h) remain relevant in the jamming limit ∆ 1 = ∆(x 1 ) → 0 or not. We examine contributions of the integrals dhT 1 (h) and dhT 2 (h) from the regime (1) h → −∞ and (2) where we took leading terms for the jamming limit ∆ 1 → 0. Similarly using Eq. (273),Eq. (283) and Eq. (269) we find for the regime (1) h ∼ 0, Collecting the above results we find the relevant contribution in the jamming limit ∆ 1 → 0 is given by regime(1) dhT 2 (h). It means that we must satisfy, ∞ −∞ dzP 1 (z) (π (z)) 2 − 2 c κ 1 p (π 1 (z)) 3 + (π (z)) 2 = 0 (293) Again we find the apparent p dependence can be formally eliminated by the replacement Eq. (287). Based on the above analysis we can conclude that the critical exponents and the scaling functions P 1 (z), and π 1 (z) does not dependent on p, i. .e. super-universal.

Divergence of the pressure
The pressure can be expressed as Eq. (194) which reads, where x can be chosen arbitrary. Using the scaling ansatz Eq. (271) and Eq. (269) at the core x → x 1 and jamming ∆ 1 → 0 we find contribution from largely negative region of h becomes Similarly we can analyze contribution from the region h ∼ 0 using Eq. (273), and Eq. (269) If a < 1, which is the case, the latter gives a only sub-dominant contribution. To sum up we find, the 'cage size' ∆ 1 vanishing in the jamming limit Π → ∞ as,

Distribution of gap
For the hardcore model the distribution of the gap g(r) within the k-RSB ansatz Eq. (187) reads, 1. For fixed finite r, sending q k → 1, we find, where we used lim X→∞ Θ(X) = 1. This is a generalization of the RS (k = 0) result Eq. (235). In the k → ∞ limit, the scaling behavior of P (x, h) close to the core x → x 1 as described by Eq. (271) and Eq. (273) in the region vanishing in the jamming limit ∆ 1 → 0 implies development of a delta peak δ(r). On the other hand, we have the scaling behavior P (x, h) ∼ h −α for fixed h ∼ 0 + as given by Eq. (281) with α = a/b Eq. (282). These observations implies, where c nt is some numerical factor.
2. In the vanishing region around r = 0 parametrized as r = (1 − q p k )λ we find, Assuming q k ∼ 1 we find for r > 0, This is a generalization of the RS (k = 0) result Eq. (235). Now in the k → ∞ limit we have the scaling behavior P (x, h) ∼ ∆ −c/κ P 0 (h∆ −c/κ ) for h < 0 Eq. (271). Using this for x → x 1 we find, where we used c = κ − 1 and 1 − q p k p∆ 1 for ∆ 1 → 0. Using the above result we can evaluate the fraction of interactions or contacts which is closed. For any small but finite we have, Thus in the jamming limit, the faction of closed contact Eq. (18) ca be expressed as,

Isostaticity
Let us consider whether isostaticity discussed in sec. II C holds in the jamming limit in the present model. The condition of isostaticity Eq. (20) becomes in M → ∞ limit with α = c/M fixed at jamming ∆ 1 → 0 becomes, where we used Eq. (304). Actually using the exact identity Eq. (211) which holds for the continuous RSB together with the scaling behavior Eq. (271) in the h < 0 region and the relationΛ(x) pG(x) given by the 2nd equation of Eq. (269) which hold close to the core x → x 1 at jamming ∆ 1 → 0 we find, Thus we see that the isostaticity holds at jamming. Note that the term which is proportional to p − 1 apparently violates the isostaticity but it scales as ∆ 2c/κ 1 and becomes irrelevant in the jamming limit ∆ 1 → 0 as long as c/κ > 0.
D. Detailed analysis on the p = 2 case Let us take here the p = 2 case and study the model more in detail to work out the phase diagram and behavior of physical quantities. We have found that the system exhibit continuous transition to antiferromagnetic phase for δ < 0 at T c given by Eq. (57). Thus the following analysis is valid for the disorder-model in the range δ > 0. For the other side δ < 0, we have to eliminate the antiferromagnetic phase using the completely disordered model with λ/ √ M = 0 in order to realize the glassy phases.

RS solution
For the p = 2 case, we find from Eq. (230) that the paramagnetic solution q = 0 becomes unstable at the critical point α c (δ). Then we are naturally led to examine the possibility of the q = 0 solution. Within the RS ansatz, it must solve G(q) = 0 (see Eq. (132)), where the function G(q) for the hardcore model is given by Eq. (225). Expanding G(q) up to order O(q 2 ) we find, where The above equation can be solved for q to find, Thus we find that q = 0 solution emerges at the critical point α = α c (δ), where the q = 0 solution becomes unstable, and the EA order parameter q grows continuously increasing α.
In Fig. 6 we show the phase diagram for the p = 2 hardcore model within the RS ansatz. The glass transition line α = α c (δ) is given by Eq. (230). The jamming line α = α j (δ) is given by Eq. (232). δ c ∼ 0.51 (determined by α c (δ c ) = 4, see Fig. 6) and increases by decreasing δ and saturates to q = 1 approaching the jamming point. There we also show the behavior of the pressure Eq. (233) which diverges approaching the jamming and evolution of g(r) Eq. (234) which develops a diverging contact peak at r = 0 approaching the jamming. Finally let us examine the stability of this solution. From the result reported in appendix B 1 c we find the replicon eigenvalue λ R as, In the last equation we made an expansion in series of q which can be obtained by using r (x) = −2xr(x) − r 2 (x) which follows from Eq. (226). Comparing the function G(q) in Eq. (307) and R(q) in Eq. (312) we notice that they are identical up to O(q) but different in the O(q 2 ) terms. Using Eq. (309) in Eq. (312) we find up to O( 2 ), It turns out that A(x) is a monotonically decreasing function of x. Using the asymptotic behavior of the function r(x) given in Eq. (227) one can find lim x→−∞ A(x) = −1/4. Thus we find that the replicon eigenvalue is definitely negative meaning that the RS solution is unstable for α > α c (δ). We also checked numerically, solving G(q) = 0 for q and evaluating λ R (Eq. (311)) that this is indeed the case in the whole regime of α > α c (δ). Thus the replica symmetry must be broken for α > α c (δ). Remarkably the situation is very similar to the Sherrington-Kirkpatrick (SK) model for the spin glasses [42,46,51]. To summarize we find the liquid solution described by the q = 0 RS solution which becomes unstable approaching the critical point α c (δ) where all eigenvalues of the Hessian matrix vanish. The line α = α c (δ) is the equivalent of the d'Almeida-Thouless (AT) line [46]. It immediately means divergence of the so called spin-glass susceptibility negative divergence of non-linear compressibility d 2 p/dδ 2 much as the spinglass transition of the SK model [52,53]. Beyond the transition point, going into the glass phase, we have to consider breaking of the replica symmetry. [43,54,55].

RSB solution
Finally let us study the glass phase of the p = 2 hardcore model using the RSB ansatz. Here we use the softcore potential Eq. (13) and extend the analysis to finite temperatures. Using Eq. (139) and evaluating the integral numerically using the softcore potential we can easily find the plane α = α c (δ, T ) where the AT instability λ R = 0 occurs. The result is shown in Fig. 8 where the AT plane is indicated as 'AT'. The zero temperature limit of it agrees with the AT line of the hardcore model shown in Fig. 6 as it should be. The AT plane separates the liquid phase (paramagnet) with q = 0 on its left hand side and glass the glass phase on the right hand side.
Next let us examined the 1RSB ansatz on the glass side of the AT plane. We solved the 1RSB equations numerically following the scheme explained in sec. VIII E 2 to obtain q 1 assuming q 0 = 0. Note that q 0 = 0 always solves the saddple point equation for p > 1 as we mentioned in sec. VIII E. We found q 1 emerges continuously starting from the AT plane, as expected. Then we evaluated the complexity Σ * (m 1 ) numerically (see sec. VIII B 2) and determined m 1 where the complexity vanishes. We examined the stability of the 1RSB solution by evaluating the replicon eigenvalue λ R of the 1RSB solution Eq. (213). As shown in Fig. 8. we have two planes indicated as 'G1' and 'G2' where the replicon eigenvalue vanishes suggesting the Gardner's transition [45]. We found 'G1' plane merges with the AT plane at higher temperatures as can be seen in the figure. The 1RSB solution is stable above these Gardner planes but become unstable below them.
Below the Gardner planes and on the right on the glass side of the AT plane we naturally expect continuous RSB. Indeed we obtain the full RSB solution (approximated by k = 200 RSB) as shown in Fig. 9 where we show some examples of the q(x) functions obtained at T = 0 (hardcore limit). We obtained the result using the scheme explained in sec. VIII E 3 together with the inputs for the hardcore case shown in sec. X B 1. As can be seen in Fig. 9 a), we found the continuous RSB solution with nonzero q(x) function emerges continuously starting from the AT line α = α c (δ) as expected.
Using the scheme explain in sec. X B 4 we obtained the jamming line α = α J (δ) of the hardcore model and the result is displayed in Fig. 8 at the bottom. It is also shown in Fig. 6 where we can see that the RS ansatz overestimates the jamming line. Quite interestingly we find in Fig. 8 that the two Gardner planes 'G1' and 'G2' merges on the the jamming line in the zero temperature limit. The geometry of the phase diagram is very different from that of the hardspheres [56] where there is only one Gardner line. We analyzed the criticality of jamming q(1) → 1 of the hardcore model ( → ∞), i. e. α → α − J (δ) at T = 0. As shown in Fig. 9 b), we find power law behavior 1 − q(x) ∝ x −κ with κ = 1.4157... This confirms the scaling argument presented in sec. X C 1 establishing that the present model belong to the same universality class as the hardspheres [11].
The liquid phase q = 0 at T = 0 can be regarded as an easy SAT region where the space of the solutions to satisfy the hard constraints ( → ∞), i. e. the manifold of the ground states are continuously connected. The glass phase at T = 0 with q(1) < 1 can be regarded as hard SAT phase where the manifold of the ground states splits into clusters. The major difference with respect to the case of usual discrete coloring [26] is that the transition is continuous and the clustering is hierarchical reflecting the continuous RSB. The region α > α J (δ) is the UNSAT region where the constrain cannot be satisfied.  The straight line represents the power law fit ax −κ with κ = 1.4157, the same exponent as that for the hardspheres [11].

XI. CONCLUSIONS
In the present paper developed a family of exactly solvable large M -component vectorial spin systems with pbody interactions which exhibit glass transitions by the self-generated randomness. We also established a connection between the disorder-free model and a completely disordered spin glass model by constructing a model which interpolates the two limits. We showed that the supercooled paramagnetic states and glassy states are locally stable against crystallization under certain conditions, namely either 1) p > 2 or 2) the interaction potential V (x) has a flatness. We developed a replica formalism to solve the problem exactly. We applied the scheme to study models with non-linear potentials: the quadratic and soft/hardcore potentials. With the non-linear potential, we found continuous RSB also with p = 2 continuous spins which does not exhibit RSB in the case of the linear potential.
In the hardcore model we found the system exhibits jamming. Using the continuous RSB solution we found the isostaticity holds at jamming and the universality of it turned out to be the same as hardspheres for all values of p establishing the superniversality. This observation extend the result on the perceptron [28,29] which corresponds to p = 1.
Although we limited ourselves to the case M → ∞ in the present paper, a systematic 1/M expansion is possible. Such an approach has been conducted in the case of p = 2 continuous spin model with the linear potential where RSB does not take place [57,58]. This would be an alternative, analytically tractable scheme to analyze systems on tree-like lattices of finite connectivity (Bethe lattice) where mean-field approach should remain valid. So far such systems remained hard to be analyzed by the the replica approach [59] so that the cavity approach is usually preferred which however is limited to 1RSB at the moment [60]. An advantage of the 1/M expansion approach is that one can analyze the system almost as easily as the globally coupled systems so that one can construct continuous RSB explicitly as we have done in M → ∞ limit, which may become necessary deep in the glassy phases [45]. It will be interesting to study, for example, how the nature of jamming change if M becomes finite. To what extent such mean-field results would remain useful for finite dimensional systems where the lattices are no more like trees, is of course remains as an important problem.
We expect our results provide a useful basis to formulate theoretical approaches to study glass transitions of rotational degree of freedoms. For example it is natural to study glass transitions of particulate systems with rotational degrees of freedom such as patchy colloids and ellipsoids (Fig. 1 c)) extending the present work . Another interesting problem is to study the apparently disorder-free spinglass transitions realized in some frustrated magnets [23,24] (Fig. 1 b)). Continuous constrained satisfaction problems such as the continuous coloring problem (Fig. 1 a)) and related statistical inference problems can also be studied in a similar framework.
The free-energy F is obtained by minimizing the variational free-energy functional F[ρ(S)]. The 1st term on the r.h.s of Eq. (A5) represents the entropic (paramagnetic) part of the free-energy. The 2nd term is the 1st virial correction due to interactions. The reason for the absence of the higher order terms, all of which are represented as 1 particle irreducible (1PI) diagrams such as a triangle, a square, e.t.c. [34], is the tree-like geometry of the lattices that we consider.

Replicated spin liquid
In principle all stable and metastable states of the system, including liquid (paramagnetic) state ρ liq (S), crystalline state ρ crystal (S) and glassy states ρ α (S) (α = 1, 2, . . . ,), would be found as local minima of the free-energy functional Eq. (A5). In the present paper we focus on the properties of glassy states which emerge from supercooled paramagnetic state.
A useful way to analyze the properties of glassy states is the replica approach. We consider replicated spin liquid of n replicas labeled as a = 1, 2, . . . , n obeying the Hamiltonian,

Glass order parameter functional
We look for glassy metastable states which keep the statistical rotational invariance of the liquid (paramagnetic) state. To this end it is natural to consider the overlap matrix Eq. (69) as the order parameter, for a, b = 1, 2, . . . , n. Note that Q aa = 1 due to the normalization of the spins |S a i | 2 = M .

a. Variational free-energy
Based on the above discussion we expect that ρ(S) of the the glassy states, which keeps the statistical rotational invariance of the liquid, is parametrized solely by the overlap matrixQ, Similarly we anticipate that the replicated Mayer function can be parametrized as, where dQ = a<b dQ ab . Here J(Q) is the Jacobian (see below) and the parameter λ in the last term of Eq. (A15) is a Lagrange multiplier to enforce the normalization of the spin density. Note that in the 2nd integral on the r. h. s. of Eq. (A15) we assumed a simply factorized Jacobian p l=1 J(Q) disregarding possible cross-correlations of spins at different sites l = 1, 2, . . . , p. We comment on the validity of this assumption later.
The Jacobian J is defined as Here we used the fact that in the M → ∞ limit, different components of the spins S µ become independent from each other. This can be checked by introducing integral representation of the δ function in Eq. (A34) which can be evaluated by the saddle point method in M → ∞ limit.
To sum up we find the replicated Mayer function in the M → ∞ limit as, The last equation is obtained by repeating integrations by parts.
Collecting the above results we find the thermodynamic free-energy Eq. (A9)as, with the variational free-energy (more precisely free-entropy), for all a = b.

c. Gaussian ansatz
Finally let us note that one can check that the above result can be reproduced by assuming an Gaussian ansatz for the replicated spin density, The situation is essentially the same as that of hardspheres in large dimensional limit [9][10][11].
Appendix B: Eigenvalues of the stability matrix Here we analyze the Hessian matrix M a =b,c =d of the free-energy around the saddle points. It is a matrix of size n(n − 1) × n(n − 1) defined as,