# The Replica Symmetry Broken States of some Glass Models

J. Yeo  Department of Physics, Konkuk University, Seoul 05029, Korea  M. A.
Moore  Department of Physics and Astronomy, University of Manchester,
Manchester M13 9PL, United Kingdom.

###### Abstract

We have studied in detail the $M$-$p$ balanced spin glass model, which is a
candidate for being a model for structural glasses. Such models possess two
kinds of broken replica states; those with one-step replica symmetry breaking
(1RSB) and those with full replica symmetry breaking (FRSB). To determine
which arises requires studying the Landau expansion to quintic order. There
are 9 quintic order coefficients, and 5 quartic order coefficients, whose
values we determine for this model. We show that it is only for $2\leq
M<2.4714\cdots$ that the transition at mean-field level is to a state with
FRSB, while for larger $M$ values there is either a continuous transition to a
state with 1RSB (when $M\leq 3$) or a discontinuous transition for $M>3$. The
Gardner transition from a 1RSB state at low temperatures to a state with FRSB
also requires the Landau expansion to be taken to quintic order. Our result
for the form of FRSB in the Gardner phase is similar to that found when $2\leq
M<2.4714\cdots$, but differs from that given in the early paper of Gross et
al. [Phys. Rev. Lett. 55, 304 (1985)]. Finally we discuss the effects of
fluctuations on our mean-field solutions using the scheme of Höller and Read
[Phys. Rev. E 101, 042114 (2020)] and argue that such fluctuations will remove
the continuous 1RSB transition in dimension $d$ when $8>d\geq 6$ leaving just
the FRSB continuous transition (and possibly also the discontinuous 1RSB
transition). We suggest values for $M$ and $p$ which might be used in
simulations to resolve the outstanding question of whether fluctuation
corrections can remove the discontinuous 1RSB transition.

##  I Introduction

Spin models of the $p$-spin or Potts glass variety Gross  _et al._ (1985);
Gardner (1985) played an important role in the development of one of the
current theories of structural glasses, the Random First Order Transition
(RFOT) picture Kirkpatrick  _et al._ (1989); Kirkpatrick and Thirumalai
(2015); Lubchenko and Wolynes (2007); Cavagna (2009); Biroli and Bouchaud
(2010). These models have been primarily studied in the infinite
dimensionality limit, which is equivalent to mean-field theory. Of course what
is really wanted is an understanding of what happens in the physical realm of
two and three dimensions, and for these dimensions simulations Franz and
Parisi (1999); Campellone  _et al._ (1998) of models of the type studied in
this paper have revealed that they behave completely differently from what is
predicted by the mean-field calculations. In particular in the simulations
there is no sign of the random first-order transition which is one of the
central features of RFOT theory. Below the ideal glass transition there is
supposed to exist the ideal glass state, a state of low configurational
entropy but with a high stability due to the assumed paucity of glass states.
This state in replica language has one-step replica symmetry breaking (1RSB).
The transition temperature to this state is identified as the Kauzmann
temperature in RFOT theory, which is the temperature at which the entropy of
the glass state becomes equal to that of the crystalline state Kauzmann
(1948). While a discontinuous transition was not seen in the simulations,
evidence was found for the existence of long correlation lengths, which is
also the behavior found in real-space renormalization group (RG) calculations
Yeo and Moore (2012a, b) of $p$-spin models in three dimensions.

That simulations in three dimensions lead to a picture quite different to that
which arises from mean-field calculations has largely been ignored: Work has
continued apace using the large $d$ limit and mean-field techniques. We have
therefore begun a program of trying to understand why the mean-field picture
does not extend to three dimensions Yeo and Moore (2020). For one particular
$p$-spin model, the $M$-$p$ spin glass model with $p=6$, we were able to give
an argument that the 1RSB state of that model was unstable in any finite
dimension due to the excitation of droplets of flipped spins whose interface
free energy are very small Moore (2006). That argument is specific to glass
models with a particular form of time reversal symmetry which gives rise to a
field theory in which the cubic term $w_{2}$ is zero (see Eq. (25)).
Unfortunately the generic field theories thought relevant to glasses have
$w_{2}$ non-zero and it is these which we study in this paper. The 1RSB phase
for $p=6$ spin glasses is destroyed by non-perturbative droplet excitations.
For generic glass models with $w_{2}$ non-zero, we can only find perturbative
arguments. They are strong enough to lead us to the conclusion that the
continuous phase transition to a state with 1RSB will not exist for dimensions
$d$ less than $8$ and will be replaced by a continuous transition to a state
with FRSB. We shall suggest that fluctuation corrections to the coupling terms
in Eq. (25) might also drive the system away from having a discontinuous
transition to a 1RSB state to a continuous transition to a state with full
replica symmetry breaking (FRSB), but we do not know whether the fluctuation
corrections are large enough to bring that about. We suspect that this
question will only be resolved by simulations and values of $p$ and $M$ which
might be appropriate for such simulations are suggested in Sec. III.

Our procedure is based upon the old idea Rudnick and Nelson (1976) of using
the renormalization group recursion relations for the coupling constants of
the field theory to map the coefficients of the critical field theory into a
region where the correlation lengths are small and Landau theory (i.e. mean-
field theory) with small fluctuation corrections can be employed. This program
has also been used by Höller and Read Höller and Read (2020) on the problem of
the de Almeida-Thouless transition of the Ising spin glass in a field de
Almeida and Thouless (1978). It has a field theory identical to that of the
$M$-$p$-spin glass models discussed in this paper, i.e. that of Eq. (25), but
with different numerical values for the coefficients. (To discuss finite
dimensions a gradient term of the form $\int d^{d}r\sum_{a,b}(\nabla
q_{ab}(r))^{2}$ would need to be included in Eq. (25).) The program therefore
requires us to understand in detail the stationary solutions i.e. mean-field
solutions of Eq. (25), and the bulk of this paper is devoted to this task.
Because Höller and Read discussed the RG aspects of the calculations in great
detail, we shall treat those briefly, just focussing on the implications of
numerical studies which were carried out after their paper was written
Aguilar-Janita  _et al._ (2023).

In Sec. II we introduce the balanced $M$-$p$ models and the replica procedure
which was used to average their free energy over disorder. The balanced
$M$-$p$ spin models are very convenient to study with simulations as they are
readily extended to finite dimensions on a $d$-dimensional lattice. When this
is done the resulting field theory acquires the already mentioned gradient
squared term. One of the attractions of the balanced version of these models
is the absence of “hard modes”, which are just usually cast aside (as in the
paper of Caltagirone et al. Caltagirone  _et al._ (2011)), but this leaves the
subsequent calculations of uncertain accuracy. We shall focus on the case
$p=4$ and regard the number of types of Ising spins $M$ as a variable which
can take non-integer values. The simulations of Campellone et al. Campellone
_et al._ (1998) which failed to find a discontinuous 1RSB transition were in
fact done for a closely related model with $p=4$ and $M=4$ in three
dimensions. At cubic order there are two coupling constants, $w_{1}$ and
$w_{2}$, at quartic order, there are five coupling constants,
$y_{1},\cdots,y_{5}$ and at quintic order, there are nine coupling constants,
$z_{1},\cdots,z_{9}$. The quadratic term $\tau$ vanishes as usual at the mean-
field transition temperature $T_{c}$ and is negative when $T<T_{c}$. We
calculate the “bare” value of all these coefficients in Appendix A for the
case $p=4$. Fluctuation corrections will modify the bare values. In studying
the model at non-integer values of $M$ we are anticipating that the
fluctuation corrections can modify the bare coefficients. Studying the field
theory of Eq. (25) for general values of the coefficients would be a good
idea, but there are so many of these coefficients that we have limited our
study to those values which can be reached by varying $M$ in the bare values.
In Sec. III we discuss what we believe will be the likely consequences of
fluctuation effects on the coupling constants.

In Sec. II.1 we determine the free energy of the system in the high-
temperature or paramagnetic phase where the order parameter $q_{ab}$ is
independent of $a$ and $b$, that is, replica symmetric. At mean-field level
$q_{ab}=0$, (but fluctuation corrections would leave it replica symmetric but
non-zero). If the transition is continuous, so that $q_{ab}$ is small just
below the transition, then the expansion of the Landau-Ginzburg free energy
functional in powers of $q_{ab}$ should be useful and we give its form in Sec.
II.2. Most workers have stopped at the quartic terms, but we have continued up
to the quintic terms. This is necessary for two reasons. The difference in
free energy between the 1RSB free energy and the FRSB free energy is of
$O(\tau^{5})$, (see for example, Ref. Aspelmeier  _et al._ (2008)). Thus one
needs to worry about the quintic terms when working out whether the state
which forms at the continuous transition is of 1RSB type or is of FRSB type.
Fortunately, we can show that the borderline value of $M$, $M^{**}\approx
2.47140$ between these types is not dependent on the quintic terms. (For
$2\leq M<M^{**}$ the continuous transition is to a state with FRSB, while for
$M^{**}<M<3$, the continuous transition is to a state with 1RSB.) The second
reason relates to studies of the Gardner transition Gross  _et al._ (1985);
Gardner (1985). The Gardner transition is the transition from a state with
1RSB to a state with FRSB as the temperature is lowered. Right from the
beginning it was realized that the quintic terms are needed for its study
Gross  _et al._ (1985). We shall find though that our actual FRSB solution is
quite different to that of Ref. Gross  _et al._ (1985). This is discussed in
Sec. II.5.

A feature of the FRSB solutions is a singularity first noticed by Goldbart and
Elderfield Goldbart and Elderfield (1985). They found that the FRSB solution
for $q(x)$ at quartic level could have an unphysical singularity in the
interval $0<x<1$ which would imply that the probability of two states having
an overlap $q$ would be negative, which is impossible. This problem was
studied in some detail by Janiš and colleagues using a non-standard approach
to replica symmetry breaking Janiš  _et al._ (2013). We find in Sec. II.5 that
the singularity at quartic level in fact determines the value of $M^{**}$ and
that one avoids the singularity at $M>M^{**}$ by simply being in the state
with 1RSB. At the Gardner transition the quintic terms remove the quartic
level singularities. However, similar singularities are to be found also at
quintic level. Right at the Gardner transition temperature $T_{G}$ , just
where the free energies of the FRSB state and the 1RSB state are equal, the
Goldbart-Elderfield singularity is at the lower breakpoint $x_{1}$. This
causes the derivative of $q(x)$ at $x=x_{1}$ to be infinite. However for
temperatures $T$ less than $T_{G}$, the singularity is below $x_{1}$ and the
derivative stays finite.

In Sec. II.3 we derive the free energy at mean-field level for the 1RSB state.
For $M>3$, when $w_{2}/w_{1}>1$, the transition from the high temperature
normal phase to a state with 1RSB is a discontinuous transition which takes
place at a transition temperature above $T_{c}$. We suspect that this behavior
would be seen for all values of $M>3$. However, if one truncates the free
energy to quartic level terms, as is commonly done, the 1RSB state only exists
in the interval $3<M\lesssim 6.64$. With the inclusion of the quintic terms,
the 1RSB forms at a discontinuous transition when $14.41\gtrsim M\gtrsim 3.98$
and $3.27\gtrsim M>3$. Thus with the quintic form the 1RSB state persists up
to larger values of $M$. We believe that if all terms were kept then the
discontinuous transition to the 1RSB state would exist for all $M>3$. In Sec.
II.4 we describe the simplifications which arise in the large $M$ limit.
Truncation leads to spurious features as the Landau expansion cannot be
expected to be accurate when $q_{ab}$ is not small. Another spurious feature
of truncation is the apparent phase transition at low temperatures from the
1RSB state to the replica symmetric state with $q_{ab}$ non-zero. In the large
$M$ limit we can solve without truncation and such a transition does not arise
(see Sec. II.3).

The form of the FRSB solutions at both quartic and quintic level, together
with the Gardner transitions, is in Sec. II.5. In Sec. III we discuss how
fluctuation corrections to the coupling constants used in the mean-field
solution will change the continuous 1RSB transition into the continuous FRSB
solution, using extensions of the approach of Höller and Read Höller and Read
(2020). We suspect that the discontinuous 1RSB transition might also suffer
the same fate, based on the results of simulations in low dimensions Franz and
Parisi (1999); Campellone  _et al._ (1998), but we cannot support this
possibility with analytical arguments. We finally conclude with suggestions of
the kinds of model which could be studied numerically to resolve these issues,
and also to resolve the question of whether the FRSB state can exist for
dimensions $d<6$.

##  II The balanced $M$-$p$ model in the fully connected limit

In this section, we study the $M$-$p$ spin glass model in the fully connected
limit, where one has $M$ different types of Ising spins, $S_{i}(x)$,
$i=1,2,\cdots,M$ at each site $x$ coupled with spins on other sites via
$p$-body interactions. Here we focus on the so-called balanced model
introduced in Ref. Yeo and Moore (2020) for even $p$, where only the coupling
between two sets of $p/2$ spins on two different sites is considered. It
amounts to considering only the soft mode in a more general $M$-$p$ model,
where all the couplings between $k$ spins and $p-k$ spins are included for
$k=1,2,\cdots,p-1$.

In this paper, we focus on the $p=4$ case. For $p=4$, the balanced model is
given by four-spin interactions between a pair of two spins on two different
sites. Each site has ${M\choose 2}$ different two-spin combinations.
Therefore, for given pair of sites, there are ${M\choose 2}^{2}$ terms in the
Hamiltonian. The Hamiltonian is given by

$\displaystyle H=-\frac{1}{2}\sum_{x\neq y}$
$\displaystyle\Big{[}\sum_{i_{1}<i_{2}}^{M}\sum_{j_{1}<j_{2}}^{M}J^{(i_{1},i_{2}),(j_{1},j_{2})}_{x,y}$
$\displaystyle\times S_{i_{1}}(x)S_{i_{2}}(x)S_{j_{1}}(y)S_{j_{2}}(y)\Big{]},$
(1)

where each $J^{(i_{1},i_{2}),(j_{1},j_{2})}_{x,y}$ is drawn from the Gaussian
distribution with zero mean and the variance

$\frac{J^{2}}{NM^{p-1}}=\frac{J^{2}}{NM^{3}}.$ (2)

We will set $J=1$ for convenience. After neglecting the terms of subleading
order in $N$, we can write the replicated partition function averaged over the
disorder as

$\displaystyle\overline{Z^{n}}=$
$\displaystyle\mathrm{Tr}\exp\Big{[}\frac{\beta^{2}}{4NM^{3}}$ (3)
$\displaystyle\times\sum^{n}_{a,b}\Big{\\{}\sum^{N}_{x}\sum_{i_{1}<i_{2}}^{M}S^{a}_{i_{1}}(x)S^{a}_{i_{2}}(x)S^{b}_{i_{1}}(x)S^{b}_{i_{2}}(x)\Big{\\}}^{2}\Big{]}.$

The diagonal terms ($a=b$) in the replica indices give a factor
$\exp[nN\beta^{2}C]$ where

$\displaystyle C=\frac{1}{4M^{3}}{M\choose 2}^{2}=\frac{(M-1)^{2}}{16M}.$ (4)

For $a\neq b$, following the convention used in Ref. Caltagirone  _et al._
(2011), we introduce the delta functions enforcing

$\displaystyle
q_{ab}=\frac{1}{NM^{2}}\sum^{N}_{x}\sum_{i_{1}<i_{2}}^{M}S^{a}_{i_{1}}(x)S^{a}_{i_{2}}(x)S^{b}_{i_{1}}(x)S^{b}_{i_{2}}(x)$
(5)

in the replicated partition function. Using the integral representation of the
delta function, we can write

$\displaystyle\overline{Z^{n}}=e^{nN\beta^{2}C}\int\prod_{a<b}dq_{ab}d\mu_{ab}\;\exp[-NG(\underline{q},\underline{\mu})],$
(6)

where

$G(\underline{q},\underline{\mu})=-\frac{M}{4}\beta^{2}\sum_{a\neq
b}q^{2}_{ab}+\frac{M}{2}\sum_{a\neq b}\mu_{ab}q_{ab}-\ln L(\underline{\mu})$
(7)

and

