Exact finite-size scaling spectra in the quenched and annealed Sherrington-
Kirkpatrick spin glass

Ding Wang1, 2 and Lei-Han Tang1, 3$\star$

1 Department of Physics, Hong Kong Baptist University, Kowloon Tong, Hong Kong
SAR, China

2 Department of Materials Science and Engineering, Southern University of
Science and Technology, Shenzhen, China

3 Center for Interdisciplinary Studies, Westlake University, Hangzhou 310024,
Zhejiang, China

$\star$<EMAIL_ADDRESS>
## Abstract

Fine resolution of the discrete eigenvalues at the spectral edge of an
$N\times N$ random matrix is required in many applications. Starting from a
finite-size scaling ansatz for the Stieltjes transform of the maximum
likelihood spectrum, we show that the scaling function satisfies a first-order
ODE of the Riccati type. Further transformation yields a linear second-order
ODE for the characteristic function, whose nodes yield the exact locations of
leading eigenvalues. With this new technique, we examine in detail the
spectral crossover of the annealed Sherrington-Kirkpatrick (SK) spin glass
model where a gap develops below a certain critical temperature. Our analysis
yields analytic predictions for the finite-size scaling of the spin
condensation phenomenon in both quenched and annealed SK models, which are
then compared and validated in Monte Carlo simulations of the annealed model.
More generally, rescaling of the spectral axis, adjusted to the distance of
neighboring eigenvalues, offers a powerful strategy to deal with singularities
that arise in the infinite size limit.

Copyright attribution to authors.
This work is a submission to SciPost Physics.
License information to appear upon publication.
Publication information to appear upon publication. Received Date
Accepted Date
Published Date

###### Contents

  1. 1 Introduction
  2. 2 The SK model and its spectral representation
    1. 2.1 Quenched, annealed, Mattis models and the Nishimori line
    2. 2.2 Random matrix formulation and spectral representation
    3. 2.3 Condensate formation at low temperatures
  3. 3 The maximum likelihood spectrum
    1. 3.1 Stieltjes transform and solution at infinite size
    2. 3.2 Solution at large but finite size
      1. 3.2.1 Scaling analysis at the spectral edge
      2. 3.2.2 Largest eigenvalues in the quenched model
      3. 3.2.3 Spectral crossover in the annealed model
  4. 4 Finite-size scaling of collective spin components
    1. 4.1 Maximum likelihood amplitudes
      1. 4.1.1 Spin condensation in the quenched model
      2. 4.1.2 Spin condensation in the annealed model
    2. 4.2 Monte Carlo simulations of the annealed SK model
  5. 5 Summary and conclusions

##  1 Introduction

A recent paper by Foini and Kurchan[1] rekindled interest in the annealed
spin-glass and related models where interactions between the spins, for
example, are allowed to evolve according to the rules of equilibrium
thermodynamics. At sufficiently low temperatures, spins and interactions
conspire to yield self-planted states whose spin fluctuations in each instance
resemble those of the Mattis model, a ferromagnet in disguise[2, 3, 4].
However, the translational or permutation symmetry of the annealed model
implies an exponential number of such states, with no discernible order from
spin configurations alone. As the interactions drift over time, each self-
planted spin state with “hidden Mattis order” evolves accordingly, with memory
of its previous configurations determined by the rate of change of the
coupling matrix. As such, these annealed models may be considered as
prototypes of certain glass forming systems that lack identifiable structural
regularity yet behave dynamically as a solid or an extremely viscous fluid. In
biological systems, interactions among proteins and other cellular components
can be reset via chemical modifications such as phosphorylation and
acetylation, etc [5, 6]. These could then change the way proteins bind each
other to form condensates or other spatial structures through quasi-
equilibrium processes.

The mean-field Sherrington-Kirkpatrick (SK) model offers a concrete example
where the above-mentioned phenomenon can be examined quantitatively in the
random matrix formulation. Here, onset of the Mattis-type order is identified
with gap opening at the edge of the coupling matrix spectrum. The associated
principal eigenvector then specifies the instantaneous hidden order in the
spin configurations. By integrating out the spin degrees of freedom, spectral
properties of the coupling matrix can be analyzed in a one-dimensional Coulomb
gas formulation. Interestingly, the annealed or Coulomb gas formulation can
also be used to learn about various statistical properties of a class of
matrix models, where the distribution of the elements are correlated with each
other in some way[1]. By linking the annealed problem to the theory of random
matrices, many powerful tools developed in the latter community can be
employed to generate precise mathematical results.

The aim of this paper is to extend the pedagogical discussions by Foini and
Kurchan to systems of finite size, particularly in the neighborhood of the
critical point for gap opening. This turns out to be less trivial than we
initially thought and requires an extended analysis of the Stieltjes transform
of the spectral density which, to our knowledge, has not been reported
previously in the literature[7]. By introducing a scaling ansatz, we arrive at
a first order ODE of the Riccati type which can be turned into a second order
linear ODE. In the case of the original random matrix, its solution is given
by the Airy function whose nodes yield position of the successive eigenvalues
near the spectral edge. In general, the edge spectrum of the matrix can be
determined through a numerically exact procedure. These results provide a
complete description of the onset of the hidden Mattis order in the SK model
in a large but finite system.

Besides treating finite-size effects in the annealed SK model, our approach
may also shed light on the discussion of finite-size effects in more general
matrix models. Random Matrix Theory (RMT) has emerged as a foundational tool
in understanding the statistical properties of complex systems. Initially
developed to model the energy levels of heavy nuclei [8], RMT has since
broadened its scope significantly. It now finds diverse applications across
quantum physics [9, 10], biology [11, 12], and finance [13, 14]. A common
aspect among these applications is the use of finite-sized matrices, where
size effects are critical; they can substantially influence theoretical
predictions and practical interpretations.