$\displaystyle
L(\underline{\mu})=\underset{\\{S_{i}^{a}\\}}{\mathrm{Tr}}\exp\Big{[}\frac{1}{2M}\sum_{a\neq
b}\mu_{ab}\sum^{M}_{i<j}S^{a}_{i}S^{a}_{j}S^{b}_{i}S^{b}_{j}\Big{]}.$ (8)

In the large-$N$ limit, the integral is dominated by the saddle points which
are determined by

$\displaystyle\mu_{ab}=\beta^{2}q_{ab}$ (9)

and

$\displaystyle
q_{ab}=\frac{1}{M^{2}}\left\langle\sum^{M}_{i<j}S^{a}_{i}S^{a}_{j}S^{b}_{i}S^{b}_{j}\right\rangle_{L},$
(10)

where $\langle\cdots\rangle_{L}$ is evaluated with respect to $L$ in Eq. (8).
The free energy $F$ is then given by

$\displaystyle\frac{\beta F}{N}$ $\displaystyle=-\frac{1}{N}\lim_{n\to
0}\frac{1}{n}\ln\overline{Z^{n}}=-C\beta^{2}+\lim_{n\to
0}\frac{1}{n}G(\underline{q},\underline{\mu}).$ (11)

###  II.1 Replica Symmetric Solution

We first look for the saddle point solutions in the replica symmetric (RS)
form $q_{ab}=q$ and $\mu_{ab}=\mu$ for all $a\neq b$. We have

$\displaystyle\lim_{n\to
0}\frac{1}{n}G(q,\mu)=\frac{M}{4}\beta^{2}q^{2}-\frac{M}{2}\mu q-\lim_{n\to
0}\frac{1}{n}\ln L(\mu).$ (12)

Using

$\displaystyle\sum_{a\neq
b}S^{a}_{i}S^{a}_{j}S^{b}_{i}S^{b}_{j}=\left(\sum_{a}S^{a}_{i}S^{a}_{j}\right)^{2}-n$
(13)

in Eq. (8) and the Hubbard-Stratonivich transformation on the first term, we
can rewrite Eq. (8) as

$\displaystyle L(\mu)=$ $\displaystyle
e^{-n\mu(K/2M)}\;\underset{\\{S_{i}^{a}\\}}{\mathrm{Tr}}\;\int D^{K}\bm{y}$
$\displaystyle\times\exp\left[\sqrt{\frac{\mu}{M}}\sum_{a}\sum_{i<j}^{M}y_{(i,j)}S^{a}_{i}S^{a}_{j}\right],$
(14)

where

$K\equiv{M\choose 2}$ (15)

and the integral over the $K$-dimensional vector
$\bm{y}=(y_{1},y_{2},\cdots,y_{K})\equiv(y_{(1,2)},y_{(1,3)},\cdots,y_{(M-1,M)})$
is defined as

$\int
D^{K}\bm{y}\equiv\prod_{\alpha=1}^{K}\left(\int_{-\infty}^{\infty}\frac{dy_{\alpha}}{\sqrt{2\pi}}e^{-y^{2}_{\alpha}/2}\right).$
(16)

We therefore have

$\displaystyle\lim_{n\to 0}\frac{1}{n}\ln L(\mu)=-\frac{K}{2M}\mu+M\ln 2+\int
D^{K}\bm{y}\;\ln\zeta(\bm{y},\mu),$ (17)

where

$\displaystyle\zeta(\bm{y},\mu)\equiv\frac{1}{2^{M}}\underset{\\{S_{i}\\}}{\mathrm{Tr}}\;\exp\left[\sqrt{\frac{\mu}{M}}\bm{y}\cdot\bm{\Psi}\right]$
(18)

with the $K$-dimensional vector
$\bm{\Psi}=(\Psi_{1},\Psi_{2},\cdots,\Psi_{K})=(S_{1}S_{2},S_{1}S_{3},\cdots,S_{M-1}S_{M})$.
The RS free energy is then given by

$\displaystyle\frac{\beta F_{\rm RS}}{N}=$
$\displaystyle-C\beta^{2}+\frac{M}{4}\beta^{2}q^{2}-\frac{M}{2}\mu
q+\frac{K}{2M}\mu$ $\displaystyle-M\ln 2-\int
D^{K}\bm{y}\;\ln\zeta(\bm{y},\mu).$ (19)

By varying the free energy with respect to $q$ and $\mu$, respectively, we
have saddle point equations,

$\displaystyle\mu=\beta^{2}q$ (20)

and

$\displaystyle q=$ $\displaystyle\frac{1}{M^{2}}\int
D^{K}\bm{y}\;\frac{1}{\zeta^{2}(\bm{y},\mu)}$
$\displaystyle\times\sum_{\alpha=1}^{K}\left\\{\frac{1}{2^{M}}\underset{\\{S_{i}\\}}{\mathrm{Tr}}\;\Psi_{\alpha}\exp\left[\sqrt{\frac{\mu}{M}}\bm{y}\cdot\bm{\Psi}\right]\right\\}^{2}.$
(21)

At high temperatures, the RS solutions are given by $q=\mu=0$. In that case,
$\zeta=1$ and the corresponding free energy is

$\frac{\beta F_{\rm RS}}{N}=-C\beta^{2}-M\ln 2.$ (22)

The entropy $S=-\partial F/\partial T$ for this phase is

$\frac{S_{\rm RS}}{N}=-C\beta^{2}+M\ln 2.$ (23)

This becomes negative below

$T_{*}=\sqrt{\frac{C}{M\ln 2}}=\frac{M-1}{4M\sqrt{\ln 2}}.$ (24)

Some values of $T_{*}$ are $T_{*}$=0.20019 for $M=3$, 0.22521 for $M=4$,
0.25023 for $M=6$ and 0.25738 for $M=7$. It keeps increasing with $M$ and
approaches 0.30028 in the $M\to\infty$ limit.

###  II.2 Landau Expansion of Free Energy

In order to study a possible continuous transition, we expand the free energy,
Eq. (11) for small values of the order parameter. We first expand Eq. (8) to
$O(\mu^{5})$ and take the trace over the spins. The detailed steps are given
in Appendix A. Now using Eqs. (7), (9) and (11), we can write the free energy
as

$\displaystyle\frac{\beta F}{N}$ $\displaystyle=-C\beta^{2}-M\ln 2+\lim_{n\to
0}\frac{1}{n}\Big{[}\tau\sum_{a,b}q^{2}_{ab}$ (25) $\displaystyle-
w_{1}\sum_{a,b,c}q_{ab}q_{bc}q_{ca}-w_{2}\sum_{a,b}q^{3}_{ab}-y_{1}\sum_{a,b}q^{4}_{ab}$
$\displaystyle-
y_{2}\sum_{a,b,c}q^{2}_{ab}q^{2}_{bc}-y_{3}\sum_{a,b,c}q^{2}_{ab}q_{bc}q_{ca}-y_{5}\sum_{a,b,c,d}q_{ab}q_{bc}q_{cd}q_{da}$
$\displaystyle-
z_{1}\sum_{a,b}q^{5}_{ab}-z_{2}\sum_{a,b,c}q^{3}_{ab}q^{2}_{bc}-z_{3}\sum_{a,b,c}q^{3}_{ab}q_{bc}q_{ca}$
$\displaystyle-
z_{4}\sum_{a,b,c}q^{2}_{ab}q^{2}_{bc}q_{ca}-z_{5}\sum_{a,b,c,d}q^{2}_{ab}q_{bc}q_{cd}q_{da}$
$\displaystyle-
z_{6}\sum_{a,b,c,d}q^{2}_{ab}q_{bc}q_{cd}q_{db}-z_{7}\sum_{a,b,c,d}q^{2}_{ab}q_{bc}q^{2}_{cd}$
$\displaystyle-
z_{8}\sum_{a,b,c,d}q_{ab}q_{bc}q_{cd}q_{da}q_{ac}-z_{9}\sum_{a,b,c,d,e}q_{ab}q_{bc}q_{cd}q_{de}q_{ea}\Big{]},$

where $q_{aa}=0$, $q_{ab}=q_{ba}$, and all the sums over replica indices are
without any restriction.

The coefficient of the quadratic term is given by

$\displaystyle\tau=\frac{M}{4}\beta^{2}\left(1-\frac{K}{M^{3}}\beta^{2}\right)=\frac{M}{4}\beta^{4}\left(T^{2}-T^{2}_{c}\right),$
(26)

where

$T_{c}\equiv\sqrt{\frac{K}{M^{3}}}=\frac{1}{M}\sqrt{\frac{M-1}{2}}.$ (27)

This expression coincides with Eq. (27) of Ref. Caltagirone  _et al._ (2011).
Some values of $T_{c}$ are 0.33333 for $M=3$, 0.30619 for $M=4$, 0.26352 for
$M=6$ and 0.24744 for $M=7$. Note that $T_{c}$ decreases with $M$ and becomes
zero in the $M\to\infty$ limit. Note also that $T_{c}>T_{*}$ for
$M=2,3,\cdots,6$ and $T_{*}>T_{c}$ for $M\geq 7$.

The coefficients of the cubic terms are given by

$\displaystyle
w_{1}=\frac{\beta^{6}K}{6M^{3}},~{}~{}w_{2}=\frac{\beta^{6}K}{6M^{3}}(M-2).$
(28)

The quartic and quintic coefficients are given in Appendix A as functions of
$M$. It is known Gross  _et al._ (1985); Caltagirone  _et al._ (2011) that if
the ratio of the cubic terms $w_{2}/w_{1}$, which in our model is equal to
$M-2$, is greater than one, a discontinuous transition to the one-step replica
symmetry breaking phase (1RSB) occurs. When $M=2$, our model reduces to the
Ising spin glass and we can check that the cubic and quartic coefficients
coincide with those for the Ising spin glass except for the multiplicity
factor of $2^{3}$ for $w_{i}$ and $2^{4}$ for $y_{i}$.

###  II.3 The 1RSB Solution

We now consider the case where $q_{ab}$ and $\mu_{ab}$ take the one step
replica symmetry breaking (1RSB) form taking values $q_{1}$ and $\mu_{1}$ on
$n/m_{1}$ diagonal blocks (labelled by $B_{k}$, $k=1,2,\cdots,n/m_{1}$ of size
$m_{1}$ and $q_{0}$ and $\mu_{0}$ outside the blocks. We then have the terms
in Eq. (7) as

$\displaystyle\sum_{a\neq
b}q^{2}_{ab}=n[(m_{1}-1)q^{2}_{1}+(n-m_{1})q^{2}_{0}],$ (29)
$\displaystyle\sum_{a\neq
b}\mu_{ab}q_{ab}=n[(m_{1}-1)\mu_{1}q_{1}+(n-m_{1})\mu_{0}q_{0}].$ (30)

We will focus on the 1RSB solutions with $q_{0}=\mu_{0}=0$. By writing

$\displaystyle\frac{1}{2M}\sum^{M}_{i<j}\sum_{a\neq
b}\mu_{ab}S^{a}_{i}S^{a}_{j}S^{b}_{i}S^{b}_{j}$
$\displaystyle=\frac{\mu_{1}}{2M}\sum_{k=1}^{n/m_{1}}\sum^{M}_{i<j}\left\\{\left[\sum_{a\in
B_{k}}S^{a}_{i}S^{a}_{j}\right]^{2}-m_{1}\right\\}$ (31)

in Eq. (8) and by using the Hubbard-Stratonovich transformation, we have

$\displaystyle\underset{\\{S_{i}^{a}\\}}{\mathrm{Tr}}\,\exp\left[\frac{1}{2M}\sum^{M}_{i<j}\sum_{a\neq
b}\mu_{ab}S^{a}_{i}S^{a}_{j}S^{b}_{i}S^{b}_{j}\right]$ (32) $\displaystyle=$
$\displaystyle\exp\left[-n\frac{\mu_{1}K}{2M}\right]$
$\displaystyle\times\Big{[}\int
D^{K}\bm{y}\;\Big{\\{}\underset{\\{S_{i}\\}}{\mathrm{Tr}}\,\exp\Big{[}\sqrt{\frac{\mu_{1}}{M}}\sum_{i<j}^{M}y_{(i,j)}S_{i}S_{j}\Big{]}\Big{\\}}^{m_{1}}\Big{]}^{n/m_{1}}.$

Therefore we have

$\displaystyle\lim_{n\to 0}\frac{1}{n}\ln L(\underline{\mu})=$
$\displaystyle-\frac{K}{2M}\mu_{1}+M\ln 2$
$\displaystyle+\frac{1}{m_{1}}\ln\int
D^{K}\bm{y}\;\zeta^{m_{1}}(\bm{y},\mu_{1}),$ (33)

where $\zeta$ is defined in Eq. (18). Using Eqs. (29), (30) and (33) in Eq.
(11),

$\displaystyle\frac{\beta F_{\rm 1RSB}}{N}=$
$\displaystyle-C\beta^{2}-\frac{M}{4}\beta^{2}(m_{1}-1)q^{2}_{1}$
$\displaystyle+\frac{M}{2}(m_{1}-1)\mu_{1}q_{1}+\frac{K}{2M}\mu_{1}-M\ln 2$
$\displaystyle-\frac{1}{m_{1}}\ln\int
D^{K}\bm{y}\;\zeta^{m_{1}}(\bm{y},\mu_{1}).$ (34)

Varying the free energy with respect to $q_{1}$ and $\mu_{1}$, respectively,
we have

$\mu_{1}=\beta^{2}q_{1}.$ (35)

and

$\displaystyle q_{1}$ $\displaystyle=\frac{1}{M^{2}}\frac{1}{\int
D^{K}\bm{y}\;\zeta^{m_{1}}(\bm{y},\mu_{1})}$ (36) $\displaystyle\times\int
D^{K}\bm{y}\;\zeta^{m_{1}-2}\sum_{\alpha=1}^{K}\left\\{\frac{1}{2^{M}}\underset{\\{S_{i}\\}}{\mathrm{Tr}}\;\Psi_{\alpha}\exp[\sqrt{\frac{\mu_{1}}{M}}\bm{y}\cdot\bm{\Psi}]\right\\}^{2},$

Now varying the free energy with respect to $m_{1}$, we have

$\displaystyle\frac{M}{4}\beta^{2}q^{2}_{1}+\frac{1}{m_{1}^{2}}\ln\int
D^{K}\bm{y}\;\zeta^{m_{1}}(\bm{y},\mu_{1})$
$\displaystyle-\frac{1}{m_{1}}\frac{\int
D^{K}\bm{y}\;\zeta^{m_{1}}(\bm{y},\mu_{1})\ln\zeta(\bm{y},\mu_{1})}{\int
D^{K}\bm{y}\;\zeta^{m_{1}}(\bm{y},\mu_{1})}=0.$ (37)

In summary, Eqs. (35), (36), and (37) are the saddle point equations one has
to solve for the 1RSB state.

Note that when $m_{1}=1$, we can explicitly evaluate

$\displaystyle\int
D^{K}\bm{y}\;\zeta(\bm{y},\mu_{1})=\exp\left[\frac{K}{2M}\mu_{1}\right].$ (38)

From Eq. (34), we see that when $m_{1}=1$, the 1RSB free energy is equal to
the RS one:

$\frac{\beta F_{\rm 1RSB}}{N}\underset{m_{1}\to
1}{\rightarrow}-C\beta^{2}-M\ln 2=\frac{\beta F_{\rm RS}}{N}.$ (39)

To determine the transition temperature $T_{c}^{\rm 1RSB}$ to the 1RSB state,
we set $m_{1}=1$ in Eqs. (35), (36) and (37) and solve for $\beta$. For
$m_{1}=1$, we can combine these three equations into one equation,
$f_{M}(\sigma)=0$ for the parameter

$\displaystyle\sigma\equiv\sqrt{\frac{\mu_{1}}{M}},$ (40)

where

$\displaystyle f_{M}(\sigma)$ $\displaystyle\equiv e^{-K\sigma^{2}/2}\int
D^{K}\bm{y}\;\Big{[}\zeta(\bm{y},\mu_{1})\ln\zeta(\bm{y},\mu_{1})$ (41)
$\displaystyle-\frac{\sigma^{2}}{4}\frac{\sum_{\alpha=1}^{K}\left\\{2^{-M}\mathrm{Tr}\;\Psi_{\alpha}\exp\left[\sigma\bm{y}\cdot\bm{\Psi}\right]\right\\}^{2}}{\zeta(\bm{y},\mu_{1})}\Big{]}-\frac{K}{2}\sigma^{2}.$

Note that $\zeta(\bm{y},\mu_{1})$ is a function of $\sigma$. If there exists a
nonzero solution $\sigma$ to $f_{M}(\sigma)=0$, one can obtain nonzero $q_{1}$
from Eq. (36) and the transition temperature $T_{c}^{\rm 1RSB}$ from Eq. (35).

Figure 1: $f_{M}(\sigma)$ defined in Eq. (41) for $M=3$. A nonzero solution
$\sigma$ of $f_{M}(\sigma)=0$ would signal a discontinuous transition into the
1RSB state.  Figure 2: Same as Fig. 1 with $M=4$.

We solve this equation by numerically evaluating multi-dimensional integrals
in Eq. (41). In Figs. 1 and 2, $f_{M}$ is plotted as a function of $\sigma$
for $M=3$ and $M=4$. As we can see from the figures, $f_{M}(\sigma)$ starts
off very flat and increases monotonically for large values of $\sigma$. For
$M=3$, Fig. 1 clearly shows a monotonic increase as a function of $\sigma$,
thus we can conclude that the only solution to $f_{3}(\sigma)=0$ is
$\sigma=0$. From Eq. (36), we then have $q_{1}=0$ thus no discontinuous
transition in this case. For $M=4$, we have to evaluate six-dimensional
($K=6$) integrals in Eq. (41). For that, we use Monte Carlo methods, and the
results are shown in Fig. 2. The error bars come from sampling random points
in the integrands within the Monte Carlo evaluation of the integrals. We have
averaged over 30 trials for each data point. Since $f_{4}(\sigma)$ stays very
flat for small $\sigma$ before increasing to large positive values, it is
quite difficult to determine, if any, nonzero solution $\sigma$ from this plot
alone.

To understand the situation more clearly, we study the behavior of
$f_{M}(\sigma)$ for small $\sigma$. We can show (see Appendix B for details)
that for small $\sigma$, the leading order in the small-$\sigma$ expansion of
$f_{M}(\sigma)$ is $O(\sigma^{6})$. In fact, if we write
$f_{M}(\sigma)=\sum_{i=0}^{\infty}c_{i}(M)\sigma^{i}$, we find that $c_{i}=0$
for $i$ odd, $c_{0}=c_{2}=c_{4}=0$ and

$\displaystyle c_{6}(M)=-\frac{M}{24}(M-1)(M-3)$ (42)

for $M=3,4,5,\cdots$. Therefore, for $M=3$, the leading order is actually
$O(\sigma^{8})$. The next-order coefficient is given by

$\displaystyle c_{8}(M)=-\frac{M}{48}(M-1)(3M^{2}-27M+47),$ (43)

for $M\geq 3$. Some steps needed to obtain these are given in Appendix B. We
note that $c_{8}(M=3)>0$. This is consistent with the monotonic increase of
$f_{3}(\sigma)$ shown in Fig. 1. For $M>3$, $c_{6}$ becomes negative.
Combining this fact with the monotonic increase for large $\sigma$, we can
conclude that there exists a nonzero solution to $f_{M}(\sigma)=0$ and that a
discontinuous transition for $M>3$ is expected. From Eq. (43), we find that
$c_{8}(M)>0$ for $M\lesssim 6.64$, therefore for these values of $M$, we can
estimate the solution as $\sigma\simeq\sqrt{-c_{6}(M)/c_{8}(M)}$. This
program, however, fails when $c_{8}(M)<0$ for $M\gtrsim 6.64$. ($c_{6}<0$ for
$M>3$.)

We need to go to higher order to study the 1RSB transition beyond this value
of $M$. We find, however, that the method in Appendix B becomes too cumbersome
to get $c_{10}$. The Landau expansion of the free energy given in Eq.(25)
provides a more useful tool. Since $\sigma^{2}\sim\mu_{1}\sim q_{1}$,
$O(\sigma^{6})$ and $O(\sigma^{8})$ correspond to the cubic and quartic orders
in $q_{ab}$, respectively, and we need quintic order terms in $q_{ab}$ to
evaluate $c_{10}$. In Appendix C, we apply the 1RSB form directly to $q_{ab}$
in Eq. (25). When $m_{1}=1$, the saddle point equations can be combined into a
form

$\displaystyle-\frac{1}{2}(w_{2}-w_{1})q^{3}_{1}-(y_{1}-y_{3}+y_{5})q_{1}^{4}-\frac{3}{2}z_{1}^{\rm
eff}q^{5}_{1}=0,$ (44)

where

$\displaystyle z_{1}^{\rm eff}\equiv z_{1}-z_{3}-z_{4}+z_{5}+z_{8}-z_{9}.$
(45)

Recalling that $q_{1}=\mu_{1}/\beta^{2}=M\sigma^{2}/\beta^{2}$ and using the
values of $w_{i}$ and $y_{i}$ given in Appendix A, we can identify the first
two terms in Eq. (44) as the small-$\sigma$ expansion of $f_{M}(\sigma)$,
since we can rewrite

$\displaystyle c_{6}(M)=-\frac{M^{3}}{2\beta^{6}}(w_{2}-w_{1}),$ (46)

and

$\displaystyle c_{8}(M)=-\frac{M^{4}}{\beta^{8}}(y_{1}-y_{3}+y_{5}).$ (47)

It follows that the last term in Eq. (44) gives

$\displaystyle c_{10}(M)=-\frac{3M^{5}}{2\beta^{10}}z_{1}^{\rm eff}.$ (48)

The explicit expression as a function of $M$ is given in Eqs. (C) and (144) in
Appendix C.

In Fig. 3, $(y_{1}-y_{3}+y_{5})/\beta^{8}$ and $z_{1}^{\rm eff}/\beta^{10}$
are displayed as functions of $M$. We note that $y_{1}-y_{3}+y_{5}$ is
negative (and $c_{8}$ is positive) for $2.35\lesssim M\lesssim 6.64$.
Therefore, as we mentioned above, we can find the 1RSB solution for $m_{1}=1$
for $3<M\lesssim 6.64$ within the quartic theory. The result for the 1RSB
transition temperature obtained in this way is shown as a solid red line in
Fig. 4 (a). We note, however, that the result becomes unreliable as we
approach the boundary value $M\simeq 6.64$ as it shows a fictitious diverging
behavior. We now study how the quintic theory may improve this result. The
quintic contribution can be summarized by $z_{1}^{\rm eff}$, which is negative
for $4.37\lesssim M\lesssim 12.46$ (and for the narrow region $2\leq M\lesssim
2.12$). Since $c_{10}$ is positive in that interval, we have a chance to
extend the result of the quartic theory to larger values of $M$. As one can
see in Fig. 4 (a), the 1RSB transition line calculated within the quintic
theory indeed extends to large values of $M$. But, since Eq. (44) for
$q_{1}\neq 0$ becomes a quadratic equation for $q_{1}$, there are intervals of
$M$ where no real solution exists. We find that for $3.27\lesssim M\lesssim
3.98$ and for $M\gtrsim 14.41$, solutions to this equation become complex and
no 1RSB solution can be obtained. This can be seen in Fig. 4 (b), where one
can see a segment of the 1RSB transition line is missing. Also as in the
quartic theory, the transition line displays an apparent divergent behavior as
we approach the boundary value $M\simeq 14.41$. Therefore, we can conclude
that it is possible to obtain the 1RSB transition line using truncated models,
but the truncation of the free energy to a specific order produces some
unphysical features. Comparing the results of the quartic and quintic theories
in Fig. 4 (a), we expect that a systematic improvement may occur if we go to
even higher orders. We also note that the 1RSB transition temperatures
obtained in this way always stay above $T_{*}$.

The 1RSB transition line discussed above is obtained by setting $m_{1}=1$
where the 1RSB free energy coincides with that of the high-temperature RS
phase (with $q=0$). Using the results in Appendix C, we can obtain 1RSB
solutions for general values of $0\leq m_{1}\leq 1$ for the truncated model.
Rather unexpectedly, we find that for given $M$, the 1RSB solution ceases to
exist below a certain finite temperature for which $m_{1}=0$. We note that if
$m_{1}=0$, the 1RSB free energy becomes that of the RS phase with nonzero $q$
(see Eq. (132)). Therefore, below that temperature, we only have the RS
solution with nonzero $q$. This is illustrated in Fig. 5, where we plot the
free energies of both 1RSB and RS solutions calculated within a truncated
model. One can clearly see that the 1RSB solution exists only in a finite
temperature interval. Within that interval, the system is in the 1RSB phase
which has a higher free energy than the RS one with nonzero $q$. However,
below that interval, there is no 1RSB solution, so the system returns to the
RS phase. We believe that this rather unusual behavior is caused by the
truncation of the model in an arbitrary order. In the large-$M$ limit
considered in Sec. II.4, where one can find the 1RSB solutions without
truncation, we find that the 1RSB solution continues down to zero temperature
and has a higher free energy than the RS one.

Figure 3: $(y_{1}-y_{3}+y_{5})/\beta^{8}$ (dashed line) and $z_{1}^{\rm
eff}/\beta^{10}$ (solid line) as functions of $M$. In the large-$M$ limit,
they approach 1/16 and 1/20, respectively.

Figure 4: (a) Red and blue solid lines are the 1RSB transition temperatures
$T_{c}^{\rm 1RSB}$ as functions of $M$ for the $p=4$ balanced $M$-$p$ model
expanded up to quartic (red) and to quintic (blue) orders in the order
parameter. Dashed and dot-dashed lines are $T_{*}$ (Eq. (24)) and $T_{c}$ (Eq.
(27)), respectively. Two closely-spaced horizontal lines are the large-$M$
limits of $T_{*}$ (lower one, Eq. (62)) and $T_{c}^{\rm 1RSB}$ (upper one, Eq.
(70)). (b) Close-up of the same plot for $3\leq M\leq 4$. There is a gap in
the solid blue line in the interval $3.27\lesssim M\lesssim 3.98$, where no
1RSB solution exists at $m_{1}=1$ for the quintic theory. The red line
corresponds to the quartic theory, which has no gap. The dot-dashed line is
$T_{c}$. Figure 5: Dimensionless free energies per spin of the 1RSB solution
(solid line) and the RS solution with $q\neq 0$ (dashed line) as functions of
temperature calculated for the quartic $M=4$ model. For each case, the free
energy difference ($\Delta F$) from that of the high-temperature RS solution
($q=0$, Eq. (22)) is plotted. The 1RSB solution exists only in the temperature
interval $0.212\leq T\leq 0.311$.

###  II.4 The Large-$M$ Limit

In this subsection, we consider the situation where we take the limit
$M\to\infty$ from the start. In the large-$M$ limit, Eq. (8) can be rewritten
as

$\displaystyle L(\underline{\mu})$
$\displaystyle=\underset{\\{S_{i}^{a}\\}}{\mathrm{Tr}}\exp\left[\frac{1}{4M}\sum_{a\neq
b}\mu_{ab}\left\\{\left(\sum^{M}_{i}S^{a}_{i}S^{b}_{i}\right)^{2}-M\right\\}\right]$
$\displaystyle\simeq\underset{\\{S_{i}^{a}\\}}{\mathrm{Tr}}\exp\left[\frac{M}{4}\sum_{a\neq
b}\mu_{ab}\left(\frac{1}{M}\sum^{M}_{i}S^{a}_{i}S^{b}_{i}\right)^{2}\right],$
(49)

where we have neglected the subleading terms in the large-$M$ limit. We now
introduce the delta function $\delta(MQ_{ab}-\sum_{i}^{M}S_{i}^{a}S_{i}^{b})$
using the integral representation with the variable $\lambda_{ab}$. Then we
have from Eq. (6)

$\displaystyle\overline{Z^{n}}=$ $\displaystyle
e^{nN\beta^{2}C}\int\prod_{a<b}dq_{ab}d\mu_{ab}dQ_{ab}d\lambda_{ab}$
$\displaystyle\times$
$\displaystyle\exp\Big{[}-NM\Big{\\{}-\frac{1}{4}\beta^{2}\sum_{a\neq
b}q^{2}_{ab}+\frac{1}{2}\sum_{a\neq b}\mu_{ab}q_{ab}$
$\displaystyle-\frac{1}{4}\sum_{a\neq
b}\mu_{ab}Q^{2}_{ab}+\frac{1}{2}\sum_{a\neq
b}\lambda_{ab}Q_{ab}-\ln\tilde{L}(\underline{\lambda})\Big{\\}}\Big{]}$ (50)

where

$\displaystyle\widetilde{L}(\underline{\lambda})=\underset{\\{S^{a}\\}}{\mathrm{Tr}}\exp\left[\frac{1}{2}\sum_{a\neq
b}\lambda_{ab}S^{a}S^{b}\right].$ (51)

In the large-$M$ limit, the integral is dominated by the saddle points. In
particular, the saddle point equations obtained by varying $q_{ab}$ and
$\mu_{ab}$ are, respectively,

$\displaystyle\mu_{ab}=\beta^{2}q_{ab}$ (52)

and

$\displaystyle q_{ab}=\frac{1}{2}Q^{2}_{ab}.$ (53)

Inserting this into the above equation, we can rewrite Eq. (50) as

$\displaystyle\overline{Z^{n}}=e^{nN(\beta
J)^{2}C}\int\prod_{a<b}dQ_{ab}d\lambda_{ab}\;\exp[-NM\widetilde{G}(\underline{Q},\underline{\lambda})],$
(54)

with

$\widetilde{G}(\underline{Q},\underline{\lambda})=-\frac{1}{16}(\beta
J)^{2}\sum_{a\neq b}Q^{4}_{ab}+\frac{1}{2}\sum_{a\neq
b}\lambda_{ab}Q_{ab}-\ln\tilde{L}(\underline{\lambda}).$ (55)

The free energy in the large-$M$ limit is then given by

$\displaystyle\frac{\beta F}{NM}=-(\beta J)^{2}C_{\infty}+\lim_{n\to
0}\frac{1}{n}\widetilde{G}(\underline{Q},\underline{\lambda}),$ (56)

where

$\displaystyle C_{\infty}=\lim_{M\to\infty}\frac{C}{M}=\frac{1}{16}$ (57)

Note that we have restored $J^{2}$ which sets the variance in Eq. (2)
explicitly. This free energy is exactly the same as the one for the fully
connected $p$ spin glass model with $p=4$, which is given by the Hamiltonian

$\displaystyle H=-\sum_{1\leq x_{1}<\cdots<x_{p}\leq
N}J_{x_{1},x_{2},\cdots,x_{N}}S(x_{1})S(x_{2})\cdots S(x_{p}),$ (58)

for the Ising spin $S(x)$ at site $x$. The bonds
$J_{x_{1},x_{2},\cdots,x_{N}}$ are independent random variables satisfying the
Gaussian distribution with zero mean and variance

$\displaystyle\frac{p!\tilde{J}^{2}}{2N^{p-1}}.$ (59)

The free energy for this model is given exactly the same as Eq. (56) with
$\tilde{J}^{2}=J^{2}/4$. (The formula for this correspondence for general $p$
is $\tilde{J}^{2}=4C_{\infty}J^{2}$.)

We can readily use the known results for this model. The replica symmetric
phase with $\lambda=Q=0$ has the free energy per site as

$\frac{\beta F_{\rm RS}}{N}=-\frac{(\beta\tilde{J})^{2}}{4}-\ln 2.$ (60)

The entropy per site is then given by

$\frac{S_{\rm RS}}{N}=\ln 2-\frac{(\beta\tilde{J})^{2}}{4},$ (61)

which becomes negative for temperature
$T/\tilde{J}<T^{\infty}_{*}/\tilde{J}\equiv 1/(2\sqrt{\ln 2})$. Therefore in
the original unit

$\displaystyle T^{\infty}_{*}/J=\frac{1}{4\sqrt{\ln 2}}\simeq 0.30028.$ (62)

This is the same value as that obtained in the $M\to\infty$ limit of Eq. (24).

If we use the 1RSB form for $Q_{ab}$ and $\lambda_{ab}$ in Eq. (56), the free
energy becomes

$\displaystyle\frac{\beta F^{\infty}_{\rm 1RSB}}{N}=$
$\displaystyle-\frac{(\beta\tilde{J})^{2}}{4}[1+(m_{1}-1)Q_{1}^{p}]+\frac{1}{2}(m_{1}-1)\lambda_{1}Q_{1}$
$\displaystyle+$ $\displaystyle\frac{\lambda_{1}}{2}-\ln
2-\frac{1}{m_{1}}\ln\int Dy\;\cosh^{m_{1}}(\sqrt{\lambda_{1}}y).$ (63)

The saddle point equations are as follows:

$\lambda_{1}=\frac{(\beta\tilde{J})^{2}}{2}pQ_{1}^{p-1},$ (64)

and

$Q_{1}=\frac{\int
Dy\;\cosh^{m_{1}}(\sqrt{\lambda_{1}}y)\tanh^{2}(\sqrt{\lambda_{1}}y)}{\int
Dy\;\cosh^{m_{1}}(\sqrt{\lambda_{1}}y)}.$ (65)

There is another saddle point equation which is obtained by varying the free
energy with respect to $m_{1}$:

$\displaystyle\frac{(\beta\tilde{J})^{2}}{4}Q_{1}^{p}(p-1)+\frac{1}{m_{1}^{2}}\ln\int
Dy\;\cosh^{m_{1}}(\sqrt{\lambda_{1}}y)$
$\displaystyle~{}~{}-\frac{1}{m_{1}}\frac{\int
Dy\;\cosh^{m_{1}}(\sqrt{\lambda_{1}}y)\ln(\cosh(\sqrt{\lambda_{1}}y))}{\int
Dy\;\cosh^{m_{1}}(\sqrt{\lambda_{1}}y)}=0.$ (66)

Again, when $m_{1}=1$, $F_{\rm 1RSB}$ becomes equal to $F_{\rm RS}$. We
determine the temperature $T^{\infty}_{\rm 1RSB}$ by setting $m_{1}=1$. Using
$\int Dy\;\cosh(\sqrt{\lambda_{1}}y)=e^{\lambda_{1}/2}$, we can combine Eqs.
(65), (66) and (64) to get

$\displaystyle e^{-\lambda_{1}/2}\int
Dy\;\cosh(\sqrt{\lambda_{1}}y)\Big{[}\ln\cosh(\sqrt{\lambda_{1}}y)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{p-1}{2p}\lambda_{1}\tanh^{2}(\sqrt{\lambda_{1}}y)\Big{]}-\frac{\lambda_{1}}{2}=0.$
(67)

If we define

$\displaystyle\nu\equiv\sqrt{\lambda_{1}},$ (68)

then the above equation can be rewritten as $f_{\infty}(\nu)=0$ where

$\displaystyle f_{\infty}(\nu)\equiv$ $\displaystyle~{}e^{-\nu^{2}/2}\int
Dy\;\Big{[}\cosh(\nu y)\ln\cosh(\nu y)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{p-1}{2p}\nu^{2}\frac{\sinh^{2}(\nu
y)}{\cosh(\nu y)}\Big{]}-\frac{\nu^{2}}{2}.$ (69)

This is to be compared with the corresponding Eq. (41) for finite $M$. In Fig.
6, $f_{\infty}(\nu)$ is plotted for $p=4$. From the nonzero solution and from
the corresponding $Q_{1}$ in Eq. (65) and the relation Eq. (64), we obtain
$T^{\infty}_{\rm 1RSB}/\tilde{J}\simeq 0.61688$ or in the original unit

$\displaystyle T^{\infty}_{\rm 1RSB}/J\simeq 0.30844>T^{*}_{\infty}.$ (70)

For $f(\nu)$, the small-$\nu$ expansion yields

$\displaystyle f_{\infty}(\nu)=$
$\displaystyle\left(\frac{2-p}{4p}\right)\nu^{4}+\left(\frac{2p-3}{6p}\right)\nu^{6}$
$\displaystyle+\left(\frac{5(4-3p)}{24p}\right)\nu^{8}+O(\nu^{10}).$ (71)

We can see that for $p>2$, $f_{\infty}(\nu)$ has a negative slope near the
origin. For $p=2$, the leading order term is $\nu^{6}$ with a positive
coefficient.

Figure 6: $f_{\infty}(\nu)$ defined in Eq. (69) for $p=4$. There is a nonzero
solution $\nu\simeq 2.1163$ to the equation $f_{\infty}(\nu)=0$.

###  II.5 The FRSB Solution

Here we consider the FRSB solutions. We first write the free energy in terms
of the Parisi function $q(x)$ for $0\leq x\leq 1$. It is given by

$\displaystyle\frac{\beta F_{\rm FRSB}}{N}=-C\beta^{2}-M\ln 2-\tau\langle
q^{2}\rangle$ $\displaystyle-
w_{1}\int_{0}^{1}dx\;\left\\{xq^{3}(x)+3q(x)\int_{0}^{x}dy\;q^{2}(y)\right\\}+w_{2}\langle
q^{3}\rangle$ $\displaystyle+y_{1}\langle q^{4}\rangle+y_{2}\Big{\\{}\langle
q^{4}\rangle-2\langle q^{2}\rangle^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\int_{0}^{1}dx\int_{0}^{x}dy\;(q^{2}(x)-q^{2}(y))^{2}\Big{\\}}$
$\displaystyle-y_{3}\Big{\\{}2\langle q\rangle\langle
q^{3}\rangle+\int_{0}^{1}dx\;q^{2}(x)\int_{0}^{x}dy\;(q(x)-q(y))^{2}\Big{\\}}$
$\displaystyle-y_{5}\Big{\\{}\langle q^{2}\rangle^{2}-4\langle
q\rangle^{2}\langle q^{2}\rangle$ $\displaystyle~{}~{}~{}~{}~{}-4\langle
q\rangle\int_{0}^{1}dx\;q(x)\int_{0}^{x}dy\;(q(x)-q(y))^{2}$
$\displaystyle~{}~{}~{}~{}~{}-\int_{0}^{1}dx\int_{0}^{x}dy\int_{0}^{x}dz\;(q(x)-q(y))^{2}(q(x)-q(z))^{2}\Big{\\}}$
$\displaystyle+z_{1}\langle q^{5}\rangle,$ (72)

where

$\displaystyle\langle q^{k}\rangle=\int_{0}^{1}q^{k}(x)dx.$ (73)

and we have only kept the first quintic term. The FRSB expressions for the
rest of the quintic terms are given in Appendix E.

Because the equations for the stationarity equations of the FRSB functional
equations are so cumbersome we have relegated them to the Appendices D and E.
We can only make progress in solving these equations at the quintic level by
making simplifications. The full set of quintic terms is given in Appendix E
but in Eq. (72) we have reduced them from 9 terms to just one. We choose the
numerical value of that $z_{1}$ to equal $z_{1}^{\rm{eff}}$ in Eq. (C). A
second simplification was to set $y_{5}$ =0. When this is done the
differential equation of Eq. (153) can be solved analytically. With $y_{5}$
set to be zero we do not think that does much harm to the physics of the
problem. For example, the Goldbart-Elderfield singularity Goldbart and
Elderfield (1985) still arises.

Fortunately at quartic level, that is, if we set $z_{1}=0$, one can solve the
differential equation for $q(x)$, Eq. (153), analytically. There is no need to
set $y_{5}$ to zero when just working at quartic level. Because it is a first
order differential equation, its solution depends on one adjustable constant
$x_{0}$. The result is

$\displaystyle
q(x)=\frac{w_{1}y_{3}-2w_{2}y_{5}-\frac{2(y_{3}-2xy_{5})(y_{3}^{2}-4y_{1}y_{5})x_{0}}{\sqrt{y_{1}-xy_{3}+x^{2}y_{5}}}}{2(-y_{3}^{2}+4y_{1}y_{5})}.$
(74)

Figure 7: Plots of $q(x)$ for the FRSB solution (red) and the 1RSB solution
(black) at $M=2.25$ at $\tau=-0.001$. The FRSB state is the equilibrium state
as it has the higher free energy. These plots are for the quartic theory.
$x_{1}$ for the FRSB solution is where $q(x)$ goes to zero, $x_{1}\approx
0.24437$, while the upper breakpoint $x_{2}\approx 0.25870$. Figure 8: Plots
of $q(x)$ for the FRSB solution (red) and the 1RSB solution (black) at
$M=2.50$ for $\tau=-0.01$. This calculation has been done at quintic level,
with just one quintic coefficient, with $z_{1}=z_{1}^{\rm{eff}}$ and with
$y_{5}=0$, for both the FRSB and 1RSB solutions, in order to simplify the
numerical work in the FRSB case. At this value of $M$, the first transition is
to the 1RSB state at $\tau=0$, but below the Gardner transition temperature
$T_{G}$, (which corresponds to a value of $\tau_{G}\approx-0.0078$) there is a
transition to a state with FRSB. Below $T_{G}$, this FRSB state has a higher
free energy than the corresponding 1RSB state.

Physical requirements on the choice of $x_{0}$ are that for some interval
$0<x_{1}<x<x_{2}<1$, $q(x)$ is real, an increasing function of $x$, and
positive. $x_{1}$ is for the solutions discussed in this paper at the point
where $q(x_{1})=0$, and solving this equation gives us $x_{1}$ as a function
of $x_{0}$. The upper breakpoint, $x_{2}$, is where $q(x)$ takes the constant
value $q(x_{2})$ in the interval $1>x>x_{2}$. Its value as a function of
$x_{0}$ is determined by solving Eq. (152) at the value $x=x_{2}$. This
relates the value of $x_{2}$ to $x_{0}$. The value of $x_{0}$ itself can be
determined by setting the right-hand side of Eq. (145) to zero by choosing a
value for $x_{0}$, for any value of $x>x_{1}$. The FRSB solution for the case
$M=2.25$ at a value of $\tau=-0.001$ is shown in Fig. 7. It is contrasted with
the form of $q(x)$ for the 1RSB case at the same values of $M$ and $\tau$.

Note that there is an inverse square root singularity in $q(x)$ when
$x=x_{s}$, where $y_{1}-x_{s}y_{3}+x_{s}^{2}y_{5}=0$ but this singularity, the
Goldbart-Elderfield singularity, Goldbart and Elderfield (1985), causes no
problem so long as it occurs at a value of $x_{s}$ which is greater than
$x_{2}$ or less than $x_{1}$. In the limit $\tau\to 0$, $q(x)$ also goes to
zero ($\sim|\tau|$) so Eq. (152) fixes $x_{2}\to w_{2}/w_{1}=(M-2)$. Hence a
FRSB solution can only exist if $x_{2}<x_{s}$, which translates to $M^{**}\leq
2+\sqrt{2}/3\approx 2.47140$. The free energy difference between the FRSB and
the 1RSB state differs at order $\tau^{5}$ and we have found numerically that
the coefficient of this term goes towards zero as $M\to M^{**}$. One might
have thought that one could not ignore the quintic terms when determining
$M^{**}$ as they too give a contribution of $O(\tau^{5})$. However, in the
limit when $\tau\to 0$, both the 1RSB and the FRSB solutions have their upper
breakpoints at $w_{2}/w_{1}$ and at small $\tau$ the value of $q(x)$ on the
plateau is the same for both solutions (see Fig. 7). The form of $q(x)$ for
the two solutions only differ in the interval between $x_{2}$ and $x_{1}$ and
$x_{2}-x_{1}\sim|\tau|$ itself, so in the integrals for the free energy, Eq.
(E1), the plateau regions give the contribution of $O(|\tau|^{5})$, which is
the same for both solutions, and the region of $x$ where the solutions differ
only contributes to the higher order terms in $\tau$.

For $3>M>M^{**}$ the continuous transition is to the 1RSB state. For $M>3$,
that is for $w_{2}/w_{1}>1$, the transition is discontinuous and is to the
1RSB state. We were unable to find a solution with FRSB which had a higher
free energy than the 1RSB solution at the discontinuous transition itself.

While the quintic terms are not needed to determine the value of $M^{**}$, it
was pointed out years ago that they are needed to obtain the Gardner
transition Gross  _et al._ (1985). This is the transition which arises in the
1RSB state and it is to a state with FRSB. Provided we set $y_{5}$ to zero and
just retain one of the quintic terms $z_{1}$, MATHEMATICA can analytically
solve the first order differential equation, but its explicit form is so long
that we have not included its form in this paper. In Fig. 8 we show the
resulting FRSB solution and the 1RSB solution with the same parameters when
$M=2.50$ at a temperature below the Gardner transition temperature, so that
the FRSB state has a higher free energy than the 1RSB state. Curiously the
form of the FRSB solution is nothing like that given in Ref. Gross  _et al._
(1985). They claimed that the continuously varying feature of $q(x)$ grew from
the upper plateau. However, our solution is very similar to the FRSB solution
for $M<M^{**}$, and it seems natural to us that at low enough temperature that
solution should smoothly extend into the region $M>M^{**}$ as $M$ is
increased.

A feature of the Gardner solution is that right at the critical temperature
$T_{G}$ where the Gardner state has a free energy just equal to that of the
1RSB state, its $q(x)$ is such that its derivative $dq(x)/dx$ is infinite
right at the lower break point $x_{1}$. This is because at $T_{G}$ the
Goldbart-Elderfield singularity of the quintic order solution is just at
$x_{1}$. As the temperature is reduced below $T_{G}$, this singularity occurs
below $x_{1}$, and $dq(x)/dx$ is finite at $x_{1}$ (as in Fig. 8). For
$T>T_{G}$, the FRSB solution ceases to exist.

##  III Discussion of Fluctuation Corrections and Behavior in Finite
Dimensions

Most of this paper has been concerned with calculations at mean-field level.
Our motivation to study these was because we wished to move towards the
inclusion of fluctuations about the mean-field solutions by using RG equations
to renormalize the numerous coupling constants, ($\tau$, $w_{1}$, $w_{2}$,
$y_{1}$, $\cdots$, $y_{5}$, $z_{1}$, $\cdots$, $z_{9}$) until they lie in the
region where fluctuations have become small and mean-field theory becomes
accurate. This is the same program as followed by Höller and Read Höller and
Read (2020) for the de Almeida-Thouless (AT) transition de Almeida and
Thouless (1978). This is the transition of the Ising spin glass in a field
$h$, and in the $h-T$ phase diagram there is a line, the de Almeida-Thouless
line which separates the high-temperature paramagnetic phase replica symmetric
phase from a state with some version of replica symmetry breaking. The field
theory of our problem, Eq. (25) is identical to theirs and the reader should
consult their paper for details. However, since their paper was written new
simulations have suggested a possible extension of their approach, which we
describe. We begin by briefly summarizing some of their results and
procedures.

For the quartic coefficients below $d<8$ the coefficients $y_{1}$, $y_{2}$,
$y_{3}$, $y_{4}$ and $y_{5}$ are dominated by the “box” diagrams for
dimensions $8>d>6$ and their bare values become negligible compared to the
contribution of the box diagrams, which can be expressed in terms of the
values of $w_{1}$ and $w_{2}$. For $d>8$, a good approximation to their values
is provided by the bare values of these coefficients. The important
combination of coefficients

$\displaystyle\tilde{y}(x)=Y(x)=y_{1}-xy_{3}+x^{2}y_{5},$ (75)

at the value of $x$ corresponding to the upper break point $x_{2}$ (which in
the limit $\tau\to 0$ has the value $w_{2}/w_{1}$) plays a key role in
determining the nature of the state below the transition. When
$\tilde{y}(\rho)$ is positive, (where $\rho=w_{2}/w_{1})$), the transition is
to a state with FRSB, but if it is negative the transition is to a state with
1RSB. (This is how the value of $M^{**}$ was determined in the mean-field
calculations by setting $x=\rho=M-2$ and solving $Y(x)=0$ for $M$). Höller and
Read found from the box diagrams that