The paper is organized as follows. In Section 2, we introduce three variants
of the SK model and present a heuristic discussion of their thermodynamic
relations in terms of a generic phase diagram for spin glasses. Mapping to the
spectral representation of the quenched and annealed models is briefly
reviewed, with their respective Coulomb gas energy functions. Section 3
explores the maximum likelihood solution of the spectral distribution for
systems of finite size, with particular emphasis on the largest eigenvalues at
the edge of the spectrum. Their precise location is resolved through a
perturbative procedure in $1/N$. In Section 4, we present finite-size scaling
expressions for the principal and sub-dominant spin fluctuation amplitudes for
both quenched annealed models across the transition temperature. These scaling
relations are verified through Monte Carlo simulations of the annealed SK
model. Noticeable discrepancy in scaling amplitudes against the maximum
likelihood estimations is observed in the critical region.

##  2 The SK model and its spectral representation

###  2.1 Quenched, annealed, Mattis models and the Nishimori line

The Sherrington-Kirkpatrick spin-glass model [15] is originally introduced for
$N$ Ising spins $S_{i}=\pm 1$ with pair-wise interactions given by the
Hamiltonian,

$H(\\{S_{i}\\},\\{J_{ij}\\})=-\sum_{i<j}J_{ij}S_{i}S_{j},$ (1)

where the coupling constants $J_{ij}=J_{ji}=\pm 1/\sqrt{N}$ with equal
probability. The free energy of the quenched model at temperature $T=1/\beta$
is given by $F_{\rm q}=-T\ln Z_{\rm q}$ where