$\displaystyle\tilde{y}(\rho)=K_{d}w_{1}^{4}\rho^{2}(22-48\rho-32\rho^{2}-8\rho^{3}+\rho^{4})/(8-d),$
(76)

where $K_{d}=2/(\Gamma(d/2)(4\pi)^{d/2})$, (provided $\rho<1$). Höller and
Read studied in particular the RG flow equations in dimensions $d=6+\epsilon$,
where they could employ the Bray and Roberts Bray and Roberts (1980) RG
recursion relations. Using these recursion relation, one finds that under the
RG transforms $w_{1}$ and $w_{2}$ scale down towards zero as
$\exp[-\frac{1}{2}\epsilon l]$. As $l\to\infty$ both $w_{1}$ and $w_{2}$
approach their fixed point value, (which is $0$) but their ratio
$\rho=w_{2}/w_{1}$ approaches a constant as the RG scale parameter $l$ goes to
infinity.

It is the numerical value of this constant in the large $l$ limit which
determines whether $\tilde{y}(\rho)$ is positive or negative. The polynomial
in Eq. (76) is such that $\tilde{y}(\rho)$ is positive provided $\rho<0.8418$.
Höller and Read left the ratio $\rho$ undetermined, supposing that its value
was related in some way to the bare values of the constants.

We shall argue that its value is universal and that $\rho=\frac{1}{2}$. Then
as $\frac{1}{2}<0.8418$, the state formed will have FRSB and so is in the
universality class of the Ising spin glass in a field.

The key to understanding this is the real space RG calculation of Angelini and
Biroli Angelini and Biroli (2015). This suggested that the transition at the
AT line in high dimensions might be controlled by a zero-temperature fixed
point. They found that in a simple real-space RG approximation that in high
enough dimensions, the RG flows of $h$ and $J$, the standard deviation of the
bond distribution, which are initially close to their values on the AT line at
some non-zero temperature flowed close to their value on the AT line at zero
temperature, but then veer away up the $h$-axis at $T=0$. Then the flow is
away from the fixed point at $T=0$ and $h=h_{AT}$, where $h_{AT}$ is the value
of the field $h$ on the AT line at $T=0$. In other words the RG flow is
controlled by a zero temperature fixed point. Because their RG procedure (the
Migdal-Kadanoff approximation) works well only in low dimensions it was
uncertain whether their zero-temperature fixed point scenario in high
dimensions should be trusted. However, we believe that the recent simulation
in six dimensions in Ref. Aguilar-Janita  _et al._ (2023) strongly suggests
that it should be believed. These simulations showed that in six dimensions
that the renormalized vertices related to the “bare” couplings $w_{1}$ and
$w_{2}$ were such that their ratio was close to $1/2$. But this is the same
value (i.e. $1/2$) as was found at $T=0$ in the mean-field like Bethe lattice
calculation of the same renormalized vertices in Ref. Parisi  _et al._ (2014).
We therefore shall take it that the renormalized value of $\rho$ which should
be inserted into Eq. (76) is $1/2$. As a consequence the continuous transition
from the high-temperature phase should be to a state with FRSB, and for $d<8$
the continuous 1RSB transition should no longer occur.

The same line of argument will also apply to the AT transition of spin glasses
in a field. These have been extensively studied by simulations and the most
recent of these is that of Bharadwaj et al. Vedula  _et al._ (2023). They
found numerical evidence that the AT line might not exist below six
dimensions. The absence of the AT line below six dimensions was argued for in
Ref. Moore (2012), where it was suggested that as $d\to 6$,
$h_{AT}^{2}\sim(d-6)$, where $h_{AT}$ is the AT field at $T=0$. If this is
correct then in three dimensions there would be no phase transition to a state
with replica symmetry breaking, but there could be long length scales
according to the droplet picture of spin glasses in a field McMillan (1984);
Bray and Moore (1986); Fisher and Huse (1988) and the Imry-Ma argument Imry
and Ma (1975), especially if the field is small. That structural glasses might
behave as the Ising spin glass in a field was suggested many years ago Moore
and Yeo (2006).

RG calculations are only useful when there exist long length scales in the
physics. For $w_{2}/w_{1}>1$ the transition to the 1RSB state at mean-field
level is via a discontinuous transition at which there are no long length
scales. Then the question arises as to what the fluctuation corrections then
do. Our suspicion is that the effects of such fluctuations is to drive the
system to at least smaller values of $w_{2}/w_{1}$ and possibly into the
region where the transition might be continuous when it exists. Certainly
there is no sign of a discontinuous transition in the real space RG
calculations such as Ref. Yeo and Moore (2012b). But at present we cannot
really exclude the possibility of a discontinuous transition in physical
dimensions but we note once more that the simulations of Ref. Campellone  _et
al._ (1998) found no evidence for such a transition at $M=4$ in three
dimensions.

The chief omission of our work is therefore a stronger conclusion on the
possible existence of a discontinuous transition and its dependence on the
dimensionality $d$ of the system. The only way forward for investigating this
question, especially in high dimensions close to or above $d=6$ would seem to
be simulations on the one-dimensional proxy models. In these proxy models the
form of the long range interactions between the spins can be tuned to mimic
behavior in $d$ dimensions. Indeed for the case $p=3$, $M=2$ that has already
been done Larson  _et al._ (2010). Alas at mean-field level this model has
$w_{2}/w_{1}<1$ and so it would not be expected to have a discontinuous
transition and indeed there was no sign of such in the simulation. The case
when $p=3$ and $M=3$ has $w_{2}/w_{1}=2$ Caltagirone  _et al._ (2011) and so
might be a good model to simulate as it should have a clear discontinuous
transition. The model of the type studied in this paper, $p=4$ but with $M=4$
could also be a good model to simulate using the one-dimensional proxy model:
It has also $w_{2}/w_{1}=2$.

###### Acknowledgements.

We would like to thank Jairo de Almeida for sharing his notes dating from the
seventies on the quintic terms in the presence of FRSB for the Ising spin
glass.

##  Appendix A Expansion of the free energy to the quintic order in order
parameter

We expand Eq. (8) to $O(\mu^{5})$. We first write $L\equiv 2^{nM}L^{\prime}$,
where

$\displaystyle
L^{\prime}\equiv\mathrm{Tr}^{\prime}_{\\{S^{a}_{i}\\}}\exp\left[\frac{1}{2M}\sum_{(a,b)}\mu_{ab}f_{ab}\right].$
(77)

Here $\mathrm{Tr}^{\prime}\equiv 2^{-nM}\mathrm{Tr}$ satisfies
$\mathrm{Tr}^{\prime}_{\\{S^{a}_{i}\\}}1=1$, and we define

$\displaystyle
f_{ab}\equiv\sum_{\alpha=1}^{K}\Psi^{a}_{\alpha}\Psi^{b}_{\alpha},$ (78)

where
$\bm{\Psi}^{a}=(S^{a}_{1}S^{a}_{2},S^{a}_{1}S^{a}_{3},\cdots,S^{a}_{M-1}S^{a}_{M})$
is a $K$-dimensional vector for each replica index $a$ with components
$\Psi^{a}_{\alpha}$, $\alpha=1,2,\cdots,K\equiv M(M-1)/2$. The expansion of
$L^{\prime}$ to $O(\mu^{5})$ has the following structure:

$\displaystyle
L^{\prime}=1+\tilde{t}_{2}\sum_{(a,b)}\mu^{2}_{ab}+\tilde{w}_{1}\sum_{(a,b,c)}\mu_{ab}\mu_{bc}\mu_{ca}+\tilde{w}_{2}\sum_{(a,b)}\mu^{3}_{ab}$
$\displaystyle+\tilde{y}_{1}\sum_{a,b}\mu^{4}_{ab}+\tilde{y}_{2}\sum_{(a,b,c)}\mu^{2}_{ab}\mu^{2}_{bc}+\tilde{y}_{3}\sum_{(a,b,c)}\mu^{2}_{ab}\mu_{bc}\mu_{ca}$
$\displaystyle+\tilde{y}_{5}\sum_{(a,b,c,d)}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{da}+\tilde{d}_{1}\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu^{2}_{cd}$
$\displaystyle+\tilde{z}_{1}\sum_{(a,b)}\mu^{5}_{ab}+\tilde{z}_{2}\sum_{(a,b,c)}\mu^{3}_{ab}\mu^{2}_{bc}+\tilde{z}_{3}\sum_{(a,b,c)}\mu^{3}_{ab}\mu_{bc}\mu_{ca}$
$\displaystyle+\tilde{z}_{4}\sum_{(a,b,c)}\mu^{2}_{ab}\mu^{2}_{bc}\mu_{ca}+\tilde{z}_{5}\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{da}$
$\displaystyle+\tilde{z}_{6}\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{db}+\tilde{z}_{7}\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu^{2}_{cd}$
$\displaystyle+\tilde{z}_{8}\sum_{(a,b,c,d)}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{da}\mu_{ac}+\tilde{z}_{9}\sum_{(a,b,c,d,e)}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{de}\mu_{ea}$
$\displaystyle+\tilde{d}_{2}\sum_{(a,b,c,d)}\mu^{3}_{ab}\mu^{2}_{cd}+\tilde{d}_{3}\sum_{(a,b,c,d,e)}\mu^{2}_{ab}\mu_{cd}\mu_{de}\mu_{ec}.$
(79)

Here $(a,b)$, $(a,b,c)$, $(a,b,c,d)$ etc. indicate that the sums are over all
distinct replica indices. The coefficients are obtained by taking the trace of
the spins as we explain below.

In order to calculate the free energy, we have to take the logarithm of
$L^{\prime}$ and expand $\ln(1+x)$ to $O(\mu^{5})$. there are three
contributions to this order coming from $-(1/2)x^{2}$ part. They are

$\displaystyle-\frac{1}{2}\tilde{t}^{2}_{2}\sum_{(a,b)}\mu^{2}_{ab}\sum_{(c,d)}\mu^{2}_{cd}$
(80) $\displaystyle=$
$\displaystyle-\frac{1}{2}\tilde{t}^{2}_{2}\left[2\sum_{(a,b)}\mu^{4}_{ab}+4\sum_{(a,b,c)}\mu^{2}_{ab}\mu^{2}_{bc}+\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu^{2}_{cd}\right],$
$\displaystyle-\frac{1}{2}\cdot
2\tilde{t}_{2}\tilde{w}_{1}\sum_{(a,b)}\mu^{2}_{ab}\sum_{(c,d,e)}\mu_{cd}\mu_{de}\mu_{ec}$
(81) $\displaystyle=$
$\displaystyle-\tilde{t}_{2}\tilde{w}_{1}\Big{[}6\sum_{(a,b,c)}\mu^{3}_{ab}\mu_{bc}\mu_{ca}+6\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{db}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{(a,b,c,d,e)}\mu^{2}_{ab}\mu_{cd}\mu_{de}\mu_{ec}\Big{]},$

and

$\displaystyle-\frac{1}{2}\cdot
2\tilde{t}_{2}\tilde{w}_{2}\sum_{(a,b)}\mu^{2}_{ab}\sum_{(c,d)}\mu^{3}_{cd}$
(82) $\displaystyle=$
$\displaystyle-\tilde{t}_{2}\tilde{w}_{2}\Big{[}2\sum_{(a,b)}\mu^{5}_{ab}+4\sum_{(a,b,c)}\mu^{2}_{ab}\mu^{3}_{bc}+\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu^{3}_{cd}\Big{]}.$

Note that the last terms in Eqs. (80), (81) and (82) as well as the terms in
Eq. (79) with coefficients, $\tilde{d}_{i}$, $i=1,2,3$ have disconnected
parts. When we take the trace over the spins, we have to keep in mind that the
Ising spins must be paired to give nonvanishing contribution. For example, we
have $\mathrm{Tr}^{\prime}f_{ab}=0$ for $a\neq b$. We evaluate the first few
sets of coefficients as follows.

$\displaystyle\tilde{t}_{2}=\frac{1}{2!}\frac{2}{(2M)^{2}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}=\frac{1}{2!}\frac{1}{(2M)^{2}}2K=\frac{K}{4M^{2}},$
(83)
$\displaystyle\tilde{w}_{1}=\frac{1}{3!}\frac{8}{(2M)^{3}}\;\mathrm{Tr}^{\prime}f_{ab}f_{bc}f_{ca}=\frac{1}{3!}\frac{1}{(2M)^{3}}8K=\frac{K}{6M^{3}},$
(84)
$\displaystyle\tilde{w}_{2}=\frac{1}{3!}\frac{4}{(2M)^{3}}\;\mathrm{Tr}^{\prime}f^{3}_{ab}=\frac{1}{3!}\frac{1}{(2M)^{3}}4M(M-1)(M-2)$
$\displaystyle~{}~{}~{}=\frac{K}{6M^{3}}(M-2),$ (85)

and

$\displaystyle\tilde{d}_{1}$
$\displaystyle=\frac{1}{4!}\frac{12}{(2M)^{4}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}f^{2}_{cd}=\frac{1}{4!}\frac{1}{(2M)^{4}}12K^{2},$
(86) $\displaystyle\tilde{d}_{2}$
$\displaystyle=\frac{1}{5!}\frac{80}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f^{3}_{ab}f^{2}_{cd}$
$\displaystyle=\frac{1}{5!}\frac{1}{(2M)^{5}}80KM(M-1)(M-2),$ (87)
$\displaystyle\tilde{d}_{3}$
$\displaystyle=\frac{1}{5!}\frac{160}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f_{ab}f_{bc}f_{ca}f^{2}_{de}=\frac{1}{5!}\frac{1}{(2M)^{5}}160K^{2},$
(88)

Here all replica indices are distinct. One can see that
$\tilde{d}_{1}=\tilde{t}^{2}_{2}/2$,
$\tilde{d}_{2}=\tilde{t}_{2}\tilde{w}_{2}$ and
$\tilde{d}_{3}=\tilde{t}_{2}\tilde{w}_{1}$. Therefore all the disconnected
terms in $\ln L^{\prime}$ vanish.

We therefore have

$\displaystyle\ln L^{\prime}=$
$\displaystyle~{}\tilde{t}_{2}\sum_{(a,b)}\mu^{2}_{ab}+\tilde{w}_{1}\sum_{(a,b,c)}\mu_{ab}\mu_{bc}\mu_{ca}+\tilde{w}_{2}\sum_{(a,b)}\mu^{3}_{ab}$
$\displaystyle+$
$\displaystyle\left(\tilde{y}_{1}-\tilde{t}^{2}_{2}\right)\sum_{a,b}\mu^{4}_{ab}+\left(\tilde{y}_{2}-2\tilde{t}^{2}_{2}\right)\sum_{(a,b,c)}\mu^{2}_{ab}\mu^{2}_{bc}+\tilde{y}_{3}\sum_{(a,b,c)}\mu^{2}_{ab}\mu_{bc}\mu_{ca}+\tilde{y}_{5}\sum_{(a,b,c,d)}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{da}$
$\displaystyle+$
$\displaystyle\left(\tilde{z}_{1}-2\tilde{t}_{2}\tilde{w}_{2}\right)\sum_{(a,b)}\mu^{5}_{ab}+\left(\tilde{z}_{2}-4\tilde{t}_{2}\tilde{w}_{2}\right)\sum_{(a,b,c)}\mu^{3}_{ab}\mu^{2}_{bc}+\left(\tilde{z}_{3}-6\tilde{t}_{2}\tilde{w}_{1}\right)\sum_{(a,b,c)}\mu^{3}_{ab}\mu_{bc}\mu_{ca}+\tilde{z}_{4}\sum_{(a,b,c)}\mu^{2}_{ab}\mu^{2}_{bc}\mu_{ca}$
$\displaystyle+$
$\displaystyle\tilde{z}_{5}\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{da}+\left(\tilde{z}_{6}-6\tilde{t}_{2}\tilde{w}_{1}\right)\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{db}+\tilde{z}_{7}\sum_{(a,b,c,d)}\mu^{2}_{ab}\mu_{bc}\mu^{2}_{cd}$
$\displaystyle+$
$\displaystyle\tilde{z}_{8}\sum_{(a,b,c,d)}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{da}\mu_{ac}+\tilde{z}_{9}\sum_{(a,b,c,d,e)}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{de}\mu_{ea}.$
(89)

The first quartic coefficient is given by