$Z_{\rm q}(\\{J_{ij}\\})=\sum_{\\{S_{i}\\}}\exp[-\beta
H(\\{S_{i}\\},\\{J_{ij}\\}]$ (2)

is the partition function for a given configuration of the $J_{ij}$’s. As is
well known, the model has a phase transition at $T_{\rm q}=1$, below which the
spin-glass order develops.

Compared to the quenched model, the "annealed partition function"

$Z_{\rm a}=\sum_{\\{S_{i}\\}}\bigl{\langle}\exp[-\beta
H(\\{S_{i}\\},\\{J_{ij}\\}]\bigr{\rangle}_{\\{J_{ij}\\}}=2^{N}[\cosh(\beta
N^{-1/2})]^{N(N-1)/2}$ (3)

is much easier to calculate. The corresponding free energy free energy $F_{\rm
a}=-T\ln Z_{\rm a}$, on the other hand, is a non-monotonic function of $T$ and
gives rise to a negative entropy $S_{\rm a}=-\partial F_{\rm a}/\partial T$ at
sufficiently low temperatures. On the other hand, if we consider the
$J_{ij}$’s as another set of $N_{b}=N(N-1)/2$ state variables, the “entropy
crisis” disappears, with the proper partition function given by
$\tilde{Z}_{\rm a}=2^{N(N-1)/2}Z_{\rm a}$. In fact, the ground state is now
$2^{N}$ fold degenerate with a zero temperature entropy $\tilde{S}_{\rm
a}(T=0)=N\ln 2$.

In 1976, Mattis[2] considered the case $J_{ij}={1\over
N}\sigma_{i}\sigma_{j}$, i.e., the coupling constants are not chosen
independently, but instead constructed from a “hidden” spin configuration
$\\{\sigma_{i}\\}$. A simple gauge transformation $S_{i}\rightarrow
S_{i}\sigma_{i}$ then brings Eq. (1) to the Hamiltonian of the ferromagnetic
Ising model, which again has a transition at $T_{\rm M}=1$. (Note the
difference in how the strength of the coupling scales with $N$.) The ground
states of both the annealed and the Mattis models are unfrustrated, but the
annealed case has extensive degeneracy. In fact, any spin configuration can be
made a ground state through a suitable choice of the $J_{ij}$’s.

The quenched SK model can be extended to an uneven distribution of
ferromagnetic and antiferromagnetic bonds. Let $p$ be the percentage of bonds
with $J_{ij}=J>0$, and $1-p$ the percentage of bonds with $J_{ij}=-J$.
Nishimori discovered a special line on the $p$-$T$ plane where configuration-
dependent quantities, such as the internal energy, acquire the same value when
averaged over the equilibrium distributions in the quenched ensemble, or
obtained from the annealed model at the same temperature[16, 17]. In the
present case, the line is given by

$e^{2J/T}={p\over 1-p}.$ (4)

When $p$ is significantly larger than $1/2$, we expect the ground state of the
quenched model to be ferromagnetic. An estimate for the stability boundary of
this phase can be obtained from the following scaling argument. The mean
energy for a single spin is given by $e_{\rm f}\simeq(2p-1)NJ$, and its
standard deviation $\delta e_{\rm f}\simeq[4p(1-p)]^{1/2}N^{1/2}J$. Setting
$e_{\rm f}=\delta e_{\rm f}$ yields a critical value $p_{c}$ at the onset of
the fluctuation-dominated spin-glass phase. In terms of the scaled variable
$r=(2p-1)N^{1/2}$, we have $r_{c}\simeq 1$. The temperature scale for thermal-
induced transition is set by $e_{\rm f}=\hat{p}_{c}N^{1/2}J$. Under our
convention $J=N^{-1/2}$, the critical temperature $T_{c}$ for the spin-glass
transition is of order 1.

Combining with Eq. (4), we show in Fig. 1 the generic phase diagram for the
extended SK model, where the spin glass transition expands to a line from P to
N[18, 17]. As temperature decreases, the annealed SK model transits from the
paramagnetic phase to the “ferromagnetic” or Mattis phase along the Nishimori
line, passing through the multicritical point N. The nature of this transition
is our main focus here.

Figure 1: Schematic phase diagram of the extended SK model. The spin glass
phase has a finite width in the scaled variable $r$. The Nishimori point N is
a multicritical point signifying the onset of Mattis order in the annealed
model.

###  2.2 Random matrix formulation and spectral representation

Given the quadratic form of Eq. (1), it is tempting to rewrite the Hamiltonian
in normal coordinates $s_{k}=\vec{S}\cdot\vec{V}_{k}$,

$H=-\frac{1}{2}\vec{S}\mathbb{J}\vec{S}^{T}=-\frac{1}{2}\sum_{k}^{N}{\lambda_{k}s_{k}^{2}}$
(5)

Here $\vec{S}=(S_{1},S_{2},...,S_{N})$ represents the spin state, and
$\vec{V}_{k},k=1,\ldots,N$, are the normalized eigenvectors of the $N\times N$
dimensional matrix $\mathbb{J}=\\{J_{ij}\\}$, with corresponding eigenvalues
$\lambda_{1}<\lambda_{2}<\ldots<\lambda_{N}$ arranged in ascending order. For
the sake of convenience, we refer to $\vec{V}_{N}$, the eigenvector associated
with the largest eigenvalue $\lambda_{N}$, as the principal eigenvector, and
$s_{N}$ the principal eigenmode.

Under the spectral decomposition (5), summation over the spins in Eqs. (2) and
(3) can be readily carried out, leaving behind an effective energy function
governing the statistics of the eigenvalues. Computation of the phase diagram
shown in Fig. 1 is then reduced to the analysis of the eigenvalue spectrum
that enables various universal features and exact results to be derived thanks
to the extensive literature on the random matrix theory. For completeness, we
sketch the main steps of this approach below and refer the reader to Ref. [1]
for additional technical details.

The real symmetric matrix $\mathbb{J}=\\{J_{ij}\\}$ of independently
distributed coupling constants forms the well-known Gaussian-orthogonal
ensemble (GOE). Their eigenvalues are distributed according to

$P_{\rm q}(\vec{\lambda})\sim\exp[-Nf_{\rm q}(\vec{\lambda})],$ (6)

where $\vec{\lambda}=(\lambda_{1},\lambda_{2},\ldots,\lambda_{N})$ and

$f_{\rm q}(\vec{\lambda})={1\over 2}\sum_{k=1}^{N}V_{0}(\lambda_{k})-{1\over
N}\sum_{1\leq j<k\leq N}\ln|\lambda_{k}-\lambda_{j}|.$ (7)

Here the potential $V_{0}(\lambda)$ is related to the bare distribution of the
$J_{ij}$’s, e.g., whether it is continuous or taking only discrete values.
However, when the central limit theorem applies, we may write
$V_{0}(\lambda)=\lambda^{2}/(2N\sigma^{2})$ [19], where $\sigma$ is the
standard deviation of the $J_{ij}$’s which is equal to $1/\sqrt{N}$ in our
case.

In the annealed SK model, collective spin fluctuations modify the distribution
of the eigenvalues. Previous studies [1] and our own numerical simulations
discussed in Sec.4 show that the discreteness of Ising spins does not play an
important role in determining $P_{\rm a}(\vec{\lambda})$ when $N$ is
sufficiently large. Moving to the spherical model of continuous spins under
the constraint $\sum_{k=1}^{N}s_{k}^{2}=N$, we may write,

$P_{\rm a}(\vec{\lambda})=A\int[d\vec{s}]P_{\rm
q}(\vec{\lambda})e^{{\beta\over
2}\sum_{k}\lambda_{k}s_{k}^{2}}\delta\Bigl{(}\sum_{k}s_{k}^{2}-N\Bigr{)},$ (8)

where $A$ is a normalization constant. Integration over the spin components
can be readily carried out with the help of the identity

$\delta(x)={1\over 2\pi i}\int_{z_{0}-i\infty}^{z_{0}+i\infty}e^{zx}dz,$ (9)

where $z_{0}$ is an arbitrary real number. We further make the substitution
$z\rightarrow\beta z/2$ and obtain,

$P_{\rm a}(\vec{\lambda})={\rm const}\times\int dz\exp[-Nf_{\rm
a}(\vec{\lambda},z)],$ (10)

where

$f_{\rm a}(\vec{\lambda},z)=f_{\rm q}(\vec{\lambda})+{1\over
2N}\sum_{k=1}^{N}\ln(z-\lambda_{k})-{\beta\over 2}z.$ (11)

The parameter $z$ in Eq. (11) needs to be greater than the largest eigenvalue
$\lambda_{N}$.

A physical interpretation of the expression (11) was given, i.e. the energy of
a linear array of $N$ negative charges located at the $\lambda_{k}$’s and a
positive charge at $z$, interacting with a logarithmic potential. Furthermore,
the negative charges experience a harmonic confining potential while the
positive charge, confined to the right of all the negative charges,
experiences a linear potential whose strength increases with $\beta$.

Figure 2: One-dimensional Coulomb gas of $N$ identical charges in a confining
potential $V(\lambda)$ and an extra half-charge of opposite sign at the right
end.

###  2.3 Condensate formation at low temperatures

For a given coupling matrix $\mathbb{J}$ and a set of eigenvalues
$\\{\lambda_{k}\\}$, equipartition theorem under the Hamiltonian (5) yields
the thermal average

$\langle s_{k}^{2}\rangle={T\over z-\lambda_{k}},$ (12)

where $z$ (the saddle point solution) is chosen to satisfy the spherical
constraint

$\sum_{k=1}^{N}\langle s_{k}^{2}\rangle=\sum_{k=1}^{N}{T\over
z-\lambda_{k}}=N.$ (13)

As the temperature decreases, $z$ moves closer to the spectral edge
$\lambda_{N}$. Below a critical temperature $T_{c}$, $\langle
s_{N}^{2}\rangle$ grows to a magnitude of order $N$, i.e., condensation of
spin configurations onto the principal eigenvector takes place.

Previous work has established that $T_{c}=1$ for both quenched[18] and
annealed[3, 4, 1] SK models. In the low temperature phase, the amplitude of
spin fluctuations onto the principal eigenvector is predicted to be,

$\langle s_{N}^{2}\rangle=\begin{cases}N(1-T),&\text{quenched};\\\
N(1-T^{2}),&\text{annealed}.\end{cases}\qquad\qquad\text{(low temperature
phase)}$ (14)

Thus $z-\lambda_{N}$ needs to be of order $N^{-1}$ so as to be consistent with
Eq. (12). This is much smaller than the typical size of the gap
$\lambda_{N}-\lambda_{N-1}$ which is of order $N^{-2/3}$[9] or greater. To
resolve the condensation process in the critical region, one needs to examine
carefully granularity of the spectrum at the spectral edge.

##  3 The maximum likelihood spectrum

###  3.1 Stieltjes transform and solution at infinite size

In the limit of large $N$, the probability distributions Eqs. (6) and (10) are
sharply peaked at their respective maximum likelihood eigenvalues that satisfy

${\partial f_{\rm q,a}(\vec{\lambda})\over\partial\lambda_{k}}={1\over
2}V^{\prime}(\lambda_{k})-{1\over N}\sum_{j\neq
k}{1\over\lambda_{k}-\lambda_{j}}=0.\qquad(k=1,\ldots,N)$ (15)

Here $V(\lambda)=V_{0}(\lambda)$ in the quenched case and
$V(\lambda)=V_{0}(\lambda)+N^{-1}\ln(z-\lambda)$ in the annealed case.

To explore analytic properties of the maximum likelihood spectrum, it is
customary to introduce Stieltjes transform

$g_{N}(x)={1\over N}\sum_{k=1}^{N}{1\over x-\lambda_{k}}$ (16)

where the $\lambda_{k}$’s are simple poles of $g_{N}(x)$ on the complex plane.
Under Eq. (15), $g_{N}(x)$ obeys the nonlinear first order ordinary
differential equation[19],

$V^{\prime}(x)g_{N}(x)-\Pi_{N}(x)=g_{N}^{2}(x)+{1\over N}g_{N}^{\prime}(x),$
(17)

with

$\Pi_{N}(x)\equiv{1\over
N}\sum_{j=1}^{N}{V^{\prime}(x)-V^{\prime}(\lambda_{j})\over x-\lambda_{j}}.$
(18)

For $V_{0}(\lambda)=\lambda^{2}/2$, simple algebra gives

$\Pi_{N}(x)=\begin{cases}1,&\text{quenched};\\\
1+N^{-1}g_{N}(z)/(x-z),&\text{annealed}.\end{cases}$ (19)

The parameter $z$ satisfies Eq. (13) which now takes the form

$g_{N}(z)=\beta.$ (20)

In the limit $N\rightarrow\infty$, Eq. (17) reduces to a quadratic equation
for $g_{\infty}(x)$ whose solution is given by

$g_{\infty}(x)={x-\sqrt{x^{2}-4}\over 2}.$ (21)

As required by Eq. (16), $g_{\infty}(x)$ decreases monotonically with $x$ on
the $x>2$ side, starting from $g_{\infty}(2)=1$. The corresponding spectral
density of eigenvalues satisfies the well-known semi-circle law[19],

$\displaystyle\rho(\lambda)=\frac{\sqrt{4-\lambda^{2}}}{2\pi}.$ (22)

As discussed in detail in Ref. [1], Eqs. (20) and (21) are compatible with
each other only when $\beta<1$ or $T>1$. In this regime, the two equations
combine to give

$z(\beta)=\beta+\beta^{-1}.$ (23)

For $\beta>1$, the principal eigenvalue $\lambda_{N}$ merges with $z$ such
that its contribution to $g_{N}(z)$ becomes of order 1 and needs to be
isolated. Excluding the term $k=N$ from the sum in (16), we rewrite Eq. (20)
as,

$g_{N-1}(z)+{1\over N}{1\over z-\lambda_{N}}=\beta.$ (24)

In the quenched case, replacing $g_{N-1}(z)$ with $g_{\infty}(z)$ yields

$z_{\rm q}(\beta)=\lambda_{N}+{1\over N(\beta-1)}.\qquad\qquad(\beta>1)$ (25)

In the annealed case, $z$ and $\lambda_{N}$ are obtained by solving Eq. (24)
together with Eq. (15) at $k=N$. The latter is given by,

${\lambda_{N}\over 2}-{1\over 2N}{1\over z_{\rm
a}-\lambda_{N}}-g_{N-1}(\lambda_{N})=0.$ (26)

In the large $N$ limit, we have $\lambda_{N}\simeq z_{\rm
a}(\beta)=\beta+\beta^{-1}$ with

$z_{\rm a}(\beta)-\lambda_{N}=\beta/[N(\beta^{2}-1)].\qquad\qquad(\beta>1)$
(27)

Equations (25) and (27) give the results in Eq. (14).

###  3.2 Solution at large but finite size

####  3.2.1 Scaling analysis at the spectral edge

To gain a better understanding of the spectral evolution across the transition
at $T_{c}=1$, we examine Eq. (17) more systematically, treating
$\epsilon\equiv 1/N$ as a small parameter. Specializing on
$V_{0}(\lambda)=\lambda^{2}/2$, we write Eq. (17) explicitly in the form,

$g_{N}^{2}(x)-xg_{N}(x)+1=-\epsilon
g_{N}^{\prime}(x)+\epsilon{g_{N}(x)-g_{N}(z)\over x-z}.\qquad\text{(annealed
model)}$ (28)

The quenched case is recovered by dropping the last term on the right-hand-
side of Eq. (28).

Starting from the zeroth order solution (21) and focusing on the edge of the
spectrum at $x\simeq 2$, we introduce the scaling ansatz

$g_{N}(x)={x\over
2}-\epsilon^{1/3}\varphi\bigl{(}(x-2)\epsilon^{-2/3}\bigr{)}.$ (29)

The expanded spectral axis on the scale of the gap between leading eigenvalues
allows one to resolve individual poles of $g_{N}(x)$. Substituting (29) into
Eq. (28) and keeping terms up to order $\epsilon^{2/3}$, we obtain,

$\varphi^{\prime}(u)=\varphi^{2}(u)-u+{\varphi(u)-\varphi(c)\over
u-c}.\qquad\text{(annealed model)}$ (30)

Here $u=(x-2)N^{2/3}$ and $c=(z-2)N^{2/3}$ are the scaled variables. To be
consistent with Eq. (21), we demand

$\varphi(u)\simeq u^{1/2}$ (31)

at large and positive $u$.

At a given temperature $T=1/\beta$, we introduce a scaled variable

$\Delta\equiv{1\over 2}N^{1/3}\Bigl{(}{1\over\beta}-\beta\Bigr{)}$ (32)

to measure deviations from the critical point at $\Delta=0$. Note that
$(\beta^{-1}-\beta)/2\simeq T-T_{c}$. Equation (20) with (23) yields, to the
leading order in $\epsilon$,

$\varphi(c)=\Delta.$ (33)

Apart from the nonlocal term on the right-hand-side, Eq. (30) is known as a
first order ordinary differential equation of the Riccati type. The
substitution

$\varphi(u)=-\frac{Q^{\prime}(u)}{Q(u)}$ (34)

brings it to a linear second order equation,

$Q^{\prime\prime}(u)=uQ(u)+\frac{Q^{\prime}(u)+\Delta Q(u)}{u-c},$ (35)

with the condition $Q^{\prime}(c)=-\Delta Q(c)$. The asymptotic form (31) sets
the boundary condition

$\lim_{u\rightarrow\infty}Q(u)=0.$ (36)

The transformation (34) converts poles of $\varphi(u)$ at $u_{1}>u_{2}>\ldots$
(starting from the largest one on the right) into nodes of $Q(u)$, where
$Q(u_{k})=0$. The eigenvalues at the leading edge are then obtained, to the
leading order in $\epsilon$,

$\lambda_{N-k+1}=2+u_{k}N^{-2/3}.\qquad\qquad(k=1,2,\ldots)$ (37)

####  3.2.2 Largest eigenvalues in the quenched model

For the quenched model, spectrum of the coupling matrix is not affected by
spin fluctuations. The equation governing $Q(u)$ can be obtained from (35) by
taking the limit $c\rightarrow\infty$, i.e.,

$Q^{\prime\prime}(u)=uQ(u).\qquad\qquad\text{(quenched)}$ (38)

Solutions to Eq. (38) can be expressed as a linear superposition of the Airy
functions ${\rm Ai}(u)$ and ${\rm Bi}(u)$, whose asymptotic behavior at $u\gg
1$ are given by $\exp(-2u^{3/2}/3)$ and $\exp(2u^{3/2}/3)$, respectively. The
requirement (36) selects $Q(u)={\rm Ai}(u)$.

On the $u<0$ side, Eq. (38) resembles that of a harmonic oscillator but with a
spring constant that grows with $|u|$. This observation suggests approximate
solutions of the form

$Q(u)\simeq A\sin(|u|^{3/2})+B\cos(|u|^{3/2}).$ (39)

Therefore the nodes are approximately located at

$u_{k}\simeq u_{0}-(k\pi)^{2/3},\qquad(k=1,2,\ldots)$ (40)

where $u_{0}$ is an offset. The density of nodes thus increases as
$|u_{k}|^{1/2}$, in agreement with the square-root singularity of Eq. (22) at
the spectral edge.

####  3.2.3 Spectral crossover in the annealed model

We now consider the more general case of Eq. (35), where the point $c$ needs
to be determined self-consistently at a given $\Delta$. Far away from the
origin, the last term on the right-hand-side of the equation is smaller than
the other two terms. Consequently, the asymptotic behavior of the solution is
the same as that of ${\rm Ai}(u)$.

The spurious singular behavior of Eq. (35) at $u=c$ can be dealt with by
considering a series solution in $v=u-c$,

$Q(c+v)=\sum_{n=0}^{\infty}{a_{n}v^{n}}.$ (41)

Substituting into Eq. (35), we obtain

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle-\Delta a_{0},$ (42)
$\displaystyle a_{n}$ $\displaystyle=$ $\displaystyle\frac{1}{n(n-2)}(\Delta
a_{n-1}+\Delta^{2}a_{n-2}+a_{n-3}).\qquad(n\geq 3)$ (43)

At $n=2$, terms containing $a_{2}$ cancel and one is left with

$c=\Delta^{2}.\qquad\qquad\text{(annealed model)}$ (44)

The coefficient $c_{2}$ needs to be chosen to match the asymptotic behavior
$Q(u)\sim\exp({-2u^{3/2}/3})$ at large positive $u$.

We have computed the function $Q(u)$ from the series solution. Nodes
$\\{u_{k}\\}$ can also be computed using a simple algorithm such as the
bisection method. Figure 3 shows our results for three different values of
$\Delta$. Also shown is the Airy function ${\rm Ai}(u)$ for the quenched model
(shaded curve). In general, the $k$th node $u_{k}(\Delta)$ increases
monotonically as $\Delta$ decreases, with $u_{k}(\infty)=u_{k}^{0}$ and
$u_{k}(-\infty)=u_{k-1}^{0}>u_{k}^{0}$. Here $\\{u_{k}^{0}\\}$ are nodes of
${\rm Ai}(u)$, with values of the first five nodes given in Table 1. As a
convention, we set $u_{0}^{0}=+\infty$.

$k$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$u_{k}^{0}$ | $-2.3381074$ | $-4.0879494$ | $-5.5205598$ | $-6.7867080$ | $-7.9441336$
$u_{k}^{c}$ | $-1.0187930$ | $-3.2481977$ | $-4.8200993$ | $-6.1633075$ | $-7.3721773$
Table 1: First five nodes $u_{k}^{0}$ and $u_{k}^{c}$ of ${\rm Ai}(u)$ and
$Q(u)$ at the critical point $\Delta=0$, respectively, computed to eight
significant digits using the series expansion. Figure 3: Edge of the maximum
likelihood spectrum with the largest eigenvalues $\lambda_{N-k+1}\simeq
2+u_{k}N^{-2/3}$. As the parameter $\Delta$ decreases across the critical
point $\Delta=0$, the nodes $u_{k}$ of $Q(u)$ move to the right within each of
the intervals defined by successive nodes of the Airy function ${\rm Ai}(u)$
(vertical dashed lines). The first node $u_{1}$ moves to infinity as $\Delta$
decreases further into the low temperature phase.

Figure 4 shows the first (black) and second (red) gap at the edge of the
spectrum as well as $z-\lambda_{N}=[\Delta^{2}-u_{1}(\Delta)]N^{-2/3}$ (blue)
against $\Delta$ across the critical region. On the high temperature side
($\Delta>0$), the two gaps approach their respective asymptotic values that
can be computed from Table 1). On the low temperature side,
$u_{1}(\Delta)\simeq c=\Delta^{2}$ and $u_{2}(\Delta)\rightarrow u_{1}^{0}$,
so that the first gap follows the diverging curve $\Delta^{2}-u_{1}^{0}$
indicated by the dashed line. The latter also gives the asymptotic behavior of
$(z-\lambda_{N})N^{2/3}$ on the high temperature side, in agreement with the
data shown by the blue curve. The low temperature side of the blue curve
approaches the function $-(2\Delta)^{-1}$ predicted by Eq. (27) (dash-dotted
line). Together they provide a detailed description of the spectral crossover
in the transition region for a large but finite system.

As a check of the analysis based on Stieltjes transform, we also performed
numerical solution of Eqs. (15) together with Eq. (23) for systems sizes up to
$N=512$. Our results for the scaled gaps at two selected sizes are included in
Fig. (4) for comparison. It is clear that corrections to the leading order
results are already quite small at $N=32$ and decrease further as $N$
increases.

Figure 4: Gap scaling in the critical region showing crossovers in a system of
large but finite size. On the low temperature side, the principal eigenvalue
$\lambda_{N}$ separates from the rest of the spectrum. Also shown are the
numerical solutions of Eqs. (15) together with Eq. (23) for $N=32$ (crosses)
and $N=256$ (circles).

##  4 Finite-size scaling of collective spin components

###  4.1 Maximum likelihood amplitudes

With the maximum likelihood spectrum obtained in Sec. 3.2, we are ready to
establish finite-size scaling expressions for the spin components $\langle
s_{k}^{2}\rangle$ across the transition under Eq. (12).

####  4.1.1 Spin condensation in the quenched model

In the quenched model, the eigenvalue spectrum does not change with
temperature. To predict $\langle s_{k}^{2}\rangle$ in the critical region, we
need to determine $z$ from Eq. (13). In terms of the scaled variables
$\Delta=N^{1/3}(T-T^{-1})/2$ and $c=(z-2)N^{2/3}$, Eq. (33) now takes the
form,

$\Delta=-{{\rm Ai}^{\prime}(c)\over{\rm Ai}(c)}.\qquad\qquad\text{(quenched
model)}$ (45)

Figure 5 shows $c$ against $\Delta$ (blue line) obtained from numerical
evaluation of the right-hand-side of Eq. (45). Also shown is Eq. (44) for the
annealed case (red line). The two curves approach each other on the high
temperature side. On the low temperature side, $c_{\rm q}(\Delta)$ for the
quenched model approaches the first node of ${\rm Ai}(u)$ at $u_{1}^{0}$.
Noting that this corresponds to a first order node of the Airy function, we
have

$c_{\rm q}(\Delta)\simeq u_{1}^{0}-\Delta^{-1}.\qquad\qquad(\Delta\ll-1)$ (46)

Figure 5: Relation between $c$ and $\Delta$ for quenched (blue) and annealed
(red) models. The dashed line gives the minimal value for $c$. A larger $c$
suppresses amplitudes of spin fluctuations apart from the principal component.

With $c_{\rm q}(\Delta)$ given above, we obtain spin components in the
critical region in scaling form,

$\langle s_{N-k+1}^{2}\rangle={TN^{2/3}\over c_{\rm
q}(\Delta)-u_{k}^{0}}\equiv TN^{2/3}\psi_{k}^{\rm
q}\Bigl{(}N^{1/3}{T^{2}-1\over 2T}\Bigr{)}.\qquad(k=1,2,\ldots)$ (47)

Here $\psi_{k}^{\rm q}(u)=[c_{q}(u)-u_{k}^{0}]^{-1}$ are scaling functions
associated with the leading modes.

In the low temperature phase, substituting Eq. (46) into Eq. (47), we obtain

$\begin{split}\langle s_{N}^{2}\rangle&={T\over z-\lambda_{N}}=N(1-T),\\\
\langle s_{N-k+1}^{2}\rangle&={TN^{2/3}\over
u_{1}^{0}-u_{k}^{0}}.\qquad\qquad(k=2,3,\ldots)\end{split}\qquad(T<T_{c})$
(48)

Equation (14) for the quenched case is reproduced.

Far away from the critical region on the high temperature side, $c_{\rm
q}(\Delta)\simeq\Delta^{2}\gg u_{k}^{0}$. Equation (47) then reduces to

$\langle s_{N-k+1}^{2}\rangle\simeq{T\over
z(\beta)-2}={T^{2}\over(T-1)^{2}},\qquad\qquad(T>T_{c})$ (49)

for the leading modes, where Eq. (23) is used.

####  4.1.2 Spin condensation in the annealed model

The scaling behavior of $\langle s_{k}^{2}\rangle$ in the annealed model can
be established in a similar way as in the quenched case. Quite generally,
$c=\Delta^{2}$ and hence the maximum likelihood estimation (MLE) of the
spectrum yields

$\langle s_{N-k+1}^{2}\rangle_{\rm
MLE}={TN^{2/3}\over\Delta^{2}-u_{k}(\Delta)}\equiv TN^{2/3}\psi_{k}^{\rm
a}\Bigl{(}N^{1/3}{T^{2}-1\over 2T}\Bigr{)}.\qquad(k=1,2,\ldots)$ (50)

The scaling functions $\psi_{k}^{\rm a}(u)$ differ from $\psi_{k}^{\rm q}(u)$
in two ways. One is the relation between $c$ and $\Delta$ as illustrated in
Fig. 5. The second is that $u_{k}(\Delta)$ approaches $u_{k}^{0}$ only when
$\Delta$ becomes large and positive (see Fig. 3). Therefore spin fluctuations
become nearly identical in the two cases only in the high temperature phase.

Outside the critical region and into the low temperature phase, Eqs. (23) and
(27) together give,

$\begin{split}\langle s_{N}^{2}\rangle&={T\over z-\lambda_{N}}=N(1-T^{2}),\\\
\langle s_{N-k+1}^{2}\rangle&={T\over
z-\lambda_{N-k+1}}\simeq{T^{2}\over(1-T)^{2}}.\qquad(k=2,3,\ldots)\end{split}\qquad(T<T_{c})$
(51)

Again we recover the annealed result in Eq. (14). Note that fluctuation
amplitudes in the annealed model are much weaker when compared to the quenched
case, except for the principal component.

###  4.2 Monte Carlo simulations of the annealed SK model

To gain a quantitative understanding of fluctuation effects not accounted for
by the maximum likelihood spectrum, we carried out Monte Carlo simulations of
the annealed SK model. In our implementation of the Metropolis algorithm for
both spin and $J_{ij}$ updates, a sweep of all $N$ spins (i.e., one Monte
Carlo step) is performed following an attempted update of a randomly chosen
bond. The energy involved in a single bond flip is $\Delta E_{J}=\pm
2/\sqrt{N}$, hence the acceptance rate is quite high except at very low
temperatures. The energy change $\Delta E_{S}$ for a spin flip, on the other
hand, is of order 1. In a typical simulation at a given temperature $T$, we
equilibrate an initially random configuration of spins and bonds for $5N^{2}$
Monte Carlo steps, followed by $20N^{2}$ Monte Carlo steps for data
collection. This is found to be sufficient for the computation of equilibrium
properties of systems up to $N=512$.

Figure 6 presents mean-square fluctuations of the first four components of
spin configurations obtained from our simulations. Away from the critical
region, the data agree well with the predicted behavior given by Eqs. (49) and
(51) (dashed and solid lines).

Figure 6: Mean-square amplitudes of the first four spin components against
temperature $T$ from simulations of the $N=512$ annealed SK model. Onset of
condensation is seen from the rise of the principal component below $T=1$.

To check the finite-size scaling predictions, we plot in Fig. 7 the scaled
mean-square amplitudes of the principal and second component against $\Delta$
for three different sizes. Indeed, data at different system sizes collapse
well in the whole temperature range. The dashed line in the figure is obtained
by replacing the average over the eigenvalue spectrum with the maximum
likelihood values obtained in Sec. 3. Their difference is more evident near
the transition point $\Delta=0$.

Figure 7: Finite-size scaling of the mean-square amplitude of the first
($k=N$) and second ($k=N-1$) components of spin fluctuations across the
transition. The difference between the maximum likelihood estimation (MLE) and
the full annealed SK model becomes more pronounced at the transition point
$\Delta=0$.

Figure 8 shows the ratio between the mean-square amplitudes and their MLE
values for the first four components of spin fluctuations at the critical
point. In agreement with Fig. 7, fluctuations of the principal component is
weaker than predicted based on the MLE. On the other hand, the second, third
and fourth components have stronger fluctuations in the full model as compared
to the MLE predictions. Further work is needed to understand these
observations in terms of collective fluctuations of the eigenvalue spectrum.

Figure 8: Ratio of spin fluctuation amplitudes at $T=T_{c}$ (or equivalently,
$\Delta=0$) against values predicted by the maximum likelihood estimation
(MLE). Spectral fluctuations appear to reduce the amplitude of the principal
component while enhance all others.

##  5 Summary and conclusions

This paper is a follow-up to recent observations made by Foini and Kurchan[1]
on the fascinating properties of annealed models. We introduced a perturbative
scheme to solve the exact equation satisfied by the Stieltjes transform of the
maximum likelihood spectrum of the coupling matrix at finite $N$. The discrete
set of eigenvalues at the leading edge of the spectrum can then be resolved
exactly to order $N^{-2/3}$, which appears to have been an open problem before
the reported work.

With this technique, we analyzed the spectral crossover of the annealed SK
spin glass model, which exhibits various singular behaviors at the transition
temperature $T_{c}$. Finite-size scaling forms of principal spin fluctuations
were obtained in the critical region, connecting smoothly to the known results
in the low and high temperature phases. These predictions, based on the
maximum likelihood spectrum, were then compared with Monte Carlo simulation
results. While the scaling exponents are borne out by the simulation data,
discrepancies in the scaling amplitudes were observed in the critical region.
These are presumably due to fluctuations in the eigenvalue spectrum in the
annealed model, which deserve further study.

## Acknowledgements

We would like to thank Leticia Cugliandolo and Jorge Kruchan for insightful
discussions. This work is supported by Center for Computational Science and
Engineering at Southern University of Science and Technology.

##### Author contributions

L-H T. designed the project and performed most of the analytic calculations.
D. W. carried out numerical solution of the maximum likelihood spectrum and
performed Monte Carlo simulations of the annealed SK model. Both participated
in writing up the work.

##### Funding information

The work is supported in part by the Research Grants Council of the HKSAR
under grants 12304020 and 12301723.

## References

  * [1] L. Foini and J. Kurchan,  _Annealed averages in spin and matrix models_ ,  SciPost Phys. 12, 080 (2022),  10.21468/SciPostPhys.12.3.080.
  * [2] D. C. Mattis,  _Solvable spin systems with random interactions_ ,  Physics Letters A 56, 421 (1976),  https://doi.org/10.1016/0375-9601(76)90396-0.
  * [3] Y. Kasai and A. Okiji,  _Hidden mattis phase in annealed system_ ,  Progress of Theoretical Physics 69, 20 (1983),  10.1143/PTP.69.20.
  * [4] Y. Matsuda, Y. Kasai and A. Okiji,  _Hidden mattis phase in the sherrington-kirkpatrick model_ ,  Progress of Theoretical Physics 71, 1091 (1984),  10.1143/PTP.71.1091.
  * [5] T. Hunter,  _The age of crosstalk: Phosphorylation, ubiquitination, and beyond_ ,  Molecular Cell 28, 730 (2007),  https://doi.org/10.1016/j.molcel.2007.11.019.
  * [6] C. Choudhary, C. Kumar, F. Gnad, M. L. Nielsen, M. Rehman, T. C. Walther, J. V. Olsen and M. Mann,  _Lysine acetylation targets protein complexes and co-regulates major cellular functions_ ,  Science 325, 834 (2009),  10.1126/science.1175371.
  * [7] D. Barbier, P. H. de Freitas Pimenta, L. F. Cugliandolo and D. A. Stariolo,  _Finite size effects and loss of self-averageness in the relaxational dynamics of the spherical sherrington–kirkpatrick model_ ,  Journal of Statistical Mechanics: Theory and Experiment, p. 073301 (2021),  10.1088/1742-5468/ac0900.
  * [8] E. P. Wigner,  _Characteristic vectors of bordered matrices with infinite dimensions_ ,  Annals of Mathematics 62, 548 (1955),  10.2307/1970079.
  * [9] M. L. Mehta,  _Random matrices_ , vol. 142,  Elsevier (2004).
  * [10] C. W. J. Beenakker,  _Random-matrix theory of quantum transport_ ,  Rev. Mod. Phys. 69, 731 (1997),  10.1103/RevModPhys.69.731.
  * [11] F. Luo, J. Zhong, Y. Yang and J. Zhou,  _Application of random matrix theory to microarray data for discovering functional gene modules_ ,  Phys. Rev. E 73, 031924 (2006),  10.1103/PhysRevE.73.031924.
  * [12] R. M. May,  _Will a large complex system be stable_ ,  Nature 238, 413 (1972),  10.1038/238413a0.
  * [13] M. Potters, J. P. Bouchaud and L. Laloux,  _Financial applications of random matrix theory: Old laces and new pieces_ ,  10.48550/arXiv.physics/0507111 (2005).
  * [14] V. Plerou, P. Gopikrishnan, B. Rosenow, L. A. N. Amaral, T. Guhr and H. E. Stanley,  _Random matrix approach to cross correlations in financial data_ ,  Phys. Rev. E 65, 066126 (2002),  10.1103/PhysRevE.65.066126.
  * [15] D. Sherrington and S. Kirkpatrick,  _Solvable model of a spin-glass_ ,  Phys. Rev. Lett. 35, 1792 (1975),  10.1103/PhysRevLett.35.1792.
  * [16] H. Nishimori,  _Internal energy, specific heat and correlation function of the bond-random ising model_ ,  Progress of Theoretical Physics 66, 1169 (1981),  10.1143/PTP.66.1169.
  * [17] H. Nishimori,  _Statistical physics of spin glasses and information processing: an introduction_ ,  Oxford University Press,  ISBN 9780198509417,  10.1093/acprof:oso/9780198509417.001.0001 (2001).
  * [18] J. M. Kosterlitz, D. J. Thouless and R. C. Jones,  _Spherical model of a spin-glass_ ,  Phys. Rev. Lett. 36, 1217 (1976),  10.1103/PhysRevLett.36.1217.
  * [19] M. Potters and J.-P. Bouchaud,  _A first course in random matrix theory_ ,  Cambridge University Press,  10.1017/9781108768900 (2020).