$\displaystyle\tilde{y}_{1}=\frac{1}{4!}\frac{8}{(2M)^{4}}$
$\displaystyle\;\mathrm{Tr}^{\prime}f^{4}_{ab}$
$\displaystyle=\frac{1}{4!}\frac{8}{(2M)^{4}}$ $\displaystyle\Big{[}K+3K(K-1)$
$\displaystyle+3M(M-1)(M-2)(M-3)\Big{]}.$ (90)

This is valid for $M\geq 3$. For $2\leq M\leq 3$, there are not enough spins
whose combination makes the second term in the square bracket. Therefore, the
square bracket must be just $K+3K(K-1)$ for $2\leq M\leq 3$. The rest of them
are

$\displaystyle\tilde{y}_{2}=\frac{1}{4!}\frac{48}{(2M)^{4}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}f^{2}_{bc}=\frac{1}{4!}\frac{48}{(2M)^{4}}K^{2},$
(91) $\displaystyle\tilde{y}_{3}=$
$\displaystyle\frac{1}{4!}\frac{96}{(2M)^{4}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}f_{bc}f_{ca}$
$\displaystyle=$ $\displaystyle\frac{1}{4!}\frac{96}{(2M)^{4}}M(M-1)(M-2),$
(92)

and

$\displaystyle\tilde{y}_{5}=\frac{1}{4!}\frac{48}{(2M)^{4}}\;\mathrm{Tr}^{\prime}f_{ab}f_{bc}f_{cd}f_{da}=\frac{1}{4!}\frac{48}{(2M)^{4}}K.$
(93)

These are valid for $M\geq 2$.

We obtain the first quintic coefficient as

$\displaystyle\tilde{z}_{1}=\frac{1}{5!}\frac{16}{(2M)^{5}}$
$\displaystyle\;\mathrm{Tr}^{\prime}f^{5}_{ab}$ (94)
$\displaystyle=\frac{1}{5!}\frac{16}{(2M)^{5}}$
$\displaystyle\Big{[}10M(M-1)(M-2)K$ $\displaystyle+$ $\displaystyle
12M(M-1)(M-2)(M-3)(M-4)\Big{]}.$

This is valid for $M\geq 4$. For $2\leq M\leq 4$, the second term in the
square bracket should be dropped for the same reason as given for
$\tilde{y}_{1}$. The next coefficient is given for $M\geq 2$ as

$\displaystyle\tilde{z}_{2}$
$\displaystyle=\frac{1}{5!}\frac{320}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f^{3}_{ab}f^{2}_{bc}$
$\displaystyle=\frac{1}{5!}\frac{320}{(2M)^{5}}M(M-1)(M-2)K,$ (95)

The third and fourth quintic coefficients are given by

$\displaystyle\tilde{z}_{3}=\frac{1}{5!}\frac{320}{(2M)^{5}}$
$\displaystyle\;\mathrm{Tr}^{\prime}f^{3}_{ab}f_{bc}f_{ca}$
$\displaystyle=\frac{1}{5!}\frac{320}{(2M)^{5}}$
$\displaystyle\Big{[}K+3K(K-1)$ $\displaystyle+$ $\displaystyle
3M(M-1)(M-2)(M-3)\Big{]},$ (96)

and

$\displaystyle\tilde{z}_{4}=\frac{1}{5!}\frac{480}{(2M)^{5}}$
$\displaystyle\;\mathrm{Tr}^{\prime}f^{2}_{ab}f^{2}_{bc}f_{ca}$
$\displaystyle=\frac{1}{5!}\frac{480}{(2M)^{5}}$
$\displaystyle\Big{[}2M(M-1)(M-2)$ $\displaystyle+$ $\displaystyle
2M(M-1)(M-2)(M-3)\Big{]}.$ (97)

Again these expressions are valid only for $M\geq 3$. For $2\leq M\leq 3$, the
second terms in the square brackets in Eqs. (96) and (97) do not appear. The
remaining quintic coefficients are given by

$\displaystyle\tilde{z}_{5}$
$\displaystyle=\frac{1}{5!}\frac{960}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}f_{bc}f_{cd}f_{da}$
$\displaystyle=\frac{1}{5!}\frac{960}{(2M)^{5}}M(M-1)(M-2),$ (98)
$\displaystyle\tilde{z}_{6}=\frac{1}{5!}\frac{960}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}f_{bc}f_{cd}f_{db}=\frac{1}{5!}\frac{960}{(2M)^{5}}K^{2},$
(99)
$\displaystyle\tilde{z}_{7}=\frac{1}{5!}\frac{480}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f^{2}_{ab}f_{bc}f^{2}_{cd}=0,$
(100) $\displaystyle\tilde{z}_{8}$
$\displaystyle=\frac{1}{5!}\frac{960}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f_{ab}f_{bc}f_{cd}f_{da}f_{ac}$
$\displaystyle=\frac{1}{5!}\frac{960}{(2M)^{5}}M(M-1)(M-2),$ (101)

and

$\displaystyle\tilde{z}_{9}=\frac{1}{5!}\frac{384}{(2M)^{5}}\;\mathrm{Tr}^{\prime}f_{ab}f_{bc}f_{cd}f_{de}f_{ea}=\frac{1}{5!}\frac{384}{(2M)^{5}}K.$
(102)

These expressions are valid for all $M\geq 2$.

We now convert the summations over replica indices in Eq. (89) into those
without any restriction. We obtain

$\displaystyle\ln
L^{\prime}=t^{\prime}_{2}\sum_{a,b}\mu^{2}_{ab}+w^{\prime}_{1}\sum_{a,b,c}\mu_{ab}\mu_{bc}\mu_{ca}+w^{\prime}_{2}\sum_{a,b}\mu^{3}_{ab}$
$\displaystyle+y^{\prime}_{1}\sum_{a,b}\mu^{4}_{ab}+y^{\prime}_{2}\sum_{a,b,c}\mu^{2}_{ab}\mu^{2}_{bc}+y^{\prime}_{3}\sum_{a,b,c}\mu^{2}_{ab}\mu_{bc}\mu_{ca}$
$\displaystyle+y^{\prime}_{5}\sum_{a,b,c,d}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{da}+z^{\prime}_{1}\sum_{a,b}\mu^{5}_{ab}+z^{\prime}_{2}\sum_{a,b,c}\mu^{3}_{ab}\mu^{2}_{bc}$
$\displaystyle+z^{\prime}_{3}\sum_{a,b,c}\mu^{3}_{ab}\mu_{bc}\mu_{ca}+z^{\prime}_{4}\sum_{a,b,c}\mu^{2}_{ab}\mu^{2}_{bc}\mu_{ca}$
$\displaystyle+z^{\prime}_{5}\sum_{a,b,c,d}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{da}+z^{\prime}_{6}\sum_{a,b,c,d}\mu^{2}_{ab}\mu_{bc}\mu_{cd}\mu_{db}$
$\displaystyle+z^{\prime}_{7}\sum_{a,b,c,d}\mu^{2}_{ab}\mu_{bc}\mu^{2}_{cd}+z_{8}^{\prime}\sum_{a,b,c,d}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{da}\mu_{ac}$
$\displaystyle+z^{\prime}_{9}\sum_{a,b,c,d,e}\mu_{ab}\mu_{bc}\mu_{cd}\mu_{de}\mu_{ea},$
(103)

where $t^{\prime}_{2}=\tilde{t}_{2}$, $w^{\prime}_{1}=\tilde{w}_{1}$ and
$w^{\prime}_{2}=\tilde{w}_{2}$. The first two quartic coefficients are

$\displaystyle y^{\prime}_{1}$
$\displaystyle=\tilde{y}_{1}-\tilde{t}^{2}_{2}-\left(\tilde{y}_{2}-2\tilde{t}^{2}_{2}\right)+\tilde{y}_{5}$
(104) $\displaystyle=\left(\frac{K}{24M^{4}}\right)\begin{cases}2,&\text{if}\
2\leq M\leq 3\\\ (3M^{2}-15M+20),&\text{if}\ M\geq 3\end{cases}$

and

$\displaystyle
y^{\prime}_{2}=\tilde{y}_{2}-2\tilde{t}^{2}_{2}-2\tilde{y}_{5}=-\left(\frac{K}{4M^{4}}\right).$
(105)

The rest of them are the same as when the summations are restricted.

$\displaystyle
y^{\prime}_{3}=\tilde{y}_{3},~{}~{}~{}y^{\prime}_{5}=\tilde{y}_{5}.$ (106)

The quintic coefficients are given by

$\displaystyle z^{\prime}_{1}$
$\displaystyle=\tilde{z}_{1}-2\tilde{t}_{2}\tilde{w}_{2}-\left(\tilde{z}_{2}-4\tilde{t}_{2}\tilde{w}_{2}\right)+\tilde{z}_{5}+\tilde{z}_{7}$
(107)
$\displaystyle=\left(\frac{K}{10M^{5}}\right)\begin{cases}5(M-2),&\text{if}\
2\leq M\leq 4\\\ (M-2)(M^{2}-7M+17),&\text{if}\ M\geq 4\end{cases}$
$\displaystyle z^{\prime}_{2}$
$\displaystyle=\tilde{z}_{2}-4\tilde{t}_{2}\tilde{w}_{2}-2\tilde{z}_{5}-2\tilde{z}_{7}$
$\displaystyle=-\left(\frac{K}{M^{5}}\right)(M-2),$ (108) $\displaystyle
z^{\prime}_{3}$
$\displaystyle=\tilde{z}_{3}-6\tilde{t}_{2}\tilde{w}_{1}-2\left(\tilde{z}_{6}-6\tilde{t}_{2}\tilde{w}_{1}\right)+5\tilde{z}_{9}$
(109) $\displaystyle=\left(\frac{K}{6M^{5}}\right)\begin{cases}2,&\text{if}\
2\leq M\leq 3\\\ (3M^{2}-15M+20),&\text{if}\ M\geq 3\end{cases}$
$\displaystyle z^{\prime}_{4}$
$\displaystyle=\tilde{z}_{4}-\tilde{z}_{7}-\tilde{z}_{8}$ (110)
$\displaystyle=\left(\frac{K}{2M^{5}}\right)\begin{cases}0,&\text{if}\ 2\leq
M\leq 3\\\ (M-2)(M-3),&\text{if}\ M\geq 3\end{cases}$

and

$\displaystyle
z^{\prime}_{6}=\tilde{z}_{6}-6\tilde{t}_{2}\tilde{w}_{1}-5\tilde{z}_{9}=-\left(\frac{K}{2M^{5}}\right).$
(111)

The other coefficients are unchanged, namely,

$\displaystyle z^{\prime}_{i}=\tilde{z}_{i}$ (112)

for $i=5,7,8$ and $9$.

Finally, the free energy is now given by Eq. (11) with Eq. (7). One of the
saddle point equations gives $\mu_{ab}=\beta^{2}q_{ab}$. Inserting this
relation into Eq. (11), we obtain the free energy in the form given in Eq.
(25) with

$\displaystyle
w_{i}\equiv\beta^{6}w^{\prime}_{i},~{}~{}~{}~{}~{}y_{j}\equiv\beta^{8}y^{\prime}_{j},~{}~{}~{}~{}~{}z_{k}\equiv\beta^{10}z^{\prime}_{k},$
(113)

for $i=1,2$, $j=1,2,3,5$ and $k=1,2,\cdots,9$.

##  Appendix B Small-$\sigma$ behavior of $f_{M}(\sigma)$

Here we present some steps leading to the small-$\sigma$ expansion of
$f_{M}(\sigma)$ defined in Eq. (41). As mentioned in the main text, we expand
$f_{M}(\sigma)$ up to $O(\sigma^{8})$. There are numerous terms to be
evaluated. In the following, for brevity, we only list the quantities needed
for the calculation of the $O(\sigma^{6})$-coefficient. We first write

$\displaystyle\zeta(\bm{y},\mu_{1})\equiv\frac{1}{2^{M}}\underset{\\{S_{i}\\}}{\mathrm{Tr}}\;\exp\left[\sigma\bm{y}\cdot\bm{\Psi}\right]=\sum_{j=0}^{\infty}\frac{\sigma^{j}}{j!}\zeta_{j}(\bm{y}),$
(114)

where $\sigma\equiv\sqrt{\mu_{1}/M}$. We immediately see that
$\zeta_{1}(\bm{y})=0$ since $\mathrm{Tr}\;\Psi_{\alpha}=0$. Using the fact
that $\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}=0$ for $\alpha\neq\beta$, we find
that $\zeta_{2}(\bm{y})=\sum_{\alpha}^{K}y^{2}_{\alpha}$ and
$\zeta_{3}(\bm{y})=\sum^{K}_{(\alpha,\beta,\gamma)}y_{\alpha}y_{\beta}y_{\gamma}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}$.
Higher order contributions are

$\displaystyle\zeta_{4}(\bm{y})$
$\displaystyle=\sum_{\alpha}^{K}y^{4}_{\alpha}+3\sum^{K}_{\alpha\neq\beta}y^{2}_{\alpha}y^{2}_{\beta}$
$\displaystyle+\sum^{K}_{(\alpha,\beta,\gamma,\delta)}y_{\alpha}y_{\beta}y_{\gamma}y_{\delta}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}\Psi_{\delta},$
(115) $\displaystyle\zeta_{5}(\bm{y})$
$\displaystyle=10\sum^{K}_{(\alpha,\beta,\gamma)}y^{3}_{\alpha}y_{\beta}y_{\gamma}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}$
$\displaystyle+10\sum^{K}_{(\alpha,\beta,\gamma,\delta)}y^{2}_{\alpha}y_{\beta}y_{\gamma}y_{\delta}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\beta}\Psi_{\gamma}\Psi_{\delta}$
(116)
$\displaystyle+\sum^{K}_{(\alpha,\beta,\gamma,\delta,\sigma)}y_{\alpha}y_{\beta}y_{\gamma}y_{\delta}y_{\sigma}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}\Psi_{\delta}\Psi_{\sigma},$

and

$\displaystyle\zeta_{6}(\bm{y})$
$\displaystyle=\sum_{\alpha}^{K}y^{6}_{\alpha}+15\sum^{K}_{\alpha\neq\beta}y^{4}_{\alpha}y^{2}_{\beta}+15\sum^{K}_{(\alpha,\beta,\gamma)}y^{2}_{\alpha}y^{2}_{\beta}y^{2}_{\gamma}$
(117)
$\displaystyle+20\sum^{K}_{(\alpha,\beta,\gamma,\delta)}y^{3}_{\alpha}y_{\beta}y_{\gamma}y_{\delta}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}\Psi_{\delta}$
$\displaystyle+15\sum^{K}_{(\alpha,\beta,\gamma,\delta,\sigma)}y^{2}_{\alpha}y_{\beta}y_{\gamma}y_{\delta}y_{\sigma}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\beta}\Psi_{\gamma}\Psi_{\delta}\Psi_{\sigma}$
$\displaystyle+\sum^{K}_{(\alpha,\beta,\gamma,\delta,\sigma,\mu)}y_{\alpha}y_{\beta}y_{\gamma}y_{\delta}y_{\sigma}y_{\mu}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}\Psi_{\delta}\Psi_{\sigma}\Psi_{\mu}.$

Here $(\alpha,\beta,\gamma)$, etc indicate the summation is over all distinct
indices and $K\equiv{M\choose 2}$. Performing the Gaussian integrals, we have
$\int D^{K}\bm{y}\;\zeta_{j}(\bm{y})=0$ for $j$ odd, $\int
D^{K}\bm{y}\;\zeta_{2}(\bm{y})=K$, $\int
D^{K}\bm{y}\;\zeta_{4}(\bm{y})=3K+3K(K-1)$, and

$\displaystyle\int D^{K}\bm{y}\;\zeta_{6}(\bm{y})$
$\displaystyle=15K+45K(K-1)$ $\displaystyle+15K(K-1)(K-2).$ (118)

For the calculation up to $O(\sigma^{6})$, we also need the following
quantities:

$\displaystyle\int D^{K}\bm{y}\;\zeta^{2}_{2}(\bm{y})=$ $\displaystyle
3K+K(K-1),$ (119) $\displaystyle\int D^{K}\bm{y}\;\zeta^{3}_{2}(\bm{y})=$
$\displaystyle 15K+9K(K-1)$ $\displaystyle+K(K-1)(K-2),$ (120)
$\displaystyle\int D^{K}\bm{y}\;\zeta_{2}(\bm{y})\zeta_{4}(\bm{y})=$
$\displaystyle 15K+21K(K-1)$ $\displaystyle+3K(K-1)(K-2),$ (121)
$\displaystyle\int D^{K}\bm{y}\;\zeta^{2}_{3}(\bm{y})=$ $\displaystyle
6M(M-1)(M-2).$ (122)

These expressions are valid when $K\geq 2$ or $M=3,4,5,\cdots$.

Now in the second term inside the integral in Eq. (41), we can write by
symmetry

$\displaystyle\sum_{\alpha=1}^{K}\left[\frac{1}{2^{M}}\mathrm{Tr}\;\Psi_{\alpha}\exp\left[\sigma\bm{y}\cdot\bm{\Psi}\right]\right]^{2}$
$\displaystyle=$ $\displaystyle
K\left[\frac{1}{2^{M}}\mathrm{Tr}\;\Psi_{1}\exp\left[\sigma\bm{y}\cdot\bm{\Psi}\right]\right]^{2}.$
(123)

We then define

$\displaystyle\frac{1}{2^{M}}\mathrm{Tr}\;\Psi_{1}\exp\left[\sigma\bm{y}\cdot\bm{\Psi}\right]\equiv\sum_{j=1}^{\infty}\frac{\sigma^{j}}{j!}\eta_{j}(\bm{y}).$
(124)

We find that $\eta_{1}(\bm{y})=y_{1}$,

$\displaystyle\eta_{2}(\bm{y})=\sum_{(\alpha,\beta)}y_{\alpha}y_{\beta}2^{-M}\mathrm{Tr}\Psi_{1}\Psi_{\alpha}\Psi_{\beta},$
(125)

and

$\displaystyle\eta_{3}(\bm{y})$ $\displaystyle=y_{1}+3y_{1}\sum_{\alpha\neq
1}y^{2}_{\alpha}$
$\displaystyle+\sum^{K}_{(\alpha,\beta,\gamma)}y_{\alpha}y_{\beta}y_{\gamma}\;\frac{1}{2^{M}}\mathrm{Tr}\Psi_{1}\Psi_{\alpha}\Psi_{\beta}\Psi_{\gamma}.$
(126)

For the calculation up to $O(\sigma^{6})$, we need

$\displaystyle\int D^{K}\bm{y}\;\eta^{2}_{2}(\bm{y})=4(M-2),$ (127)
$\displaystyle\int D^{K}\bm{y}\;\eta_{1}(\bm{y})\eta_{3}(\bm{y})=3K,$ (128)
$\displaystyle\int D^{K}\bm{y}\;\eta^{2}_{1}(\bm{y})\zeta_{2}(\bm{y})=K+2.$
(129)

It is now a matter of Taylor expanding the functions inside the integral in
Eq. (41) and using the above results to get the expansion coefficients in
$f_{M}(\sigma)=\sum_{j=0}^{\infty}c_{2j}(M)\sigma^{2j}$. We find that
$c_{0}=c_{2}=c_{4}=0$ and the leading order term is $O(\sigma^{6})$. We obtain

$\displaystyle c_{6}(M)=-\frac{M}{24}(M-1)(M-3).$ (130)

As mentioned in the main text, it becomes negative for $M>3$. To go up to
$O(\sigma^{8})$, we need results of more Gaussian integrals similar to Eqs.
(118)-(122) and to Eqs. (127)-(129). After a rather long calculation with the
help of symbolic algebra packages in MATHEMATICA, we obtain

$\displaystyle c_{8}(M)=-\frac{M}{48}(M-1)(3M^{2}-27M+47),$ (131)

which is valid for $K\geq 3$ or $M=3,4,5,\cdots$. We note that
$c_{8}(M=3)=7/8>0$.

##  Appendix C The 1RSB equations for the quintic Landau free energy

Here we consider the 1RSB saddle point equations corresponding to the free
energy expanded up to quintic order as given in Eq. (25). Let us assume that
$q_{ab}$ takes the 1RSB form having values $q_{1}$ on $n/m_{1}$ diagonal
blocks of size $m_{1}$ and $q_{0}=0$ outside the blocks. We can then express
the cubic, quartic, and quintic terms in $q_{ab}$ in terms of $q_{1}$ and
$m_{1}$ as we have done in Eqs. (29) and (30) for the quadratic terms. We
obtain

$\displaystyle\frac{\beta F_{\rm 1RSB}}{N}=$ $\displaystyle-C\beta^{2}-M\ln
2+\tau(m_{1}-1)q_{1}^{2}-w_{1}(m_{1}-1)(m_{1}-2)q^{3}_{1}-w_{2}(m_{1}-1)q^{3}_{1}$
$\displaystyle-
y_{1}(m_{1}-1)q^{4}_{1}-y_{2}(m_{1}-1)^{2}q^{4}_{1}-y_{3}(m_{1}-1)(m_{1}-2)q^{4}_{1}-y_{5}(m_{1}-1)(m_{1}^{2}-3m_{1}+3)q^{4}_{1}$
$\displaystyle-
z_{1}(m_{1}-1)q^{5}_{1}-z_{2}(m_{1}-1)^{2}q^{5}_{1}-z_{3}(m_{1}-1)(m_{1}-2)q^{5}_{1}-z_{4}(m_{1}-1)(m_{1}-2)q^{5}_{1}$
$\displaystyle-
z_{5}(m_{1}-1)(m_{1}^{2}-3m_{1}+3)q^{5}_{1}-z_{6}(m_{1}-1)^{2}(m_{1}-2)q^{5}_{1}-z_{7}(m_{1}-1)^{3}q_{1}^{5}$
$\displaystyle-
z_{8}(m_{1}-1)(m_{1}-2)^{2}q^{5}_{1}-z_{9}(m_{1}-1)(m_{1}-2)(m_{1}^{2}-2m_{1}+2)q^{5}_{1}.$
(132)

The saddle point equations are obtained by varying the free energy with
respect to $q_{1}$ and $m_{1}$. They are given by

$\displaystyle 2\tau q_{1}=$ $\displaystyle
3\Big{[}w_{1}(m_{1}-2)+w_{2}\Big{]}q^{2}_{1}+4\Big{[}y_{1}+y_{2}(m_{1}-1)$
$\displaystyle+y_{3}(m_{1}-2)+y_{5}(m_{1}^{2}-3m_{1}+3)\Big{]}q^{3}_{1}$
$\displaystyle+$ $\displaystyle
5\Big{[}z_{1}+z_{2}(m_{1}-1)+z_{3}(m_{1}-2)+z_{4}(m_{1}-2)$
$\displaystyle+z_{5}(m_{1}^{2}-3m_{1}+3)+z_{6}(m_{1}-1)(m_{1}-2)$
$\displaystyle+z_{7}(m_{1}-1)^{2}+z_{8}(m_{1}-2)^{2}$
$\displaystyle+z_{9}(m_{1}-2)(m_{1}^{2}-2m_{1}+2)\Big{]}q^{4}_{1}$ (133)

and

$\displaystyle\tau q^{2}_{1}=$
$\displaystyle\Big{[}w_{1}(2m_{1}-3)+w_{2}\Big{]}q^{3}_{1}+\Big{[}y_{1}+2y_{2}(m_{1}-1)$
$\displaystyle+y_{3}(2m_{1}-3)+y_{5}(3m^{2}_{1}-8m_{1}+6)\Big{]}q^{4}_{1}$
$\displaystyle+$
$\displaystyle\Big{[}z_{1}+2z_{2}(m_{1}-1)+z_{3}(2m_{1}-3)+z_{4}(2m_{1}-3)$
$\displaystyle+z_{5}(3m^{2}_{1}-8m_{1}+6)+z_{6}(3m^{2}_{1}-8m_{1}+5)$
$\displaystyle+3z_{7}(m_{1}-1)^{2}+z_{8}(3m^{2}_{1}-10m_{1}+8)$
$\displaystyle+z_{9}(4m^{3}_{1}-15m^{2}_{1}+20m_{1}-10)\Big{]}q^{5}_{1}$ (134)

Combining the above equations with the condition $q_{1}\neq 0$, we have

$\displaystyle 0=$
$\displaystyle\Big{[}-m_{1}w_{1}+w_{2}\Big{]}+2\Big{[}y_{1}-y_{3}+y_{5}m_{1}(2-m_{1})\Big{]}q_{1}$
$\displaystyle+$
$\displaystyle\Big{[}3z_{1}+z_{2}(m_{1}-1)+z_{3}(m_{1}-4)+z_{4}(m_{1}-4)$
$\displaystyle+$ $\displaystyle
z_{5}(-m^{2}_{1}+m_{1}+3)+z_{6}m_{1}(1-m_{1})-z_{7}(m_{1}-1)^{2}$
$\displaystyle+$ $\displaystyle
z_{8}(4-m^{2}_{1})+z_{9}m_{1}(-3m^{2}_{1}+10m_{1}-10)\Big{]}q^{2}_{1}$ (135)

The 1RSB transition temperature is determined by setting $m_{1}=1$ in the
above equation. We obtain

$\displaystyle(w_{2}-w_{1})+2(y_{1}-y_{3}+y_{5})q_{1}$
$\displaystyle+3(z_{1}-z_{3}-z_{4}+z_{5}+z_{8}-z_{9})q^{2}_{1}=0.$ (136)

Equivalently, we have an equation without factors of $\beta$ as

$\displaystyle(w^{\prime}_{2}-w^{\prime}_{1})+2(y^{\prime}_{1}-y^{\prime}_{3}+y^{\prime}_{5})\mu_{1}$
$\displaystyle+3(z^{\prime}_{1}-z^{\prime}_{3}-z^{\prime}_{4}+z^{\prime}_{5}+z^{\prime}_{8}-z^{\prime}_{9})\mu^{2}_{1}=0.$
(137)

From Appendix A, the coefficients are given by

$\displaystyle w_{2}-w_{1}=\frac{\beta^{6}}{12M^{2}}(M-1)(M-3),$ (138)

and

$\displaystyle y_{1}-y_{3}+y_{5}$ (139) $\displaystyle=$
$\displaystyle\begin{cases}-\frac{\beta^{8}}{48M^{3}}(M-1)(12M-29),&\text{if}\
2\leq M\leq 3\\\ \frac{\beta^{8}}{48M^{3}}(M-1)(3M^{2}-27M+47),&\text{if}\
M\geq 3.\end{cases}$

In Sec. II.5, we have defined the effective quintic coefficient $z_{1}^{\rm
eff}$ as the one that appears in the above equation, which can be calculated
from the results in Appendix A as

$\displaystyle z_{1}^{\rm eff}$ $\displaystyle\equiv
z_{1}-z_{3}-z_{4}+z_{5}+z_{8}-z_{9}$ (140)
$\displaystyle=\begin{cases}\frac{\beta^{10}}{60M^{4}}(M-1)(45M-103),\\\
-\frac{\beta^{10}}{60M^{4}}(M-1)(30M^{2}-195M+283),\\\
\frac{\beta^{10}}{60M^{4}}(M-1)(3M^{3}-57M^{2}+273M-355),\end{cases}$

In the above equation, the three cases from top to bottom correspond to the
regions, $2\leq M\leq 3$, $3\leq M\leq 4$ and $M\geq 4$, respectively. This is
related to the small-$\sigma$ expansion of $f_{M}(\sigma)$ discussed in Sec.
II.3 as follows. If we multiply Eq. (136) by $-q^{3}_{1}/2$ and use
$q_{1}=\mu_{1}/\beta^{2}=M\sigma^{2}/\beta^{2}$, Eq. (136) becomes

$\displaystyle c_{6}(M)\sigma^{6}+c_{8}(M)\sigma^{8}+c_{10}(M)\sigma^{10}=0,$
(141)

where

$\displaystyle c_{6}(M)=-\frac{M^{3}}{2\beta^{6}}(w_{2}-w_{1}),$ (142)
$\displaystyle c_{8}(M)=-\frac{M^{4}}{\beta^{8}}(y_{1}-y_{3}+y_{5}),$ (143)

and

$\displaystyle
c_{10}(M)=-\frac{3M^{5}}{2\beta^{10}}(z_{1}-z_{3}-z_{4}+z_{5}+z_{8}-z_{9}).$
(144)

##  Appendix D FRSB equations for the free energy with one quintic term

Taking a functional derivative of the free energy in Eq. (72) with respect to
$q(x)$, we have

$\displaystyle 0=$ $\displaystyle\frac{\delta}{\delta q(x)}\left(\frac{\beta
F_{\rm FRSB}}{N}\right)=-2\tau
q(x)-w_{1}\left\\{3xq^{2}(x)+3\int_{0}^{x}dy\;q^{2}(y)+6q(x)\int_{x}^{1}dy\;q(y)\right\\}+3w_{2}q^{2}(x)$
$\displaystyle+4y_{1}q^{3}(x)-4y_{2}\langle q^{2}\rangle
q(x)-y_{3}\left\\{2\langle q^{3}\rangle+6\langle q\rangle q^{2}(x)+2\langle
q^{2}\rangle
q(x)+4xq^{3}(x)-6q^{2}(x)\int_{0}^{x}dy\;q(y)-2\int_{x}^{1}dyq^{3}(y)\right\\}$
$\displaystyle-y_{5}\Bigg{\\{}4\langle q^{2}\rangle q(x)-8\langle
q\rangle^{2}q(x)-8\langle q\rangle\langle
q^{2}\rangle-4\int_{0}^{1}dx^{\prime}\;q(x^{\prime})\int_{0}^{x^{\prime}}dy\;(q(x^{\prime})-q(y))^{2}$
$\displaystyle~{}~{}~{}-4\langle
q\rangle\Big{[}3xq^{2}(x)-4q(x)\int_{0}^{x}dy\;q(y)-2\int_{x}^{1}dy\;q^{2}(y)+\int_{0}^{x}dy\;q^{2}(y)+2q(x)\int_{x}^{1}dy\;q(y)\Big{]}$
$\displaystyle~{}~{}~{}-\Bigg{[}4x^{2}q^{3}(x)-12xq^{2}(x)\int_{0}^{x}dy\;q(y)-4\int_{x}^{1}dy\;yq^{3}(y)+4xq(x)\int_{0}^{x}dy\;q^{2}(y)+4q(x)\int_{x}^{1}dy\;yq^{2}(y)$
$\displaystyle~{}~{}~{}~{}~{}~{}-4\int_{0}^{x}dyq(y)\int_{0}^{x}dz\;q^{2}(z)-4\int_{x}^{1}dy\;q(y)\int_{0}^{y}dz\;q^{2}(z)-8q(x)\int_{x}^{1}dy\;q(y)\int_{0}^{y}dz\;q(z)+8q(x)\left[\int_{0}^{x}dy\;q(y)\right]^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}+8\int_{x}^{1}dy\;q^{2}(y)\int_{0}^{y}dz\;q(z)+4q(x)\int_{x}^{1}dy\;\int_{0}^{y}dz\;q^{2}(z)\Bigg{]}\Bigg{\\}}+5z_{1}q^{4}(x).$
(145)

For $0\leq x\leq 1$ where $q^{\prime}(x)\neq 0$, we can take a derivative of
the above equation and have

$\displaystyle\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left[\frac{\delta}{\delta
q(x)}\left(\frac{\beta F_{\rm FRSB}}{N}\right)\right]=0.$ (146)

This gives us

$\displaystyle 0=$ $\displaystyle-2\tau-
w_{1}\left\\{6xq(x)+6\int_{x}^{1}dy\;q(y)\right\\}+6w_{2}q(x)$
$\displaystyle+$ $\displaystyle 12y_{1}q^{2}(x)-4y_{2}\langle q^{2}\rangle-
y_{3}\Bigg{\\{}12\langle q\rangle q(x)+12xq^{2}(x)$
$\displaystyle-12q(x)\int_{0}^{x}dy\;q(y)+2\langle
q^{2}\rangle\Bigg{\\}}-y_{5}\Bigg{\\{}4\langle q^{2}\rangle-8\langle
q\rangle^{2}$ $\displaystyle-4\langle
q\rangle\Big{[}6xq(x)-4\int_{0}^{x}dy\;q(y)+2\int_{x}^{1}dy\;q(y)\Big{]}$
$\displaystyle-\Bigg{[}12x^{2}q^{2}(x)-24xq(x)\int_{0}^{x}dy\;q(y)+4x\int_{0}^{x}dy\;q^{2}(y)$
$\displaystyle~{}~{}+4\int_{x}^{1}dy\;yq^{2}(y)-8\int_{x}^{1}dy\;q(y)\int_{0}^{y}dz\;q(z)$
$\displaystyle~{}~{}+8\left[\int_{0}^{x}dyq(y)\right]^{2}+4\int_{x}^{1}dy\int_{0}^{y}dz\;q^{2}(z)\Bigg{]}\Bigg{\\}}$
$\displaystyle+$ $\displaystyle 20z_{1}q^{3}(x)$ (147)

Taking one more derivative with respect to $x$ and divide by $q^{\prime}(x)$,
we have for $x$ with $q^{\prime}(x)\neq 0$

$\displaystyle\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left[\frac{\delta}{\delta
q(x)}\left(\frac{\beta F_{\rm FRSB}}{N}\right)\right]=0.$ (148)

This is given by

$\displaystyle 0=$ $\displaystyle-6(w_{1}x-w_{2})+24Y(x)q(x)$
$\displaystyle+12Y^{\prime}(x)\int_{x}^{1}dy\;q(y)+60z_{1}q^{2}(x),$ (149)

where

$\displaystyle Y(x)\equiv y_{1}-xy_{3}+x^{2}y_{5}.$ (150)

Taking a derivative of the above equation with respect to $x$ once again, we
have

$\displaystyle\frac{d}{dx}\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left[\frac{\delta}{\delta
q(x)}\left(\frac{\beta F_{\rm FRSB}}{N}\right)\right]=0,$ (151)

This can be written as

$\displaystyle 0=$
$\displaystyle-6w_{1}+24Y(x)q^{\prime}(x)+12Y^{\prime}(x)q(x)$
$\displaystyle+24y_{5}\int_{x}^{1}dy\;q(y)+120z_{1}q^{\prime}(x)q(x).$ (152)

Eliminating $\int_{x_{0}}^{1}dy\;q(y)$ from Eqs. (149) and (152), we have

$\displaystyle q^{\prime}(x)$ $\displaystyle=$
$\displaystyle\frac{-y_{3}w_{1}+2y_{5}w_{2}+2(-y^{2}_{3}+4y_{1}y_{5})q(x)+20z_{1}y_{5}q^{2}(x)}{4Y^{\prime}(x)(Y(x)+5z_{1}q(x))}.$
(153)

##  Appendix E FRSB expressions for all quintic terms

Here we present the expressions in terms of the Parisi function $q(x)$ for the
quintic contributions to the free energy, which is denoted by $F^{(5)}_{\rm
FRSB}$. We have

$\displaystyle\frac{\beta F_{\rm FRSB}^{(5)}}{N}$ $\displaystyle=z_{1}\langle
q^{5}\rangle-z_{2}\Big{[}-\langle q^{5}\rangle+2\langle q^{3}\rangle\langle
q^{2}\rangle+\int_{0}^{1}dx\;\int_{0}^{x}dy\;(q^{3}(y)-q^{3}(x))(q^{2}(y)-q^{2}(x))\Big{]}$
$\displaystyle-z_{3}\left[2\langle q\rangle\langle
q^{4}\rangle+\int_{0}^{1}dx\;q^{3}(x)\int_{0}^{x}dy\;(q(y)-q(x))^{2}\right]-z_{4}\left[2\langle
q^{2}\rangle\langle
q^{3}\rangle+\int_{0}^{1}dx\;q(x)\int_{0}^{x}dy\;\left(q^{2}(y)-q^{2}(x)\right)^{2}\right]$
$\displaystyle-z_{5}\Big{[}-4\langle q\rangle^{2}\langle q^{3}\rangle+\langle
q^{2}\rangle\langle q^{3}\rangle-3\langle q\rangle\langle
q^{2}h\rangle-\langle q^{3}\rangle\langle
h\rangle-\int_{0}^{1}dx\;q^{2}(x)\int_{0}^{x}dy\;(q(y)-q(x))(h(y)-h(x))\Big{]},$
$\displaystyle-z_{6}\left[-2\langle q\rangle\langle q^{2}\rangle^{2}-\langle
q^{2}\rangle\langle qh\rangle\right]-z_{7}\Big{[}2\langle q^{2}\rangle\langle
q^{3}\rangle+\langle q\rangle\langle q^{4}\rangle-4\langle q\rangle\langle
q^{2}\rangle^{2}-3\langle q^{2}\rangle\langle g\rangle+\langle q^{2}g\rangle$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-\langle
q\rangle\int_{0}^{1}dx\int_{0}^{x}dy\;\left(q^{2}(y)-q^{2}(x)\right)^{2}-\int_{0}^{1}dx\int_{0}^{x}dy\;(g(y)-g(x))(q^{2}(y)-q^{2}(x))\Big{]},$
$\displaystyle-z_{8}\left[-4\langle q^{2}\rangle\langle q^{3}\rangle-4\langle
q\rangle\langle q^{2}h\rangle-\langle
qh^{2}\rangle\right]-z_{9}\Big{[}8\langle q\rangle^{3}\langle
q^{2}\rangle-4\langle q^{2}\rangle^{2}\langle q\rangle+10\langle
q\rangle^{2}\langle qh\rangle-2\langle q^{2}\rangle\langle qh\rangle$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+2\langle q\rangle\langle
q^{2}\rangle\langle h\rangle+3\langle q\rangle\langle h^{2}\rangle+\langle
h\rangle\langle qh\rangle+2\langle
q\rangle\int_{0}^{1}dx\;q(x)\int_{0}^{x}dy\;(q(y)-q(x))(h(y)-h(x))$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\int_{0}^{1}dx\;h(x)\int_{0}^{x}dy\;dy(q(y)-q(x))(h(y)-h(x))\Big{]},$
(154)

where

$\displaystyle h(x)=\int_{0}^{x}dy\;(q(y)-q(x))^{2}$ (155) $\displaystyle
g(x)=\int_{0}^{x}dy\;\left(q^{2}(y)-q^{2}(x)\right)(q(y)-q(x))$ (156)

Stationary conditions for the free energy obtained from the quintic
contributions are quite complicated. In this Appendix, we only present

$\displaystyle
0=\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left(\frac{1}{q^{\prime}(x)}\frac{d}{dx}\right)\left[\frac{\delta}{\delta
q(x)}\left(\beta F_{\rm FRSB}^{(5)}/N\right)\right].$ (157)

This is given by

$\displaystyle 0$ $\displaystyle=60z_{1}q^{2}(x)-6z_{2}\langle q^{2}\rangle-
z_{3}\left[6\int_{0}^{x}dy\;q^{2}(y)+48q(x)\int_{x}^{1}dy\;q(y)+60xq^{2}(x)\right]$
$\displaystyle-
z_{4}\left[12\int_{x}^{1}dy\;q^{2}(y)+24q(x)\int_{x}^{1}dy\;q(y)+60xq^{2}(x)\right]-z_{5}\Big{[}-24\langle
q\rangle^{2}+6\langle q^{2}\rangle-6\langle h\rangle-72\langle q\rangle xq(x)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+36\langle
q\rangle\int_{0}^{x}dy\;q(y)-6x\int_{x}^{1}dy\;q^{2}(y)+60xq(x)\int_{0}^{x}dy\;q(y)-12\left(\int_{0}^{x}dy\;q(y)\right)^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-54x^{2}q^{2}(x)+6\int_{0}^{x}dy\;h(y)-6xh(x)\Big{]}$
$\displaystyle-z_{6}\left[-6\langle q^{2}\rangle
x\right]-z_{8}\Big{[}-24\langle q^{2}\rangle-96\langle q\rangle
xq(x)+48\langle q\rangle\int_{0}^{x}dy\;q(y)+96xq(x)\int_{0}^{x}dy\;q(y)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-24\left(\int_{0}^{x}dy\;q(y)\right)^{2}-60x^{2}q^{2}(x)-12x\int_{0}^{x}dy\;q^{2}(y)\Big{]}$
$\displaystyle-z_{9}\Big{[}12x\langle q\rangle^{2}+48x\langle
q\rangle\\{xq(x)-\int_{0}^{x}dy\;q(y)\\}+12x^{2}h(x)+6x\left(\langle
h\rangle-2\int_{0}^{x}dy\;h(y)\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+48x\left(xq(x)-\int_{0}^{x}dy\;q(y)\right)^{2}+6x\langle
h\rangle+72x\langle q\rangle\int_{0}^{x}dy\;(q(x)-q(y))-12x\langle
q^{2}\rangle+60x\langle q\rangle^{2}\Big{]}.$ (158)

## References

  * Gross  _et al._ (1985) D. J. Gross, I. Kanter, and H. Sompolinsky, “Mean-field theory of the Potts glass,” Phys. Rev. Lett. 55, 304–307 (1985).
  * Gardner (1985) E. Gardner, “Spin glasses with p-spin interactions,” Nuclear Physics B 257, 747–765 (1985).
  * Kirkpatrick  _et al._ (1989) T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, “Scaling concepts for the dynamics of viscous liquids near an ideal glassy state,” Phys. Rev. A 40, 1045–1054 (1989).
  * Kirkpatrick and Thirumalai (2015) T. R. Kirkpatrick and D. Thirumalai, “Colloquium: Random first order transition theory concepts in biology and physics,” Rev. Mod. Phys. 87, 183–209 (2015).
  * Lubchenko and Wolynes (2007) Vassiliy Lubchenko and Peter G. Wolynes, “Theory of Structural Glasses and Supercooled Liquids,” Annual Review of Physical Chemistry 58, 235–266 (2007).
  * Cavagna (2009) Andrea Cavagna, “Supercooled liquids for pedestrians,” Physics Reports 476, 51 (2009).
  * Biroli and Bouchaud (2010) G. Biroli and J. P. Bouchaud, “The Random First-Order Transition Theory of Glasses: A Critical Assessment,” in  _Structural Glasses and Supercooled Liquids: Theory, Experiment,and Applications_, edited by P. G. Wolynes and V. Lubschenko (Wiley, Singapore, 2010) Chap. 2.
  * Franz and Parisi (1999) S. Franz and G. Parisi, “Critical properties of a three-dimensional p-spin model,” The European Physical Journal B - Condensed Matter and Complex Systems 8, 417 (1999).
  * Campellone  _et al._ (1998) Matteo Campellone, Barbara Coluzzi, and Giorgio Parisi, “Numerical study of a short-range $p$-spin glass model in three dimensions,” Phys. Rev. B 58, 12081 (1998).
  * Kauzmann (1948) W. Kauzmann, “The nature of the glassy state and the behavior of liquids at low temperatures,” Chem. Reviews 43, 219 (1948).
  * Yeo and Moore (2012a) Joonhyun Yeo and M. A. Moore, “Renormalization group analysis of the $M$-$p$-spin glass model with $p=3$ and $M=3$,” Phys. Rev. B 85, 100405 (2012a).
  * Yeo and Moore (2012b) Joonhyun Yeo and M. A. Moore, “Origin of the growing length scale in $M$-$p$-spin glass models,” Phys. Rev. E 86, 052501 (2012b).
  * Yeo and Moore (2020) J. Yeo and M. A. Moore, “Possible instability of one-step replica symmetry breaking in $p$-spin Ising models outside mean-field theory,” Phys. Rev. E 101, 032127 (2020).
  * Moore (2006) M. A. Moore, “Interface Free Energies in $p$-Spin Glass Models,” Phys. Rev. Lett. 96, 137202 (2006).
  * Rudnick and Nelson (1976) Joseph Rudnick and David R. Nelson, “Equations of state and renormalization-group recursion relations,” Phys. Rev. B 13, 2208–2221 (1976).
  * Höller and Read (2020) J. Höller and N. Read, “One-step replica-symmetry-breaking phase below the de Almeida–Thouless line in low-dimensional spin glasses,” Phys. Rev. E 101, 042114 (2020).
  * de Almeida and Thouless (1978) J R L de Almeida and D J Thouless, “Stability of the sherrington-kirkpatrick solution of a spin glass model,” Journal of Physics A: Mathematical and General 11, 983 (1978).
  * Aguilar-Janita  _et al._ (2023) Miguel Aguilar-Janita, Victor Martin-Mayor, Javier Moreno-Gordo, and Juan Jesus Ruiz-Lorenzo, “Second order phase transition in the six-dimensional Ising spin glass on a field,” (2023), arXiv:2306.00569 [cond-mat.dis-nn] .
  * Caltagirone  _et al._ (2011) F. Caltagirone, U. Ferrari, L. Leuzzi, G. Parisi, and T. Rizzo, “Ising $M$-$p$-spin mean-field model for the structural glass: Continuous versus discontinuous transition,” Phys. Rev. B 83, 104202 (2011).
  * Aspelmeier  _et al._ (2008) T Aspelmeier, A Billoire, E Marinari, and M A Moore, “Finite-size corrections in the Sherrington-Kirkpatrick model,” Journal of Physics A: Mathematical and Theoretical 41, 324008 (2008).
  * Goldbart and Elderfield (1985) P Goldbart and D Elderfield, “The failure of the Parisi scheme for spin glass models without reflection symmetry,” Journal of Physics C: Solid State Physics 18, L229 (1985).
  * Janiš  _et al._ (2013) V. Janiš, A. Kauch, and A. Klíč, “Free energy of mean-field spin-glass models: Evolution operator and perturbation expansion,” Phys. Rev. B 87, 054201 (2013).
  * Bray and Roberts (1980) A J Bray and S A Roberts, “Renormalisation-group approach to the spin glass transition in finite magnetic fields,” Journal of Physics C: Solid State Physics 13, 5405 (1980).
  * Angelini and Biroli (2015) Maria Chiara Angelini and Giulio Biroli, “Spin Glass in a Field: A New Zero-Temperature Fixed Point in Finite Dimensions,” Phys. Rev. Lett. 114, 095701 (2015).
  * Parisi  _et al._ (2014) G Parisi, F Ricci-Tersenghi, and T Rizzo, “Diluted mean-field spin-glass models at criticality,” Journal of Statistical Mechanics: Theory and Experiment 2014, P04013 (2014).
  * Vedula  _et al._ (2023) Bharadwaj Vedula, M. A. Moore, and Auditya Sharma, “Study of the de Almeida–Thouless line in the one-dimensional diluted power-law $XY$ spin glass,” Phys. Rev. E 108, 014116 (2023).
  * Moore (2012) M. A. Moore, “$1/m$ expansion in spin glasses and the de Almeida-Thouless line,” Phys. Rev. E 86, 031114 (2012).
  * McMillan (1984) W. L. McMillan, “Domain-wall renormalization-group study of the three-dimensional random Ising model,” Phys. Rev. B 30, 476–477 (1984).
  * Bray and Moore (1986) A. J. Bray and M. A. Moore, “Scaling theory of the ordered phase of spin glasses,” in  _Heidelberg Colloquium on Glassy Dynamics and Optimization_, edited by L. Van Hemmen and I. Morgenstern (Springer, New York, 1986) p. 121.
  * Fisher and Huse (1988) Daniel S. Fisher and David A. Huse, “Equilibrium behavior of the spin-glass ordered phase,” Phys. Rev. B 38, 386–411 (1988).
  * Imry and Ma (1975) Yoseph Imry and Shang-keng Ma, “Random-Field Instability of the Ordered State of Continuous Symmetry,” Phys. Rev. Lett. 35, 1399–1401 (1975).
  * Moore and Yeo (2006) M. A. Moore and J. Yeo, “Thermodynamic Glass Transition in Finite Dimensions,” Phys. Rev. Lett. 96, 095701 (2006).
  * Larson  _et al._ (2010) Derek Larson, Helmut G. Katzgraber, M. A. Moore, and A. P. Young, “Numerical studies of a one-dimensional three-spin spin-glass model with long-range interactions,” Phys. Rev. B 81, 064415 (2010).