# Metastates and Replica Symmetry Breaking

C.M. Newman  Courant Institute of Mathematical Sciences, New York University,
New York, NY 10012 USA  NYU-ECNU Institute of Mathematical Sciences at NYU
Shanghai, 3663 Zhongshan Road North, Shanghai, 200062, China  N. Read
Department of Physics, Yale University, P.O. Box 208120, New Haven,
Connecticut 06520-8120, USA  Department of Applied Physics, Yale University,
P.O. Box 208284, New Haven, Connecticut 06520-8284, USA  D.L. Stein
Department of Physics and Courant Institute of Mathematical Sciences, New York
University, New York, NY 10012 USA  NYU-ECNU Institutes of Physics and
Mathematical Sciences at NYU Shanghai, 3663 Zhongshan Road North, Shanghai,
200062, China  Santa Fe Institute, 1399 Hyde Park Rd., Santa Fe, NM USA 87501

###### Abstract

In this review we define and discuss metastates, mathematical tools with
general applicability to thermodynamic systems which are particularly useful
when working with disordered or inhomogeneous short-range systems. In an
infinite such system there may be many competing thermodynamic states, which
can lead to the absence of a straightforward thermodynamic limit of local
correlation functions. A metastate is a probability measure on the infinite-
volume thermodynamic states that restores the connection between those states
and the Gibbs states observed in finite volumes. After introducing the basic
metastates and discussing their properties, we present possible scenarios for
the spin-glass phase and discuss what the metastate approach reveals about how
replica symmetry breaking would manifest itself in finite-dimensional short-
range spin glasses.

##  I Introduction

The presence of quenched disorder in condensed matter systems creates special
challenges for statistical mechanics, designed primarily to deal with thermal
disorder. A proper statistical mechanical treatment of a system whose
Hamiltonian contains quenched random variables requires averaging the free
energy, rather than the partition function, over these variables. This led to
the so-called “replica trick”, first used in the context of spin glasses by
Edwards and Anderson Edwards and Anderson (1975) (EA). Sherrington and
Kirkpatrick (SK) Sherrington and Kirkpatrick (1975) applied the replica trick
to an infinite-range spin glass model, using the “spin freezing” order
parameter proposed in Edwards and Anderson (1975), but found the solution to
be thermodynamically unstable at low temperature. The solution, found by
Giorgio Parisi in 1979 Parisi (1979), required an exotic and non-intuitive
notion of “replica symmetry breaking” (RSB); it was a stunning result that led
(in Phil Anderson’s words Anderson (1989)) to a “cornucopia” of applications,
to which this volume is dedicated.

A related problem, on which the effectiveness of theoretical statistical
mechanics depends, is the existence of a thermodynamic limit for the state in
which a system resides. When encountering the subject for the first time a
student typically learns that for a macroscopic system comprising $N\sim
10^{24}-10^{25}$ degrees of freedom with short-range interactions, the effects
of surfaces on bulk thermodynamic properties vanishes for large system size,
so that a straightforward $N\to\infty$ limit describes the thermodynamics of a
large but finite sample. This approach works well in describing the condensed
phase for systems without quenched disorder (so long as an appropriate order
parameter is identified), as well as for some systems with quenched disorder.
But when quenched disorder is present, its success cannot be guaranteed,
particularly if the low-temperature phase consists of many pure states
unrelated by any simple symmetry transformation.

This is not merely a theoretical issue. Consider for example a dilute magnetic
alloy such as CuMn, which for a range of Mn concentrations displays spin glass
behavior, and two labs which prepare their samples under similar conditions.
The locations of the magnetic impurity Mn atoms within the Cu lattice will
differ in the two samples, leading to a different distribution of
ferromagnetic and antiferromagnetic couplings; in more formal language, the
two samples have different realizations of the spin-spin couplings. It now
becomes crucial to understand what effect the difference in coupling
realizations will have on the system thermodynamics: if there is too much
dependence then universal behavior would be absent, i.e., different samples
would behave differently.

Fortunately, it can be shown that measurable global properties (i.e.,
quantities which do not depend on the configuration of any finite subset of an
infinite system) do not in general depend on coupling realization:111An
important exception is the spin or edge overlap between two states in a system
with RSB; the fact that this global quantity is “non-self-averaging” was one
of the surprises of RSB. Non-self-averaging in the context of short-range
models is discussed in Sect. IV.4.1. these global quantities include
energy/free energy per spin, magnetization per spin, transition temperature,
and so on. Thermodynamic states, however, are local: they can be considered as
enumerations of all $k$-spin correlation functions, where $k=1,2,3,\ldots$. As
such they are exquisitely sensitive to coupling realization, and the usual
straightforward approach of taking a simple thermodynamic limit can break
down. As we will see, this is especially important when dealing with systems
in which the ordering at low temperature is described by RSB.

In what follows we will discuss in further detail the breakdown in certain
cases of the straightforward thermodynamic limit (i.e., one not conditioned on
the quenched disorder), and present the concept of a metastate as a
statistical-mechanical tool designed to handle this. We will show how
metastates are useful for describing the properties of large finite volumes in
the absence of a straightforward thermodynamic limit, and investigate some of
their uses and applications. Given that the volume in which this contribution
appears is dedicated to the many applications of RSB, we will focus our
attention on how the metastate is especially useful when working with systems
in which the low-temperature ordering is or may be described by RSB, such as
EA spin glasses.

##  II Chaotic size dependence and the thermodynamic limit

We will be primarily interested in the EA Ising spin glass Edwards and
Anderson (1975) in zero field. Its (infinite-volume) Hamiltonian is given by

$H_{J}(\sigma)=-\sum_{\\{{\bf x},{\bf y}\\}}J_{{\bf{\bf xy}}}\sigma_{\bf
x}\sigma_{\bf y}$ (1)

where ${\bf x}$, ${\bf y}\ldots\in\mathbb{Z}^{d}$ are sites in the
$d$-dimensional hypercubic lattice $\mathbb{Z}^{d}$, $\sigma_{\bf x}=\pm 1$ is
an Ising spin at site ${\bf x}$ and $\\{{\bf x},{\bf y}\\}$ denotes an “edge”
in the set $\mathbb{E}^{d}$ of nearest-neighbor pairs. The couplings (or
bonds) $J_{{\bf{\bf xy}}}=J_{{\bf yx}}$ are independent, identically-
distributed continuous random variables chosen from a distribution
$\nu(dJ_{{\bf{\bf xy}}})$, with random variable $J_{{\bf xy}}$ assigned to the
edge $\\{{\bf x},{\bf y}\\}$. For simplicity we will assume that $\nu$ is
supported on the entire real line, is distributed symmetrically about zero,
and has finite variance; e.g., a Gaussian with mean zero and variance one. We
denote by $J$ a particular realization of the couplings.

Theoretically and numerically (not to mention experimentally!) one necessarily
works with finite volumes. Let $\Lambda_{L}\subset\mathbb{Z}^{d}$ denote a
cube of side $L$ centered at the origin, $|\Lambda_{L}|$ denote its volume,
i.e., the number of spins contained within $\Lambda_{L}$, and
$\partial\Lambda_{L}$ its boundary. One then considers $H_{J,L}(\sigma)$, the
Hamiltonian (1) restricted to $\Lambda_{L}$ with a specified boundary
condition on $\partial\Lambda_{L}$, usually taken to be periodic, which among
other advantages preserves the spin-flip symmetry of the Hamiltonian. The
thermodynamics within $\Lambda_{L}$ is described by the finite-volume Gibbs
distribution $\rho^{(L)}$: for any (well-behaved) function $f(\sigma)$ of some
subset of the spins within $\Lambda_{L}$, its expectation at inverse
temperature $\beta=1/T$ is given by

$\langle
f(\sigma)\rangle_{\rho^{(L)}}=\sum_{\\{\sigma_{L}\\}}f(\sigma)e^{-\beta
H_{J,L}(\sigma)}/\sum_{\\{\sigma_{L}\\}}e^{-\beta H_{J,L}(\sigma)}$ (2)

where $\sum_{\\{\sigma_{L}\\}}$ denotes a sum over all $2^{|\Lambda_{L}|}$
spin configurations within $\Lambda_{L}$. The usual procedure is to solve for
$\rho^{(L)}$ in some $\Lambda_{L}$ with $L$ large, followed by taking some
infinite sequence of volumes $\Lambda_{L}$ with $L\to\infty$. The question
then arises, under what conditions do the thermodynamic quantities under
investigation converge in the thermodynamic limit?

The question of convergence (not depending on the choice of sequence of
volumes) typically doesn’t arise in theoretical physics (outside of
mathematical physics), because it usually doesn’t present a problem. In both
the EA and SK spin glasses, most global quantities, such as those mentioned in
the Introduction, can be shown to converge Guerra and Toninelli (2002); Guerra
(2003); Talagrand (rlin) (with probability one in the coupling
realizations222Throughout this paper, all conclusions drawn about spin glasses
should be understood as occurring with probability one in the coupling
realizations unless otherwise noted.). But if we’re interested in studying the
thermodynamic states themselves (i.e., as mentioned above, the set of all
finite-spin correlation functions), then a problem may arise at low
temperatures, particularly if many pure states are present.

This was addressed in Newman and Stein (1992), which considered the EA
Hamiltonian (1) in an infinite sequence of volumes $\Lambda_{L}$. The focus
there was on “gauge-related” boundary conditions, i.e., boundary conditions
related by a gauge transformation, such as periodic and antiperiodic or any
two fixed boundary conditions (for a definition and detailed discussion of
gauge-related boundary conditions, see Sect. IX of Newman and Stein (2022)).
Several theorems were proved in Newman and Stein (1992), but both their
statements and proofs will be omitted here. For our purposes, the consequences
of the theorems can be summarized as follows.

Consider an infinite sequence of volumes with periodic boundary conditions as
described above, and suppose that there exists only a single pair of globally
spin-flip related pure states (at positive temperature) or ground states (at
zero temperature), such as would occur in either the scaling-droplet (SD)
McMillan (1984); Bray and Moore (1985); Fisher and Huse (1986, 1988) or
trivial–non-trivial (TNT) Krzakala and Martin (2000); Palassini and Young
(2000) pictures. Then any infinite sequence, chosen independently of the
couplings, will converge to a limiting thermodynamic state: namely, a mixed
state comprising the two globally reversed pure/ground states, each with
probability 1/2.

Suppose, on the other hand, that the Hamiltonian (1) supports infinitely many
incongruent Huse and Fisher (1987) (i.e., differing by relative interfaces
with dimensionality equal to the space dimension) pure/ground states, as
predicted by the RSB Parisi (1979, 1983); Mézard  _et al._ (1984a, b, 1987);
Marinari  _et al._ (2000); Read (2014) or chaotic pairs (CP) Newman and Stein
(1996a, 1997, 1998a) pictures. Then an infinite sequence of volumes chosen
independently of the couplings will generally not converge to an infinite-
volume (pure or mixed) thermodynamic state, a phenomenon referred to as
chaotic size dependence (CSD) Newman and Stein (1992, 2003). Instead, there
are infinitely many convergent subsequences of volumes, each of which
converges to a thermodynamic state different from the others. Although
compactness of the space of spin configurations ensures Chung (2001) that the
thermodynamic states of appropriately chosen subsequences of volumes will
converge, in general these subsequences must be chosen dependent on the
coupling realization.

Before proceeding, it may be worthwhile to step back and note that
nonconvergence of the above type can occur in any system, if the boundary
conditions are not chosen appropriately. Consider for example a simple uniform
Ising ferromagnet in zero field in two or higher dimensions. An infinite
sequence of volumes with periodic or free boundary conditions will of course
converge to a mixture of the positive and negatively magnetized states, each
with probability $1/2$. But suppose one chose instead fixed random boundary
conditions, where for every volume the spins at the boundary are chosen to be
$+1$ or $-1$ independently in accordance with a fair coin toss. In that case
the sequence will not converge, in the following sense: consider any fixed
volume $\Lambda_{W}$, and consider the spin configuration (at zero
temperature) or correlation functions (at positive temperature), within this
fixed volume, which we hereafter refer to as a “window”. Then at low
temperature, for a sequence of $\Lambda_{L}$’s with fixed random boundary
conditions, the thermodynamic state within the window $\Lambda_{W}$ will
continually flip between the positively and negatively magnetized states
(depending on whether the boundary condition on $\partial\Lambda_{L}$ has an
excess of plus or minus spins) as $L\to\infty$, never tending toward a limit.
In this case there are two subsequences of boundary conditions leading to
separate limits, one being the positively magnetized and the other the
negatively magnetized state.

Of course for the uniform ferromagnet we understand the nature of the broken
symmetry and the order parameter, and how to choose appropriate boundary
conditions to arrive at the desired thermodynamic limit. For the spin glass
and similar systems with quenched disorder, however, the situation is far more
complicated: even knowing the order parameter, whether the EA order parameter
Edwards and Anderson (1975) or the Parisi order parameter Parisi (1979, 1983);
Mézard  _et al._ (1984a, b, 1987), or something else, will not solve the
problem. We simply don’t know how to choose coupling-dependent boundary
conditions that lead to thermodynamic convergence, and even if we did, these
coupling-dependent boundary conditions would lead to non-measurability and
consequently no clear way of averaging over coupling realizations. Put another
way, there is likely no finite procedure for selecting convergent
subsequences, thereby generating thermodynamic states. CSD arising in this way
appears to be a fundamental property of Hamiltonians with quenched randomness,
and cannot be transformed away by any known methods.333CSD occurs also in the
SK model, though in the sense of overlap distributions; see Sect. 5.1 of
Newman and Stein (2003).

The metastate concept was introduced as a tool for handling these sorts of
situations; we turn to its definition, construction, and properties in the
next section.

##  III Metastates

###  III.1 Analogy to chaotic dynamical systems

The term “chaotic size dependence” provides more than a picturesque
description of the nonconvergence of correlation functions in a coupling-
independent sequence of volumes: there is a deep analogy to dynamical chaos
which provides a clue to resolving the difficulties posed by CSD. Consider the
chaotic orbit of a particle in a dynamical system: its behavior in time $t$ is
deterministic but effectively unpredictable, and can be treated as if it were
a random sampling from some distribution $\kappa$ on the particle’s space of
states. One way to construct $\kappa$ is through a histogram that records the
proportion of time the particle spends in each (coarse-grained) region of
states.

Similarly, for systems with quenched disorder such as spin glasses, the
behavior of correlation functions as $L$ changes is analogous to the
particle’s chaotic behavior in time $t$. Roughly speaking, the fraction of
$\Lambda_{L}$’s in which a given thermodynamic state appears can be shown to
converge as $L\to\infty$ Newman and Stein (1996b, 1998b); Newman (1997), just
as the $\kappa$ describing the chaotic behavior of a dynamical system
converges as $t\to\infty$. The resulting distribution over thermodynamic
states carries information on how often a given state appears within
$\Lambda_{L}$’s in the infinite sequence.

Strictly speaking, a thermodynamic state $\Gamma$ is an infinite-volume
quantity; by saying that it “appears” within a finite volume $\Lambda_{L}$, we
mean the following. Fix a window $\Lambda_{W}$ deep inside $\Lambda_{L}$,
i.e., with $W\ll L$; within this window all correlation functions computed
using the finite-volume Gibbs state $\rho^{(L)}$ are the same as those
computed using $\Gamma$ (with negligibly small deviations going to zero as
$L\to\infty$ with $W$ fixed).

###  III.2 Constructions of metastates

There are several kinds of metastate which carry different kinds of
information; here we focus on the simplest variety. A metastate depends on the
Hamiltonian (which we will hereafter assume to be (1)), the dimension $d$,
temperature $T$ (which can be zero Arguin  _et al._ (2016) or nonzero Arguin
_et al._ (2014)), disorder realization (corresponding here to $J$), and an
infinite sequence of volumes $\Lambda_{L}$ each with a specified boundary
condition. We will hereafter assume that every $\Lambda_{L}$ has periodic
boundary conditions unless otherwise specified. We will denote a resulting
“periodic boundary condition metastate” by $\kappa_{J}$, where dependence on
all other quantities is suppressed but understood.

There are two independent constructions of metastates, one initially
constructed for random-field magnets Aizenman and Wehr (1990) and one
initially constructed for spin glasses Newman and Stein (1996b). It was proved
in Newman and Stein (1998b) (see also Newman (1997)) that there exists an
infinite sequence of volumes, chosen independently of $J$, for which the two
constructions give the same metastate, and so either method can be used,
depending on which is more convenient to address the problem at hand. Both
constructions are sufficiently general that they can be used for a wide
variety of applications, including mean-field (MF) Curie-Weiss ferromagnets
with random couplings Külske (1998a, b), neural networks Külske (1998a); van
Enter and Schaap (2001), and other disordered systems.

The first construction, due to Aizenman and Wehr (AW) Aizenman and Wehr
(1990), uses a canonical ensemble approach based on varying the couplings
outside $\Lambda_{L}$. Let $\rho^{(L)}$ denote the finite-volume Gibbs state
(or ground state pair at zero temperature) in $\Lambda_{L}$ with periodic
boundary conditions, and consider for each $\Lambda_{L}$ the random pair
$(J_{L},\rho^{(L)})$, where $J_{L}$ is the restriction of $J$ to
$\mathbb{E}_{L}$, and take the limit (using compactness Chung (2001)) of these
finite-dimensional distributions along a $J$-independent subsequence of $L$’s.
This yields a probability distribution $\kappa$ on infinite-volume pairs
$(J,\Gamma)$, where $\Gamma$ denotes a thermodynamic state, which is
translation-invariant (because of the use of periodic boundary conditions)
under simultaneous lattice translations of $J$ and $\rho^{(L)}$. The metastate
is then the conditional distribution $\kappa_{J}$ of $\kappa$ given a fixed
$J$; it is supported on the infinite-volume thermodynamic states arising from
(sub)sequence limits of finite-volume Gibbs states.

The second construction, due to Newman and Stein (NS) Newman and Stein
(1996b), is motivated by the chaotic orbits analogy. Given an infinite
sequence of volumes $\Lambda_{L_{1}}$, $\Lambda_{L_{2}}$, …with $L_{1}\ll
L_{2}\ll\ldots$ such that $L_{k}\to\infty$ as $k\to\infty$, construct a
microcanonical ensemble $\kappa_{N}$ in which each of the finite-volume Gibbs
states (or ground states at zero temperature)
$\rho^{(L_{1})},\rho^{(L_{2})},\ldots,\rho^{(L_{N})}$ has weight $N^{-1}$. The
ensemble $\kappa_{N}$ converges to the metastate $\kappa_{J}$ as $N\to\infty$
(again, possibly with use of a subsequence) in the sense that, for every well-
behaved function $g(\cdot)$ on thermodynamic states $\Gamma$ (e.g., a function
on finitely many spins),

$\lim_{N\to\infty}N^{-1}\sum_{k=1}^{N}g(\rho^{(L_{k})})=\int g(\Gamma)\
d\kappa_{J}(\Gamma)\,.$ (3)

From (3) we see that $\int d\kappa_{J}(\Gamma)=1$ and $\kappa_{J}$ can
therefore be interpreted as a probability measure on thermodynamic states: the
finite-volume probability of any event depending on a finite set of spins
and/or couplings converges in the infinite-volume limit to the
$\kappa_{J}$-probability of that event.

###  III.3 Large finite volumes and the thermodynamic limit

A metastate reconciles how nonconvergence in the thermodynamic limit can be
used to provide information on the state of a typical macroscopically large
volume. The information contained in $\kappa_{J}$ includes the fraction of
cube sizes $L_{k}$ which the system spends in different infinite-volume
thermodynamic states  as $k\to\infty$. If there is only a single pair of
pure/ground states, as in SD McMillan (1984); Bray and Moore (1985); Fisher
and Huse (1988, 1986), then along any deterministic (i.e., not conditioned on
$J$) sequence of volumes the distribution of spin configurations in any fixed
window $\Lambda_{W}$ generated by the finite-volume Gibbs states will
eventually settle down to a fixed state, which as discussed above is the
restriction of the infinite-volume thermodynamic state to $\Lambda_{W}$.

On the other hand, if there are many infinite-volume thermodynamic states, as
in RSB Read (2014); Parisi (1979, 1983); Mézard  _et al._ (1984a, b, 1987);
Marinari  _et al._ (2000) (see below), then CSD occurs and the set of
correlation functions in $\Lambda_{W}$ never converges to a limit. Instead,
for any $\Lambda_{L_{k}}$ with $L_{k}$ sufficiently large, the set of all
correlation functions in $\Lambda_{W}$ will be identical to that of one of the
many (pure or mixed) infinite-volume thermodynamic states available for the
system to choose from, and the “chosen” state varies with $L_{k}$. Although
the correlation function values in $\Lambda_{W}$ never settle down, the
fraction of volumes $\Lambda_{L_{k}}$ in which a particular thermodynamic
states appears in $\Lambda_{W}$ does converge (along the subsequence) to a
limit, and this information is contained within $\kappa_{J}$ Newman and Stein
(1998b); Newman (1997).

###  III.4 Formal definition and covariance properties of a metastate

As mentioned earlier, there are more complex metastate constructions that
contain more information than the “simplest” metastate discussed above. A
particularly useful construct at zero temperature is a so-called excitation
metastate Newman and Stein (2000, 2001); Arguin  _et al._ (2010, 2019)), which
contains information not only on the probability of appearance of different
ground states but also all possible local excitations above each ground state.
We do not discuss these other metastates here but refer the reader to the
references for more details.

For the simple periodic boundary condition metastate, we saw in the previous
section two different but equivalent constructions. By themselves they do not
define a metastate, considering that there may be other constructions that
also lead to the same result. To arrive at a formal definition, we would like
a metastate to satisfy certain properties: it must be supported on (infinite-
volume) thermodynamic states of the Hamiltonian at positive temperature or
ground states at zero temperature, and it must be a probability distribution
on these thermodynamic states in the sense of (3). We would also like a
metastate to satisfy certain useful covariance properties, to be discussed
below.

We begin by formally defining thermodynamic states (i.e., infinite-volume
Gibbs states) as probability measures on spin configurations. Let
$\Sigma=\\{-1,+1\\}^{\mathbb{Z}^{d}}$ be the set of all infinite-volume spin
configurations and let $\mathcal{M}_{1}(\Sigma)$ be the set of (regular Borel)
probability measures on $\Sigma$. An infinite-volume Gibbs state $\Gamma$ for
the Hamiltonian $H_{J}$ (1) is an element of $\mathcal{M}_{1}(\Sigma)$ that
(at a given temperature $T$) satisfies the Dobrushin-Lanford-Ruelle (DLR)
equations Georgii (1988) for that Hamiltonian. Essentially, the DLR equations
require that the conditional probability in $\Gamma$ of the occurrence of an
event in any finite subregion $\Lambda_{L}$ (in particular, of a given
configuration of a finite set of spins ${\cal S}\subseteq\Lambda_{L}$),
conditioned on the spins outside $\Lambda_{L}$, be equal to that obtained from
the finite-volume Gibbs distribution $\rho^{(L)}$, using the boundary
condition on $\partial\Lambda_{L}$ determined by the spin configuration
outside $\Lambda_{L}$. An infinite-volume Gibbs state may be pure or mixed,
i.e., a convex combination of two or more pure states. (A pure state is a
Gibbs state that cannot be expressed as a convex combination of other Gibbs
states; see also Sec. IV.2.) We denote the set of infinite-volume Gibbs states
corresponding to the coupling realization $J$ by $\mathcal{G}_{J}$.

We can now formally define metastates $\kappa_{J}$ for the EA Hamiltonian (1)
on $\mathbb{Z}^{d}$ as follows:

Definition 3.1 A metastate $\kappa_{J}$ is a measurable mapping

$\displaystyle\mathbb{R}^{\mathbb{E}^{d}}$
$\displaystyle\to\mathcal{M}_{1}(\Sigma)$ (4) $\displaystyle J$
$\displaystyle\mapsto\kappa_{J}$

with the following properties:

  1. 1.

Support on Gibbs states. Every state sampled from $\kappa_{J}$ is a
thermodynamic state for the realization $J$:

$\kappa_{J}\Bigl{(}\mathcal{G}_{J}\Bigr{)}=1.$

  2. 2.

Coupling Covariance. For $B\subset\mathbb{Z}^{d}$ finite,
$J_{B}\in\mathbb{R}^{\mathbb{E}(B)}$ (the set of edges in $B$), and any
measurable subset $A$ of $\mathcal{M}_{1}(\Sigma)$, we define the operation
${\cal L}_{J_{B}}:\Gamma\mapsto{\cal L}_{J_{B}}\Gamma$ by its effect on
expectation $\langle\cdots\rangle_{\Gamma}$ in $\Gamma$,

$\left\langle f(\sigma)\right\rangle_{{\cal
L}_{J_{B}}\Gamma}=\frac{\left\langle f(\sigma)\exp\Bigl{(}-\beta
H_{J_{B}}(\sigma)\Bigr{)}\right\rangle_{\Gamma}}{\left\langle\exp\Bigl{(}-\beta
H_{J_{B}}(\sigma)\Bigr{)}\right\rangle_{\Gamma}}\,,$ (5)

which describes the effect of modifying the couplings within a finite subset
$B$ of ${\mathbb{Z}}^{d}$. We then require that the metastate be covariant
under local modifications of the couplings, i.e.,

$\kappa_{J+J_{B}}(A)=\kappa_{J}({\cal L}_{J_{B}}^{-1}A)$

where ${\cal L}_{J_{B}}^{-1}A=\Bigl{\\{}\Gamma\in\mathcal{M}_{1}(\Sigma):{\cal
L}_{J_{B}}\Gamma\in A\Bigr{\\}}$.

We further define a translation-covariant metastate:

Definition 3.2 A translation-covariant metastate is a metastate $\kappa_{J}$
with the additional property of

  1. 1.

Translation Covariance. For any translation $\tau$ of $\mathbb{Z}^{d}$ and any
measurable subset $A$ of $\mathcal{M}_{1}(\Sigma)$,

$\kappa_{\tau J}(A)=\kappa_{J}(\tau^{-1}A).$

These requirements, including translation covariance, are satisfied by both
the AW and NS periodic boundary condition metastate constructions. Translation
covariance follows easily in the thermodynamic limit from the torus-
translation-covariance of periodic boundary conditions. For other coupling-
independent boundary conditions, metastates can be obtained in a similar way
as for periodic boundary conditions, but are not necessarily translation
covariant; however, translation covariance can be recovered by taking an
average of the translates of the finite-volume Gibbs measures Newman and Stein
(2001); Arguin  _et al._ (2010). In practice, translation covariance is a
crucial property without which few conclusions can be drawn.

The importance of coupling covariance, i.e., covariance of a (possibly not
translation-covariant) metastate under a local modification of the couplings,
is less obvious. We note first that under any finite change in coupling
values, pure states transform to pure states (with some local changes in
correlation functions). An important consequence of coupling covariance is
that if a metastate is supported on some set of pure states, then under any
finite coupling transformation $J\to J^{\prime}$ the transformed metastate
$\kappa_{J^{\prime}}$ will be supported on the same set of pure states as
$\kappa_{J}$ (modulo finite changes in certain correlation functions due to
the coupling transformation). Coupling covariance is useful because it helps
make explicit the dependence of the metastate on any finite set of couplings.
This is crucial because numerous applications, particularly taking derivatives
with respect to coupings, require taking into account the dependence of free
energies or other thermodynamic quantities on all of the couplings. Moreover,
covariance with respect to changes of couplings ensures that taking
derivatives with respect to couplings will not generate difficult-to-handle
boundary terms induced by pure states flowing into or out of a region of
integration.

####  III.4.1 Gauge invariance of metastates

Definition 3.1, along with the equivalence of AW and NS metastates obtained
along certain deterministic sequences of volumes, leads to another covariance
property of such metastates, proved in Newman and Stein (1998a), namely that
they display gauge invariance with respect to boundary conditions. Recall that
constructing a metastate requires specification in advance of the boundary
conditions on an infinite sequence of volumes. The theorem of Newman and Stein
(1998a) shows that a metastate is unchanged by any gauge transformation along
all or any subset of boundaries $\partial\Lambda_{L}$ in the sequence. In
other words periodic and antiperiodic boundary condition metastates are
identical, as are any metastates constructed with any combination of periodic
and antiperiodic boundary conditions within or between volumes. Similarly, all
fixed boundary condition metastates are equivalent, for any way of determining
the fixed boundary conditions, so long as it is independent of $J$.

This leads to several powerful conclusions; among them is that our restriction
in this contribution to periodic boundary condition metastates incurs little
loss of generality. What is important is that Definition 3.1 almost certainly
requires the use of coupling-independent boundary conditions to generate
metastates. The use of coupling-dependent boundary conditions, though
potentially useful for some purposes, would almost certainly result in losing
the essential properties of measurability, translation-covariance, and
coupling-covariance.

Now that metastates have been defined and their basic properties enumerated,
we can turn to some applications. Our main interest in this contribution is
how the metastate construct can be used to specify what RSB can (and cannot)
mean in the EA model in finite dimensions. We turn to this question in the
next section.

##  IV Metastates and pure states

Next we will make some further definitions, and then describe distinct classes
of scenarios that are possible within the metastate framework. The entire
discussion in this section is for infinite-volume systems.

###  IV.1 Trivial versus non-trivial metastates and the metastate-average
state (barycenter)

We begin by defining what we will mean by a “trivial metastate”. We say a
metastate is trivial if $\kappa_{J}$ consists of a single atom, in other words
a point mass or $\delta$-function on a single Gibbs (i.e. DLR) state;
otherwise, it is nontrivial or “dispersed”.

An additional construction based on the metastate will be used. Given a
metastate, the Gibbs states can be averaged using the metastate, to produce
the metastate-average state (MAS) or barycenter $\rho_{J}(\sigma)$ Aizenman
and Wehr (1990); Newman and Stein (1996b, 1997, 2003), which is itself a Gibbs
state (though a rather special one), and which still depends on the sample of
$J$. That is, $\rho_{J}(\sigma)=[\Gamma(\sigma)]_{\kappa_{J}}$, where the
square bracket $[\cdots]_{\kappa_{J}}$ denotes metastate average, in other
words

$\rho_{J}(\sigma)=\int\Gamma(\sigma)d\kappa_{J}(\Gamma).$ (6)

If the metastate is trivial, then the MAS $\rho_{J}=\Gamma$, the Gibbs state
on which $\kappa_{J}$ is concentrated.

###  IV.2 Pure-state decomposition of a Gibbs state

A convex combination (or mixture) of distinct Gibbs states for the same
$H_{J}$ is again a Gibbs state; in general, a convex combination could involve
an average taken using a probability measure on Gibbs states. A Gibbs state
that cannot be expressed as such a combination of other Gibbs states is called
an extremal or pure (Gibbs) state. Any Gibbs state can be decomposed uniquely
into a mixture of pure states for the same $H$ Georgii (1988); Simon (1993),
in the form

$\Gamma(\sigma)=\sum_{\alpha}w_{\Gamma}(\alpha)\Gamma_{\alpha}(\sigma),$ (7)

where $\Gamma_{\alpha}$ are pure states, and for any given $\Gamma$ the
weights $w_{\Gamma}(\alpha)\geq 0$ sum to $1$. [In practice the decomposition
might be continuous and the sum over $\alpha$ would be replaced with an
integral using the measure $w_{\Gamma}(\alpha)$, but we will usually not show
this explicitly.] Pure states can also be characterized in other, more
intrinsic, ways; one of these is that they exhibit clustering of correlations
Georgii (1988); Simon (1993). For example, a simple consequence of clustering
is that if the Gibbs state $\Gamma$ is pure then

$\langle\sigma_{\bf x}(\tau^{-1}\sigma)_{\bf
x}\rangle_{\Gamma}\to\langle\sigma_{\bf
x}\rangle_{\Gamma}\langle(\tau^{-1}\sigma)_{\bf x}\rangle_{\Gamma}$ (8)

as $|{\bf x}-\tau({\bf x})|\to\infty$ [again, $(\tau^{-1}\sigma)_{\bf
x}=\sigma_{\tau({\bf x})}$ is a translation of $\sigma_{\bf x}$], and a
similar statement holds for any two local functions $A(\sigma)$ and
$B(\sigma)$ built from spins near ${\bf x}$ taking the place of the two
$\sigma_{\bf x}$’s. The pure states correspond to the “ordered” states, which
should be familiar in cases such as the uniform Ising ferromagnet in zero
magnetic field, which at low temperature (for dimension $\geq 2$) has (at
least) two ordered states; in one of these, the majority of spins in any
finite region in a typical configuration are $+1$, while in the other, they
are $-1$.

The Hamiltonian $H_{J}$ for a classical spin system with quenched disorder may
possess a nontrivial group of global “internal” symmetry operations that act
locally on all spins simultaneously, and leave $H_{J}$ invariant. For example,
the EA Hamiltonian in eq. (1) is invariant under a global spin flip which maps
$\sigma_{\bf x}\to-\sigma_{\bf x}$ for all ${\bf x}$. If a magnetic field were
included, by adding $-h\sum_{\bf x}\sigma_{\bf x}$ to $H_{J}$, this symmetry
would be lost. From here on we will assume that a metastate is obtained using
a boundary condition that respects the internal symmetry of the Hamiltonian,
in which case the metastate constructions preserve any such symmetries: a
Gibbs state drawn from a metastate $\kappa_{J}$ will be invariant under the
full symmetry group of the Hamiltonian, which in turn implies symmetries (Ward
identities) of correlation functions in that Gibbs state.

In the high-temperature region, there is a unique Gibbs state, which is pure,
and the metastate is unique and trivial. In this case, the pure state is
invariant under the symmetry operations. In a low-temperature phase, it may be
that a given pure state is not invariant under the full symmetry group (when
the group is nontrivial), in which case it is said to exhibit spontaneous
symmetry breaking Georgii (1988); Simon (1993). In general, the pure states
can be partitioned into orbits, defined such that the symmetry group acts
transitively on each orbit. Each orbit consists either of a single invariant
pure state, or of more than one pure state. As a Gibbs state drawn from the
metastate $\kappa_{J}$ must be invariant under the symmetry, the pure-state
decomposition of a Gibbs state drawn from $\kappa_{J}$ must consist of
symmetry orbits. We will call a Gibbs state trivial if its pure-state
decomposition consists of a single symmetry orbit, and nontrivial if it
consists of more than one. (While spontaneous symmetry breaking is certainly a
nontrivial phenomenon, and is expected to occur at low temperature in many
spin-glass systems with nontrivial symmetry groups, this terminology will be
useful for our analysis of spin-glass phases.) Thus this definition using
symmetry orbits can handle all cases, including high and low temperature, both
discrete and continuous symmetries, and zero and nonzero magnetic field.

For the MAS, we also define its pure-state decomposition

$\displaystyle\rho_{J}(\sigma)$ $\displaystyle=$
$\displaystyle\sum_{\Gamma}\kappa_{J}(\Gamma)\sum_{\alpha\in\Gamma}w_{\alpha}(\Gamma)\Gamma_{\alpha}(\sigma)$
(9) $\displaystyle=$
$\displaystyle\sum_{\alpha}\mu(\alpha)\Gamma_{\alpha}(\sigma),$ (10)

where the weights are
$\mu(\alpha)=\sum_{\Gamma}\kappa_{J}(\Gamma)w_{\Gamma}(\alpha)=[w_{\Gamma}(\alpha)]_{\kappa_{J}}$,
and implicitly depend on $J$. Here again we have written sums, though some of
these may in fact be integrals over a probability measure. The same
distinction between trivial and nontrivial decompositions of a Gibbs state can
be applied to this one.

###  IV.3 Classes of scenarios for spin glasses

The preceding definitions of trivial and nontrivial, both for metastates and
for Gibbs states, now lead to four combinations Newman and Stein (1996b,
1997); Read (2014) that provide a broad-brush classification of scenarios for
short-range spin glasses, such as the EA model. To simplify the description,
we assume the Gibbs states have the same character (i.e., either trivial or
non-trivial) for almost every Gibbs state drawn from the metastate (it is
unknown whether this must be true, however, the same conclusions would hold if
we replaced “almost surely non-trivial Gibbs states” with “nonzero probability
of non-trivial Gibbs states”). As far as is known to the authors at the time
of writing, all four classes remain open as possibilities that may occur in
some model systems, possibly in different dimensions of space, or at different
parameter values, for example different temperatures or magnetic fields. They
are:
1) Both the metastate and the Gibbs states are trivial. This is the case in
the important SD scenario McMillan (1984); Bray and Moore (1985); Fisher and
Huse (1986, 1988), as assumed either implicitly McMillan (1984); Bray and
Moore (1985), or (for the Gibbs states) explicitly, following some discussion
Fisher and Huse (1988). Scenarios in this class still allow for rich physical
phenomena at low temperature Fisher and Huse (1988). It is also the case that
occurs at high temperature.
2) Nontrivial metastate and trivial Gibbs states. NS Newman and Stein (1996a,
b, 1997) termed this class chaotic pairs (CPs), referring to the Ising EA
model at zero magnetic field (the pairs being the symmetry orbits), while the
term “chaotic singles” has been used for the same but with a magnetic field
Read (2014).
3) Trivial metastate and nontrivial Gibbs states. There is no widely accepted
name for this class of scenarios, but below we will mention its connection
with the literature. (NS Newman and Stein (1996a, b, 1997) originally termed
one scenario in this class the standard SK picture, but we discuss this point
in the following subsection.)
4) Both the metastate and the Gibbs states are nontrivial. This class is the
broadest, and its instances have the richest structure. Below we review
arguments that the scenario predicted by RSB Parisi (1979, 1983); Mézard  _et
al._ (1984a, b, 1987); Marinari  _et al._ (2000); Read (2014) is in this class
in general. (NS Newman and Stein (1996a, b, 1997) originally referred to a
scenario in this class as non-standard SK.)
It follows from this classification that in all classes, except class 1), the
MAS is a nontrivial Gibbs state.

Although we mainly discuss positive temperature, it may be useful to say
something about $T=0$ also. If the probability distribution of the bonds $J$
is continuous without atoms (so not, for example, the bimodal model, in which
$J_{{\bf xy}}=\pm J_{0}$ for $J_{0}>0$ constant), as assumed so far and as we
will continue to assume in this paragraph, then at zero temperature a pure
state is a ground state (with probability one), that is a single spin
configuration. For classical spins, ground states always form non-trivial
symmetry orbits. Some reflection on the metastate constructions shows that at
$T=0$ the Gibbs states drawn from the metastate are trivial with probability
one. Hence, at $T=0$, only classes 1) and 2) can occur. We may expect that if
3) or 4) occur at $T>0$, then as $T\to 0$, assuming that no other phase
transition intervenes, a scenario in class 3) will reduce to one in 1), and
one in 4) to one in 2).

###  IV.4 Constraints on the scenarios

Next we discuss some rigorous constraints on the behavior of the scenarios
within the preceding classification.

####  IV.4.1 Non-self-averaging of overlap distributions

Here we will need to refer to the overlaps of pure states, and to a version
Newman and Stein (1997) of Parisi’s overlap distribution function $P(q)$
Parisi (1983) for the EA model. Given two pure states $\Gamma_{\alpha}$,
$\Gamma_{\alpha^{\prime}}$, we can define the average overlap
$q_{\alpha\alpha^{\prime}}$ in the window $\Lambda_{W}$, and take
$W\to\infty$:

$q_{\alpha\alpha^{\prime}}=\lim_{W\to\infty}\frac{1}{W^{d}}\sum_{x\in\Lambda_{W}}\langle\sigma_{\bf
x}\rangle_{\alpha}\langle\sigma_{\bf x}\rangle_{\alpha^{\prime}}.$ (11)

When $(J,\Gamma)$ are drawn from the distribution $\kappa$, which is
translation invariant, and $\Gamma_{\alpha}$, $\Gamma_{\alpha^{\prime}}$ are
drawn from the distribution $w_{\Gamma}(\alpha)w_{\Gamma}(\alpha^{\prime})$ on
pairs of pure states in $\Gamma$ [and note that
$\kappa(J,\Gamma)w_{\Gamma}(\alpha)w_{\Gamma}(\alpha^{\prime})$ is translation
invariant], then the ergodic theorem (for translations) implies that the
$W\to\infty$ limit of the random variable $q_{\alpha\alpha^{\prime}}$ exists
almost surely; the limit is translation invariant. Similarly to Parisi Parisi
(1983), we then define the probability distribution of overlaps for
$(J,\Gamma)$ by

$P_{J,\Gamma}(q)=\sum_{\alpha,\alpha^{\prime}}w_{\Gamma}(\alpha)w_{\Gamma}(\alpha^{\prime})\delta(q-q_{\alpha\alpha^{\prime}}).$
(12)

We define $P_{J}(q)=[P_{J,\Gamma}(q)]_{\kappa_{J}}$ to be the metastate
average of $P_{J,\Gamma}(q)$ over $\Gamma$; thus if the metastate is trivial,
$P_{J}(q)=P_{J,\Gamma}(q)$. Finally, writing $[\cdots]_{\nu}$ for expectation
using the distribution $\nu(J)$ on $J$, we define $P(q)=[P_{J}(q)]_{\nu}$. For
the MAS $\rho_{J}$, we define $P_{J,\rho}(q)$ similarly to $P_{J,\Gamma}$,
using the pure-state decomposition of $\rho_{J}$ introduced earlier. Each of
these distributions is invariant under translations of $(J,\Gamma)$. We note
that for the EA model at zero magnetic field, the overlaps for the two members
of a non-trivial symmetry orbit with another fixed pure state have opposite
signs, and then $P_{J,\Gamma}(q)$ is an even function of $q$.

A feature of the RSB MF theory obtained for the infinite-range SK model is
that, while the thermodynamic functions such as the free energy per spin self-
average (that is, its thermodynamic limit exists for given disorder, and is
almost surely a constant without fluctuations due to the disorder) Guerra and
Toninelli (2002); Guerra (2003); Talagrand (rlin), the same is not true for
the distribution of overlaps in a given sample of the SK model Young  _et al._
(1984); Mézard  _et al._ (1984b). In contrast, for short-range models such as
the EA model, both $P_{J}(q)$ and $P_{J,\rho}(q)$ must self average, as a
consequence of translation invariance and the ergodicity of the distribution
$\nu(J)$ Newman and Stein (1996a). Consequently, the picture of RSB in short-
range models as belonging to class 3), as seemed to be assumed within the
standard interpretations of that time, was shown to be inconsistent with the
behavior that would be expected on the basis of the RSB MF theory. Given the
ergodicity of $\nu$, what was called non-self-averaging (NSA) behavior is
permitted in a short-range model only for $P_{J,\Gamma}(q)$ and only if
$\kappa$ is not ergodic, and hence only if $\kappa_{J}$ is non-trivial Newman
and Stein (1996b, 1997). Further, within the RSB scheme (discussed further
below), NSA behavior can arise only within class 4), with both a non-trivial
metastate and non-trivial Gibbs states Read (2014).

####  IV.4.2 Cardinality and structure of pure-state decompositions

The results of the preceding subsection can be refined further. Using both
coupling covariance and translation covariance, it was shown in Ref. Newman
and Stein (2006) that if the Gibbs states drawn from the metastate are
nontrivial, then the pure-state decomposition of the MAS must be uncountable,
forming a continuum without any atoms. Hence for scenarios in class 4), such
as RSB, if the Gibbs states have a decomposition into a countable number of
symmetry orbits, as suggested by RSB, then the metastate not only must be non-
trivial but also must be a continuous atomless distribution on Gibbs states.
Also, in class 3) the nontrivial Gibbs state must have a decomposition into a
continuous atomless distribution on pure states. On the other hand, for
scenarios in class 2), such as CPs, the metastate is permitted to have support
on a countable or even a finite number ($>1$) of trivial Gibbs states. We
suspect that here too for spin glasses the metastate must be supported on
infinitely many such Gibbs states, but as of now the question remains open.

####  IV.4.3 Complexity of Gibbs states and metastates

A further refinement invokes the idea of complexity of a pure-state
decomposition. In the form suggested by Palmer Palmer (1982), complexity is
the entropy of the set of weights $w_{\Gamma}(\alpha)$ occurring in the
decomposition of a given $\Gamma$; crudely, this corresponds to the logarithm
of the number of pure states, at least if they all have approximately the same
weight. This notion of entropy is well defined, but less useful if it turns
out to be infinite. In the latter case, Palmer suggested calculating it in
finite size, and that it might even be extensive, increasing proportionally to
the volume as the system size increases. But this runs into the problem that,
in finite size, pure states and hence the decomposition are not well defined.
Ignoring that difficulty for a moment, it was recognized that the complexity
cannot be extensive van Enter and van Hemmen (1984). Working with well-defined
pure states in infinite size, the idea for a given Gibbs state of counting the
number of pure states that can be distinguished when examining the spin
configuration in a finite window was proposed Newman and Stein (1996b, 1997,
2003). For a (hyper-)cubic window $\Lambda=\Lambda_{W}$ of side $W$, centered
say at the origin, it was expected in a short-range model that the logarithm
of the number of distinguishable pure states could not grow with $W$ faster
than $W^{d-1}$ Newman and Stein (1996b, 1997, 2003). For Ising spins with
nearest-neighbor interactions, this arises from counting the number of
boundary conditions on the hypercube that could produce different ground
states, and the expectation that the same bound would apply for $T>0$.

All these issues can be dealt with in a well-defined way by using some
information theory Cover and Thomas (2006). A pure state corresponds to a
boundary condition, effectively at infinity. The weights and the pure states
allow us to consider the joint distribution
$w_{\Gamma}(\alpha)\Gamma_{\alpha}(\sigma)$ of pure states $\alpha$ and spin
configurations $\sigma$ (for given $\Gamma$). If we restrict $\sigma$ to the
spin configuration $\sigma_{\Lambda}$ in the hypercube $\Lambda=\Lambda_{W}$
as before, then the mutual information $I(\sigma_{\Lambda};\alpha)\geq 0$
between $\sigma_{\Lambda}$ and $\alpha$ can be defined in terms of the joint
distribution of those random variables, and we now call this
$\Lambda_{W}$-dependent quantity the complexity $K_{\Gamma}(\Lambda_{W})$ of
the Gibbs state $\Gamma$ Höller and Read (2020); Read (2022):

$K_{\Gamma}(\Lambda_{W})\equiv
I_{\Gamma}(\sigma_{\Lambda};\alpha)=\sum_{\sigma_{\Lambda},\alpha}w_{\Gamma}(\alpha)\Gamma_{\alpha}(\sigma_{\Lambda})\ln\frac{w_{\Gamma}(\alpha)\Gamma_{\alpha}(\sigma_{\Lambda})}{w_{\Gamma}(\alpha)\Gamma(\sigma_{\Lambda})}$
(13)

[recall that
$\Gamma(\sigma_{\Lambda})=\sum_{\alpha}w_{\Gamma}(\alpha)\Gamma_{\alpha}(\sigma_{\Lambda})$].
This definition does not use any notion of distinguishability of pure states
within a window. As $W\to\infty$, $K_{\Gamma}(\Lambda_{W})$ increases
monotonically, and tends to the same value as for Palmer’s definition of
complexity. We can similarly define the complexity $K_{\rho_{J}}(\Lambda_{W})$
of the MAS $\rho_{J}$ in place of $\Gamma$, and also the complexity
$K_{\kappa_{J}}(\Lambda_{W})$ of the metastate $\kappa_{J}$ itself:
$K_{\kappa_{J}}(\Lambda_{W})$ is defined as the mutual information between the
Gibbs (not pure) state $\Gamma$ and $\sigma_{\Lambda}$, using the joint
distribution $\kappa_{J}(\Gamma)\Gamma(\sigma_{\Lambda})$ Read (2022). These
three quantities are related (for given $J$ and $\Lambda_{W}$) by

$K_{\rho_{J}}(\Lambda_{W})=\left[K_{\Gamma}(\Lambda_{W})\right]_{\kappa_{J}}+K_{\kappa_{J}}(\Lambda_{W}).$
(14)

For any of these complexities, if it (or its expectation using $\kappa_{J}$ or
$\kappa$) diverges as $W\to\infty$, then we can look for the manner in which
it diverges. If the leading behavior is $\sim W^{d-\zeta^{\prime}}$, then
$d-\zeta^{\prime}>0$ indicates an uncountable number of pure states (or Gibbs
states in the case of $K_{\kappa_{J}}$). These definitions help to quantify
the structures in the above classification, and the exponents
$d-\zeta^{\prime}$ are expected to be universal within a phase, for example,
independent of temperature.

In the nearest-neighbor EA model of Ising spins in eq. (1), a simple argument,
similar to counting boundary conditions, then produces a bound on the
$\kappa$-expectation of the complexity $K_{\Gamma}(\Lambda_{W})$ of the Gibbs
states by a constant times the surface area $W^{d-1}$ Höller and Read (2020).
More generally, a bound of the same form can be obtained for short-range
models of classical spins Read (2022), though the argument was valid only for
$T>0$. (A corresponding result for long-range models was also obtained.) The
same bounds also apply to the expectations of the other two complexities; thus
$\zeta^{\prime}\geq 1$. Within the usual interpretation of RSB theory, the
expected complexity of a Gibbs state as $W\to\infty$ is
$\lim_{W\to\infty}\left[\left[K_{\Gamma}(\Lambda_{W})\right]_{\kappa_{J}}\right]_{\nu}=\psi(1)-\psi(1-x_{1})$
Gross and Mézard (1984), where $\psi$ is the digamma function and $1-x_{1}$ is
the weight in the $\delta$-function in $P(q)$ at $q(1)$, interpreted as
$1-x_{1}=\left[\left[\sum_{\alpha}w_{\Gamma}(\alpha)^{2}\right]_{\kappa_{J}}\right]_{\nu}\leq
1$.

###  IV.5 Correlations in the MAS

Like other states, the MAS can be characterized by its spin correlations, and
it will be useful for what follows to introduce one. Due to the quenched
disorder, correlation functions in a spin glass will be random with zero mean,
and the basic quantities are expectations of squared correlations. Thus we
define for Ising spins Read (2014); Höller and Read (2020)

$\displaystyle C_{\rm MAS}({\bf x},{\bf y})$ $\displaystyle=$
$\displaystyle\left[\left(\langle\sigma_{\bf x}\sigma_{\bf
y}\rangle_{\rho_{J}}-\langle\sigma_{\bf x}\rangle_{\rho_{J}}\langle\sigma_{\bf
y}\rangle_{\rho_{J}}\right)^{2}\right]_{\nu}$ (15) $\displaystyle=$
$\displaystyle\left[\left(\left[\langle\sigma_{\bf x}\sigma_{\bf
y}\rangle_{\Gamma}\right]_{\kappa_{J}}-\left[\langle\sigma_{\bf
x}\rangle_{\Gamma}\right]_{\kappa_{J}}\left[\langle\sigma_{\bf
y}\rangle_{\Gamma}\right]_{\kappa_{J}}\right)^{2}\right]_{\nu}$ (16)

using the definition of $\rho_{J}$. When spin-flip symmetry is present this
simplifies, because then $\langle\sigma_{\bf x}\rangle_{\Gamma}=0$ for all
${\bf x}$. This expression resembles the basic correlation function $\chi$ for
a spin glass, that is, taking account of the metastate,

$\chi({\bf x},{\bf y})=\left[\left[\left(\langle\sigma_{\bf x}\sigma_{\bf
y}\rangle_{\Gamma}-\langle\sigma_{\bf x}\rangle_{\Gamma}\langle\sigma_{\bf
y}\rangle_{\Gamma}\right)^{2}\right]_{\kappa_{J}}\right]_{\nu},$ (17)

but $C_{\rm MAS}$ differs in that the square is taken after the metastate
averages. If the metastate is trivial, this of course makes no difference, so
$C_{\rm MAS}=\chi$. In the case with spin-flip symmetry, $\chi({\bf x},{\bf
y})\geq C_{\rm MAS}({\bf x},{\bf y})$. The translation invariance of $\kappa$
implies that these correlation functions depend on ${\bf x}-{\bf y}$, not
${\bf x}$, ${\bf y}$ separately.

In a spin glass phase, if the decomposition of $\Gamma$ contains more than one
pure state (due to either spontaneous symmetry breaking or $\Gamma$ being
nontrivial), then in most scenarios (at $T>0$) $\chi$ tends to a non-zero
constant as $|{\bf x}-{\bf y}|\to\infty$. For a non-trivial metastate, $C_{\rm
MAS}({\bf x},{\bf y})$ may decay to values $\leq\chi({\bf x},{\bf y})$,
especially in those cases in which the MAS has an uncountable pure-state
decomposition. There are models Newman and Stein (1994); White and Fisher
(2006); Jackson and Read (2010) in which $C_{\rm MAS}$ tends to zero and the
leading asymptotic behavior is a power law:

$C_{\rm MAS}({\bf x}-{\bf y})\sim|{\bf x}-{\bf y}|^{-(d-\zeta)}$ (18)

(times a non-universal constant in general) as $|{\bf x}-{\bf y}|\to\infty$,
and $\zeta$ is known. In the scenarios of class 3), $\rho_{J}=\Gamma$, and
$\chi=C_{\rm MAS}$ itself may behave in the fashion described here. In
general, when the power law form holds, we expect $d-\zeta$ to be a universal
constant within a phase. Heuristic arguments in the EA model, to be discussed
in the following section, further support parts of this picture. The relation
$\zeta=\zeta^{\prime}$ was proposed, based on those models and a fractal
picture of ground states Read (2014), but it is not known if it must always
hold.

##  V Metastate interpretation of replica symmetry breaking

In this Section we explicitly connect RSB and metastate concepts and results.
The final arguments are at a theoretical physics level of rigor that involve
the use of RSB MF theory and fluctuations within replica field theory in
short-range models.

We begin with a technical point: even when the magnetic field is nominally
zero, the development of RSB appears to assume that an “ordering field” is
present, that is a small magnetic field that is taken to zero as the system
size tends to infinity, but sufficiently slowly that spin-flip symmetry of the
EA model is absent in the states in the limit. It then turns out that the
overlaps $q_{\alpha\alpha^{\prime}}$ inferred from the RSB MF theory are
almost surely nonnegative. In the following, we will assume this holds.

###  V.1 $q(x)$ function

We can define a function $q(x)$ from first principles, motivated by the
interpretation of RSB results Parisi (1983); Mézard  _et al._ (1984b, 1987);
Read (2014). From $P(q)=P_{J}(q)$ we first define the cumulative distribution
$x(q)=\int_{0^{-}}^{q}P(q^{\prime})\,dq^{\prime}$ of $q$ as a function on the
support of $P(q)$. On the graph of $x(q)$ versus $q$ we now interchange the
axes to produce the graph of a function $q(x)$ for $x\in(0,1)$, with the
convention that a jump in $x(q)$ becomes a constant region of $q(x)$, while a
constant region of $x(q)$ becomes a jump discontinuity of $q(x)$; $q(x)$ is a
monotonically-increasing function for $x\in(0,1)$. [The value assigned to
$q(x)$ at a discontinuity in the open interval $(0,1)$ has no known
significance.] We further define
$q(1)=\left[\left[\sum_{\alpha}w_{\Gamma}(\alpha)q_{\alpha\alpha}\right]_{\kappa_{J}}\right]_{\nu}=\left[\sum_{\alpha}\mu(\alpha)\langle\sigma_{\bf
x}\rangle^{2}_{\alpha}\right]_{\nu}$ (for any ${\bf x}$), where the self-
overlap $q_{\alpha\alpha}$, or equivalently the average over positions of the
time autocorrelation of single spins Edwards and Anderson (1975), is the EA
order parameter $q_{\rm EA}$ in the pure state $\Gamma_{\alpha}$. (We have
recently shown that $q_{\rm EA}$ must be the same for all pure states in the
decomposition of a given $\Gamma$ drawn from the metastate Newman  _et al._ ,
but whether it is the same for all $\Gamma$ is not known for short-range
systems.) Similarly, motivated by RSB arguments Read (2014) given below, we
define
$q(0)=\left[\sum_{\alpha,\alpha^{\prime}}\mu(\alpha)\mu(\alpha^{\prime})q_{\alpha\alpha^{\prime}}\right]_{\nu}=\left[\langle\sigma_{\bf
x}\rangle^{2}_{\rho_{J}}\right]_{\nu}$. Then $q(x)$ is now defined on $[0,1]$,
and may be discontinuous at $x=0$ or $1$. For $q(1)$, this can occur if the
measure $w_{\Gamma}(\alpha)$ is atomless (i.e. continuous) for a set of
$(J,\Gamma)$ of nonzero $\kappa$ probability; for those $\Gamma$, the
probability $w_{\Gamma}(\alpha)^{2}$ of drawing a pure state $\Gamma_{\alpha}$
twice from $\Gamma$ is zero, and so the $\nu(J)\mu(\alpha)$-expected self-
overlap $q(1)$ could differ from (in particular, exceed) the supremum of the
support of $P(q)$. For $q(0)$, the pairwise overlaps in pure states in the MAS
could differ from those in the pure states in a single $\Gamma$, with nonzero
$\kappa$ probability, and so $q(0)$ could differ from (even lie below) the
infimum of the support of $P(q)$. It is not clear whether $q(x)$ as we have
defined it must necessarily be monotonic for all $x\in[0,1]$; in any case, it
contains more information than $P(q)$ does.

We remark that, defining
$\overline{q}=\int_{0}^{1}dx\,q(x)=\int_{0^{-}}^{1^{+}}dq^{\prime}\,P(q^{\prime})q^{\prime}$,
$q(x)$ obeys $q(1)-\overline{q}\geq 0$ and $\overline{q}-q(0)\geq 0$, which
follow from the definitions (both quantities are sums of variances). If
$q(1)-\overline{q}>0$ then the Gibbs states drawn from the metastate are
nontrivial, while if $\overline{q}-q(0)>0$ then the metastate is nontrivial.
If $q(x)$ were monotonic on $[0,1]$, as it is in RSB theory, then both of
these would hold whenever $q(x)$ is not constant on $(0,1)$ Read (2014). The
Cauchy-Schwarz inequality yields $q_{\alpha\alpha^{\prime}}\leq
q_{\alpha\alpha}^{1/2}q_{\alpha^{\prime}\alpha^{\prime}}^{1/2}$, which would
give monotonicity at $x=1$ if $q_{\alpha\alpha}$ were independent of $\alpha$,
so that $q_{\alpha\alpha}=q_{\alpha^{\prime}\alpha^{\prime}}$ for all
$\alpha$, $\alpha^{\prime}$.

Next we briefly review the basics of RSB (for more detail, see Refs. Parisi
(1979, 1983); Mézard  _et al._ (1984b, 1987)). At MF level, the free energy
per spin has to be obtained by maximizing a functional of the symmetric matrix
$Q_{ab}$, where $a$, $b$ index the “replicas” and run from $1$ to $n$,
$Q_{aa}=0$ for all $a$, and the functional must be evaluated in the limit
$n\to 0$ before maximizing. (The maximization rather than the usual
minimization is a feature of replicas, and is due to the $n\to 0$ limit.)
Parisi’s hierarchical scheme Parisi (1979) for breaking the symmetry $S_{n}$
of the functional under permutations of the replicas is to first divide $n$
into $n/m_{1}$ blocks of size $m_{1}$, and assign the value $q_{0}$ to
elements in the off-diagonal blocks of $Q_{ab}$. Each diagonal block is then
divided into sub-blocks of size $m_{2}$ in the same fashion (the same for each
one), with the value $q_{1}$ assigned to its off-diagonal sub-blocks. This
procedure is iterated at most a countably-infinite number of times. Finally,
when $n\to 0$, the block sizes are assumed to obey $0=n\leq m_{1}\leq
m_{2}\leq\cdots\leq 1$, giving a countable partition of $[0,1]$, and
maximization of the functional is performed with respect to all possible such
partitions. At $n=0$, for a partition, $(q_{r})_{r}$ can be viewed as a
function $q(x)$, where $q(x)=q_{r}$ when $m_{r}<x<m_{r+1}$ (where
$m_{0}=n=0$), $q(0)=q_{0}$, and $q(1)$ is defined as $q_{r_{\rm max}}$ (where
$r_{\rm max}$ is the largest $r$) for a finite partition, and as
$q(1)=\limsup_{r\to\infty}q_{r}$ for an infinite one. $q(x)$ is found to be a
non-negative monotonically increasing function of $x\in[0,1]$; in many
situations, $q(x)$ is continuous on the open interval $(0,1)$. Its
interpretation for $x>0$, similar to that above, has long been standard Parisi
(1983).

$q(x)$ as obtained from RSB is not always continuous at $x=0$ or $1$, which
brings up an important point. The Parisi functional essentially consists of
multiple integrals over $x$’s of (products of) $q(x)$, so the maximization
over functions $q(x)$ is unaffected by changes in the value of $q(x)$ on sets
of $x$ of measure zero; hence the values $q(0)$, $q(1)$ do not affect the free
energy per spin. Sometimes in a given problem, distinct stationary points
$q(x)$ are found which differ only in the value of $q(1)$ or $q(0)$, and still
have the same MF free energy per spin; we have seen that the values of these
are important for the interpretation.

As particular examples, in class 1), SD, $q(x)$ would be constant. In the
scenarios of class 2), such as CPs, which with an ordering field reduce to
chaotic singles, pairwise- and self-overlaps of pure states in the MAS should
be different, while $P(q)$ consists of a single $\delta$-function at $q_{\rm
EA}$, and then $q(0)\neq q(x)=q(1)$ for $x>0$. Similarly, scenarios in class
3) can have $q(1)\neq q(x)=q(0)$ for $x<1$ (see the next paragraph). In class
4), there are many possible forms, including $q(x)$ continuous but not
constant on $[0,1]$. Thus examples in all four classes of scenarios can be
characterized by a $q(x)$ function, and examples even of classes 2) Höller and
Read and 3) (see below) can be found at MF level in RSB. However, in the
absence of rigorous results stating that $q(x)$ is monotonic on $[0,1]$, it is
not possible to say that if $q(x)$ is non-constant on $(0,1)$ then both the
Gibbs states and the metastate are non-trivial (as is the case in RSB);
apparently $P(q)=P_{J}(q)$ could be non-trivial (not a single
$\delta$-function) due to fluctuations in only one of them.

In MF spin glass models lacking inversion symmetry (i.e. spin-flip symmetry,
in the Ising case), the form of $q(x)$ in class 3) just described was found to
occur in RSB MF theory for a range of temperatures, and to correspond to a
“dynamically-frozen” phase in a dynamical approach Kirkpatrick and Thirumalai
(1987a, b, 1988a, 1988b). It was later argued that in a short-range spin
glass, this phase, which appeared to possess extensive complexity, would be
destroyed by entropic effects Kirkpatrick  _et al._ (1989); Bouchaud and
Biroli (2004), so that there would be no transition into it on lowering the
temperature (there was still supposed to be another transition at lower
temperature, which was termed a “random first-order transition” Kirkpatrick
_et al._ (1989); Bouchaud and Biroli (2004)). However, it does not seem clear
that such effects must completely destroy such a frozen phase; instead they
could leave a similar phase, still in class 3), now having a subextensive
complexity that obeys the bounds discussed above Höller and Read (2020); Read
(2022).

###  V.2 Application of RSB to the AW metastate

Now we turn to the use of RSB to study the AW metastate for the EA model Read
(2014). We consider a finite system in the hypercube $\Lambda_{L}$, with free
or periodic boundary conditions. A sub-hypercube $\Lambda_{R}$ will be viewed
as the “inner” region, while the sites and edges in $\Lambda_{L}-\Lambda_{R}$
constitute the “outer” region. We will consider copies of the system that all
experience the same disorder (bonds) in the inner region, but some of which
experience different (independently sampled) disorder in the outer region. An
expectation over the disorder in the outer region corresponds to an AW
metastate expectation using $\kappa_{J}$, while an expectation over disorder
in the inner region corresponds to the disorder expectation using $\nu$ (when
the limits $L\to\infty$ with $R$ fixed, followed by $R\to\infty$ are taken
after all the expectations). With this we can study different types of moments
of spin correlation functions, with metastate and disorder averages taken at
separate stages; disorder expectations are taken last, after all metastate
expectations have been done.

To carry out the expectations using the replica formalism, we introduce $l$
“groups” of $n_{k}$ copies or replicas of the system for the $k$th independent
sample of disorder in the outer region ($k=1$, …, $l$); replica indices run
from $1$ to $\sum_{k}n_{k}=n$, the total number of replicas in all groups; the
replica limit $n_{k}\to 0$ for all $k$ must be taken. The replicated theory
has permutation symmetry $S_{n_{1}}\times S_{n_{2}}\times\cdots S_{n_{l}}$,
broken from $S_{n}$ because of the different disorder experienced by replicas
in different groups in the outer region. However, locally in the inner region
the theory has full $S_{n}$ symmetry, so the effect is that of a symmetry-
breaking perturbation that is at infinity once the $L$, $R\to\infty$ limits
have been taken.

At MF level, applying the hierarchical ansatz, in the outer region the
replicas that experience distinct disorder should have small mutual overlaps,
but otherwise the structure of the order-parameter matrix $Q_{ab}$ (locally in
the outer region) within each group will be the same as usual (with $n_{k}$ in
place of $n$). In the inner region, the local $Q_{ab}$ matrix will have the
same form as when only one group of replicas is used. The presence of the
additional partition into blocks of sizes $n_{k}$ (which can be taken all
equal) which $\to 0$ does not affect the free-energy functional, and the
overlaps between replicas in distinct groups will here be $q_{0}$, reflecting
the explicit breaking of the $S_{n}$ symmetry by the outer region, with
subsequent further breaking as in the usual scheme.

As a simple example, we first consider the disorder and metastate expectation
$\left[\left[\langle\sigma_{\bf
x}\rangle_{\Gamma}^{2}\right]_{\kappa_{J}}\right]_{\nu}$ for ${\bf x}$ at say
the origin. For this, we choose distinct replicas $a$, $b$ in say the first
group, and because of the distinct stationary points of the MF theory that
differ by permutations in $S_{n_{1}}$, we should average over those choices;
we obtain

$\displaystyle\left[\left[\langle\sigma_{\bf
x}\rangle_{\Gamma}^{2}\right]_{\kappa_{J}}\right]_{\nu}$ $\displaystyle=$
$\displaystyle\lim_{n_{1}\to
0}\frac{1}{n_{1}(n_{1}-1)}\sum_{a,b=1}^{n_{1}}Q_{ab}$ (19) $\displaystyle=$
$\displaystyle\int_{0}^{1}dx\,q(x)=\overline{q}.$ (20)

On the other hand, if we take the metastate expectation before squaring, we
must take two replicas from distinct groups, and then averaging within each
group has no effect; we obtain

$\left[\left[\langle\sigma_{\bf
x}\rangle_{\Gamma}\right]^{2}_{\kappa_{J}}\right]_{\nu}=q_{0}=q(0),$ (21)

the result claimed earlier Read (2014).

Similar results are found for more complicated moments. If the metastate and
disorder expectations are taken together, the results are the same as in the
old RSB forms. But when some algebraic operation such as a square occurs
between the expectations, the results differ. In particular, the so-called
non-self-averaging (NSA) of $P_{J,\Gamma}(q)$ was found in RSB by using the
former type of expectations; it occurs whenever $q(x)$ is not constant on
$(0,1)$, so NSA occurs in RSB only if both the metastate and the typical Gibbs
states are non-trivial. When the present formalism is applied to moments of
metastate averages, no NSA is found for $P_{J}(q)$, in complete agreement with
the considerations of NS Newman and Stein (1996a, b, 1997). Thus within RSB,
fluctuations of $P_{J,\Gamma}(q)$ are associated with fluctuations of
(nontrivial) $\Gamma$s due to a nontrivial metastate, not the disorder
distribution $\nu$ on $J$ Newman and Stein (1997); Read (2014).

Next we turn to correlation functions. For the correlations in which the
metastate and disorder expectations are performed together, for example for
$\chi({\bf x},{\bf y})$, the results take the traditional form; the replicas
used are in the same group. On the other hand, for the MAS correlation $C_{\rm
MAS}({\bf x},{\bf y})$ one must choose two replicas from distinct groups. In
this case it is known in RSB field theory Dominicis  _et al._ (1998);
Dominicis and Giardina (2006) (though not applied there to the MAS) that the
correlation function in the spin glass phase of the EA model at zero magnetic
field has the asymptotic form

$C_{\rm MAS}({\bf x},{\bf y})\sim|{\bf x}-{\bf y}|^{-(d-4)}$ (22)

(within a constant factor), so that $\zeta=4$ Read (2014). The calculation
here, within a statistical field theory at lowest order without any loop
corrections, should be valid for dimensions $d>6$. It can be extended to
obtain the distribution $P_{J,\rho,W}(q)$ of window overlaps, with the
overlaps defined as in (11) but with $W$ kept fixed, for the MAS $\rho_{J}$.
In leading order as $W\to\infty$, again for the EA model in zero magnetic
field, the disorder expectation $\left[P_{J,\rho,W}(q)\right]_{\nu}$ of this
distribution is found to be a Gaussian of variance $\sim W^{-(d-4)}$ for
$d>6$, and further the distribution $P_{J,\rho,W}(q)$ is found to equal its
disorder expectation at leading order as $W\to\infty$ Read (2014). Thus in the
limit $W\to\infty$, $P_{J,\rho}(q)$ is found to be a $\delta$-function at
$q=q(0)$, and it self-averages, as expected from Refs. Newman and Stein
(1996a, b, 1997).

Some authors carried out a direct numerical study of the AW MAS along these
lines, in the three-dimensional nearest-neighbor EA model at moderate sizes,
and found evidence of a non-trivial metastate Billoire  _et al._ (2017).
Another study looked at a dynamical analog White and Fisher (2006) in a one-
dimensional power-law model, and found quantitative agreement with the
exponent that corresponds to $\zeta=4$ in the short-range case for $d>6$
Wittmann and Young (2016); that work has been extended further Jensen  _et
al._ (2021).

##  VI Conclusion

In spin glass theory there are outstanding controversies surrounding the
nature of the spin glass phase: is it a single (pair of) ordered state(s), as
suggested by the scaling-droplet theory, or are there many ordered states, as
suggested by replica-symmetry breaking (RSB) mean-field theory? The notion of
distinct ordered (or “pure”) states, and hence the question itself, is not
even well-defined except in an infinite system; hence some sort of infinite-
size limit must be taken. Due to the possible presence of chaotic size
dependence, the limit is not straightforward. The metastate constructions
provide a way out of this problem. A metastate is an additional layer of
structure in the theoretical framework: it is a probability distribution on
Gibbs states for given disorder, and each Gibbs state could be a mixture of
many pure states. A metastate contains information about how the states in
finite size vary with size or depend on the disorder far from the origin,
asymptotically at large sizes. Within this framework, a number of significant
results that constrain the allowed scenarios for a spin glass have been
obtained, and some of these results were reviewed in this contribution.

RSB theory initially takes the form of a mean-field theory for thermodynamic
properties. When applied to short-range spin glasses, it makes a number of
remarkable predictions about the pure states and their dependence on the
sample of disorder. It turns out, at least from a heuristic or physical point
of view, that RSB, when properly interpreted, has so far passed all the
consistency tests available from metastate theory, and at least for some cases
makes predictions about the metastate itself.

The ultimate answers to the controversies of short-range spin glasses are
still unknown. The nature of the problem raises severe difficulties for
traditional analytical methods of theoretical physics, while numerical methods
are subject to the limitations of finite size. It may be that rigorously-
proved theorems will play a definitive role in resolving the remaining
controversies and uncovering the true behavior of these fascinating systems.

###### Acknowledgements.

NR is grateful for the support of NSF grant no. DMR-1724923.

## References

  * Edwards and Anderson (1975) S. Edwards and P. W. Anderson, J. Phys. F 5, 965 (1975).
  * Sherrington and Kirkpatrick (1975) D. Sherrington and S. Kirkpatrick, Phys. Rev. Lett. 35, 1792 (1975).
  * Parisi (1979) G. Parisi, Phys. Rev. Lett. 43, 1754 (1979).
  * Anderson (1989) P. W. Anderson, Physics Today 42(9), 9 (1989).
  * Guerra and Toninelli (2002) F. Guerra and F. L. Toninelli, Commun. Math. Phys. 230, 71 (2002).
  * Guerra (2003) F. Guerra, Commun. Math. Phys. 233, 1 (2003).
  * Talagrand (rlin) M. Talagrand,  _Spin Glasses: A Challenge for Mathematicians_ (Springer-Verlag, 2003, Berlin).
  * Newman and Stein (1992) C. M. Newman and D. L. Stein, Phys. Rev. B 46, 973 (1992).
  * Newman and Stein (2022) C. M. Newman and D. L. Stein, Phys. Rev. E (2022), to appear, available at arXiv:2110.11229.
  * McMillan (1984) W. L. McMillan, J. Phys. C 17, 3179 (1984).
  * Bray and Moore (1985) A. J. Bray and M. A. Moore, Phys. Rev. B 31, 631 (1985).
  * Fisher and Huse (1986) D. S. Fisher and D. A. Huse, Phys. Rev. Lett. 56, 1601 (1986).
  * Fisher and Huse (1988) D. S. Fisher and D. A. Huse, Phys. Rev. B 38, 386 (1988).
  * Krzakala and Martin (2000) F. Krzakala and O. C. Martin, Phys. Rev. Lett. 85, 3013 (2000).
  * Palassini and Young (2000) M. Palassini and A. P. Young, Phys. Rev. Lett. 85, 3017 (2000).
  * Huse and Fisher (1987) D. A. Huse and D. S. Fisher, J. Phys. A 20, L997 (1987).
  * Parisi (1983) G. Parisi, Phys. Rev. Lett. 50, 1946 (1983).
  * Mézard  _et al._ (1984a) M. Mézard, G. Parisi, N. Sourlas, G. Toulouse, and M. Virasoro, Phys. Rev. Lett. 52, 1156 (1984a).
  * Mézard  _et al._ (1984b) M. Mézard, G. Parisi, N. Sourlas, G. Toulouse, and M. Virasoro, J. Phys. (Paris) 45, 843 (1984b).
  * Mézard  _et al._ (1987) M. Mézard, G. Parisi, and M. A. Virasoro, eds.,  _Spin Glass Theory and Beyond_ (World Scientific, Singapore, 1987).
  * Marinari  _et al._ (2000) E. Marinari, G. Parisi, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, and F. Zuliani, J. Stat. Phys. 98, 973 (2000).
  * Read (2014) N. Read, Phys. Rev. E 90, 032142 (2014).
  * Newman and Stein (1996a) C. M. Newman and D. L. Stein, Phys. Rev. Lett. 76, 515 (1996a).
  * Newman and Stein (1997) C. M. Newman and D. L. Stein, Phys. Rev. E 55, 5194 (1997).
  * Newman and Stein (1998a) C. M. Newman and D. L. Stein, Phys. Rev. E 57, 1356 (1998a).
  * Newman and Stein (2003) C. M. Newman and D. L. Stein, J. Phys.: Cond. Mat. 15, R1319 (2003).
  * Chung (2001) K. L. Chung,  _A First Course in Probability Theory, 3rd Ed._ (Academic, San Diego, CA, 2001).
  * Newman and Stein (1996b) C. M. Newman and D. L. Stein, Phys. Rev. Lett. 76, 4821 (1996b).
  * Newman and Stein (1998b) C. M. Newman and D. L. Stein, in  _Mathematics of Spin Glasses and Neural Networks_ , edited by A. Bovier and P. Picco (Birkhauser, Boston, 1998) pp. 243–287.
  * Newman (1997) C. M. Newman,  _Topics in Disordered Systems_ (Birkhauser, Basel, 1997).
  * Arguin  _et al._ (2016) L.-P. Arguin, C. M. Newman, D. L. Stein, and J. Wehr, J. Stat. Phys. 165, 1069 (2016).
  * Arguin  _et al._ (2014) L.-P. Arguin, C. M. Newman, D. L. Stein, and J. Wehr, J. Stat. Phys. 156, 221 (2014).
  * Aizenman and Wehr (1990) M. Aizenman and J. Wehr, Commun. Math. Phys. 130, 489 (1990).
  * Külske (1998a) C. Külske, in  _Mathematics of Spin Glasses and Neural Networks_ , edited by A. Bovier and P. Picco (Birkhauser, Boston, 1998) pp. 151–160.
  * Külske (1998b) C. Külske, J. Stat. Phys. 91, 155 (1998b).
  * van Enter and Schaap (2001) A. C. D. van Enter and H. G. Schaap, J. Phys. A: Gen. Phys. 35, 2581 (2001).
  * Newman and Stein (2000) C. M. Newman and D. L. Stein, Phys. Rev. Lett. 84, 3966 (2000).
  * Newman and Stein (2001) C. M. Newman and D. L. Stein, Commun. Math. Phys. 224, 205 (2001).
  * Arguin  _et al._ (2010) L.-P. Arguin, M. Damron, C. M. Newman, and D. L. Stein, Commun. Math. Phys. 300, 641 (2010).
  * Arguin  _et al._ (2019) L.-P. Arguin, C. M. Newman, and D. L. Stein, Commun. Math. Phys. 367, 1019 (2019).
  * Georgii (1988) H. O. Georgii,  _Gibbs Measures and Phase Transitions_ (de Gruyter, Berlin, 1988).
  * Simon (1993) B. Simon,  _The Statistical Mechanics of Lattice Gases_ (Princeton University Press, Princeton, NJ, 1993).
  * Young  _et al._ (1984) A. P. Young, A. J. Bray, and M. A. Moore, J. Phys. C 17, L149 (1984).
  * Newman and Stein (2006) C. M. Newman and D. L. Stein, in  _Spin Glass Theory_ , edited by E. Bolthausen and A. Bovier (Springer, Berlin, 2006) pp. 159–175.
  * Palmer (1982) R. G. Palmer, Adv. Phys. 31, 669 (1982).
  * van Enter and van Hemmen (1984) A. C. D. van Enter and J. L. van Hemmen, Phys. Rev. A 29, 355 (1984).
  * Cover and Thomas (2006) T. M. Cover and J. A. Thomas,  _Elements of Information Theory_ (John Wiley & Sons, Hoboken, NJ, 2006) Chs. 2 and 8\.
  * Höller and Read (2020) J. Höller and N. Read, Phys. Rev. E 101, 042114 (2020).
  * Read (2022) N. Read, Phys. Rev. E (2022), to appear, available at arXiv:2108.11849.
  * Gross and Mézard (1984) D. Gross and M. Mézard, Nucl. Phys. B 240, 431 (1984).
  * Newman and Stein (1994) C. M. Newman and D. L. Stein, Phys. Rev. Lett. 72, 2286 (1994).
  * White and Fisher (2006) O. L. White and D. S. Fisher, Phys. Rev. Lett. 96, 137204 (2006).
  * Jackson and Read (2010) T. S. Jackson and N. Read, Phys. Rev. E 81, 021130 (2010).
  * (54) C. M. Newman, N. Read, and D. L. Stein, in preparation.
  * (55) J. Höller and N. Read, in preparation.
  * Kirkpatrick and Thirumalai (1987a) T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. Lett. 58, 2091 (1987a).
  * Kirkpatrick and Thirumalai (1987b) T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 36, 5388 (1987b).
  * Kirkpatrick and Thirumalai (1988a) T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 37, 5342 (1988a).
  * Kirkpatrick and Thirumalai (1988b) T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 38, 4881 (1988b).
  * Kirkpatrick  _et al._ (1989) T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys. Rev. A 40, 1045 (1989).
  * Bouchaud and Biroli (2004) J.-P. Bouchaud and G. Biroli, J. Chem. Phys. 121, 7347 (2004).
  * Dominicis  _et al._ (1998) C. D. Dominicis, I. Kondor, and T. Temesvari, in  _Spin Glasses and Random Fields_ , edited by A. P. Young (World Scientific, Singapore, 1998).
  * Dominicis and Giardina (2006) C. D. Dominicis and I. Giardina, in  _Random Fields and Spin Glasses_ , edited by A. P. Young (Cambridge University Press, Cambridge, 2006).
  * Billoire  _et al._ (2017) A. Billoire, L. A. Fernandez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, G. Parisi, F. Ricci-Tersenghi, and J. Ruiz-Lorenzo, Phys. Rev. Lett. 119, 037203 (2017).
  * Wittmann and Young (2016) M. Wittmann and A. P. Young, J. Stat. Mech. , 013301 (2016).
  * Jensen  _et al._ (2021) S. Jensen, N. Read, and A. P. Young, Phys. Rev. E 104, 034105 (2021).


Towards MM1, the averaged H13CN $J=1-0$ line is suspected to contain two
narrow velocity components, a strong and a weak one. Besides, H13CN seems to
be contaminated by strong and dense outflow from MM1. Therefore, we simply
consider a fitting model with one component for narrow and strong HFS, and one
Gaussian component for extended wings from outflow.

To fit the hyperfine structure of H13CN $J=1-0$, we make the following
assumptions: 1) all the HFS lines share the same excitation temperature; 2)
the opacity as a function of frequency (velocity) has a Gaussian profile; 3)
all the lines share the same linewidth. Since we have three individual HFS
lines (Table 8), the opacity of the $i$th component is written as,

$\tau_{i}(v)=\tau_{i}\cdot\exp\left[-4\ln 2\left(\frac{v-v_{0,i}}{\Delta
v}\right)^{2}\right],$ (19)

where $\Delta v$ is the uniform FWHM of all components. The central velocity
of component $i$ is $v_{0,i}=v_{i}+v_{\mathrm{LSR}}$, where $v_{\mathrm{LSR}}$
is the velocity of the reference component (i.e. the one with $v_{i}=0$) at
the Local Standard of Rest (LSR).

The opacity of the multiplet writes:

$\tau(v)=\tau_{\mathrm{tot}}\sum\limits_{i=1}^{N}r_{i}\times\exp\left[-4\ln
2\left(\frac{v-v_{\mathrm{LSR}}-v_{i}}{\Delta v}\right)^{2}\right]$ (20)

where $\tau_{\mathrm{tot}}$ is the sum of the optical depths of all hyperfine
lines and $r_{i}$ is the normalized relative intensity of the individual
hyperfine line in the optically thin case under the condition of
thermodynamical equilibrium (TE).

Given the opacity $\tau(v)$ and the amplitude $\mathcal{A}$, the brightness
temperature is obtained from:

$T_{b}(v)=\mathcal{A}(1-e^{-\tau(v)})$ (21)

In the optical thin regime, the brightness temperature $T_{b}(v)$ transforms
into:

$T_{b}(v)=\mathcal{A}\tau(v)=\mathcal{A}\tau_{\mathrm{tot}}\sum\limits_{i=1}^{N}r_{i}\times\exp\left[-4\ln
2\left(\frac{v-v_{\mathrm{LSR}}-v_{i}}{\Delta v}\right)^{2}\right]$ (22)

Each hyperfine component is considered as a single Gaussian profile at its
centroid and neighborhood ($\sim 3\sigma$). In other words, at the $\sim
3\sigma$ range of the hyperfine component $i$, the brightness temperature can
be written as:

$T_{b}(v)|_{-3\sigma\,\leq\,v_{0i}\,\leq\,3\sigma}=\mathcal{A}\tau_{\mathrm{tot}}r_{i}\times\exp\left[-4\ln
2\left(\frac{v-v_{\mathrm{LSR}}-v_{i}}{\Delta v}\right)^{2}\right]$ (23)

and it is not possible to determine $\mathcal{A}$ and $\tau_{\mathrm{tot}}$
degenerately. To stabilize the numerical calculation of HFS fitting, we
therefore define the product of the amplitude $\mathcal{A}$ and the total
optical depth $\tau_{\mathrm{tot}}$ as a variable
$T_{\mathrm{ex}}\tau=\mathcal{A}\tau_{\mathrm{tot}}$ (Levshakov et al., 2013;
Liu et al., 2021). The other three variables are respectively
$V_{\mathrm{LSR}}$, $\Delta V$, and $\tau_{\mathrm{tot}}$. The HFS fitting is
performed by python package lmfit with “leastsq” method to minimize the
$\chi^{2}$. The result is shown in Figure 13.

Figure 13: Result of fitting H13CN $J=1-0$. The fitting model contains one
hyperfine structure component and one Gaussian component. The former has
narrow linewidth and strong emission. The latter has extended wings. The
fitted parameters are listed on the left. Since the optical depth of H13CN
$J=1-0$ is small, the amplitude (the function of excitation temperature
$T_{\mathrm{ex}}$) and optical depth are coupled. We use the parameter
$T_{\mathrm{ex}}\tau$ to represent their product.

##  Appendix F Optical depth of H13CO+ $J=1-0$ line

To estimate the optical depth of H13CO+ $J=1-0$, we assume local thermodynamic
equilibrium (LTE), and therefore the brightness temperature $T_{b}$ can be
derived from

$T_{b}=f\left[J_{\nu}(T_{\mathrm{ex}})-J_{\nu}(T_{\mathrm{bg}})\right](1-e^{-\tau}),$
(24)

where $J_{\nu}(T)=\frac{h\nu}{k_{\mathrm{B}}}\frac{1}{\exp{(h\nu/kT)}\,-1}$.
The beam filling factor $f$ is taken to be 1 since ALMA can well resolve the
structure. The cosmic background radiation $T_{\mathrm{bg}}$ is 2.73 K.

HCO+ $J=1-0$ is taken to be optically thick with $\tau\gg 1$. We simplify it
as $T_{b}=J_{\nu}(T_{\mathrm{ex}})-J_{\nu}(T_{\mathrm{bg}})$, with which we
estimate the excitation temperature $T_{\mathrm{ex}}$ for HCO+ $J=1-0$. We
simply assume that the two molecules share the same $T_{\mathrm{ex}}$ and
return to Equation 24. For $\tau\ll 1$, the original equation becomes
$T_{b}=\tau\left[J_{\nu}(T_{\mathrm{ex}})-J_{\nu}(T_{\mathrm{bg}})\right]$.
Substituting the parameters of H13CO+ $J=1-0$, $\tau$(H13CO+) can be derived.
As seen from Figure 14, the optical depth $\tau$(H13CO+) is universally thin,
even for the densest part near the MM1.

We also test the optically thick assumption for HCO+ $J=1-0$. Based on the
radiation transfer equation 24, we have,

$\begin{split}\frac{T_{b,\mathrm{HCO}^{+}}(v)}{T_{b,\mathrm{H}^{13}\mathrm{CO}^{+}}(v)}&=\frac{1-\exp(-\tau_{v,12})}{1-\exp(-\tau_{v,13})}\\\
&>\frac{1-\exp(-\tau_{v,12})}{\tau_{v,13}}=\frac{1-\exp(-\tau_{v,12})}{\tau_{v,12}}\frac{X_{12}}{X_{13}}\end{split}$
(25)

where $\tau_{v,12}$ and $\tau_{v,13}$ are the optical depth at $v$ of the HCO+
$J=1-0$ and H13CO+ $J=1-0$, and $X_{12}/X_{13}$ = [HCO+]/[H13CO+] is the
abundance ratio of HCO+ to H13CO+, which is estimated to be 15–20 from Mopra
observations and 1D non-LTE RATRAN radiative transfer modeling (Peretto et
al., 2013).

Since the linewidth of HCO+ $J=1-0$ is always larger than that of H13CO+
$J=1-0$, we only consider the peak brightness temperature at
$v=v_{\mathrm{peak}}$. Among the whole map with solid data (both two lines
should have SNR>5),
$T_{b,\mathrm{HCO}^{+}}(v=v_{\mathrm{peak}})/T_{b,\mathrm{H}^{13}\mathrm{CO}^{+}}(v=v_{\mathrm{peak}})$
is 2.7 on average and 9 at maximum. By substituting these into equation 25, we
estimate the lower limit of $\tau_{12}$ at $v=v_{\mathrm{peak}}$ to be 7.4 on
average, and only is as low as $\sim$2 in extreme cases. In other words, the
optically thick assumption, $\tau_{12}$ > 1, is always satisfied.

Figure 14: Demonstration of the optical thinness of H13CO+ $J=1-0$. Top:
spatial distribution of optical depth $\tau$(H13CO+) overlaid with ALMA Band-3
continuum contours. The white contours of 3 mm continuum emission follows the
power-law levels of [1.0, 3.3, 7.6, 14.2, 23.4, 35.1, 49.6, 67.0] mJy beam-1.
Two massive dense cores MM1 and MM2 are shown. The pixels with SNR < 5 are
masked. The ALMA beam is shown at the lower left. Bottom: distribution of
$\tau$(H13CO+) on the $T_{b}$-$T_{\mathrm{ex}}$ plane. $\tau=$0.3, 0.5, 1, 2,
and 4 are shown with black dashed lines. The pixels within the MM1 mask
($F_{\mathrm{cont,3mm}}\geq 1.0$ mJy) are shown with gray hollow circles.

##  Appendix G SCOUSEPY decomposition of the ALMA H13CO+ $J=1-0$ map

The first release version of scousepy is scouse (short for Semi-automated
multi-COmponent Universal Spectral-line fitting Engine Henshaw et al., 2016).
scouse is an IDL package that fits complex spectral line data in a robust,
systematic and efficient way, by manually fitting the spatially-averaged
spectra into Gaussian components and then using them as input to the
subsequent fitting on each individual spectrum. Henshaw et al. (2019)
developed Python implementation scousepy999SCOUSEPY is publicly available:
https://github.com/jdhenshaw/scousepy. Alternatively, the original IDL
implementation can be downloaded here: https://github.com/jdhenshaw/scouse.
and its latest version scousepy v2.0.0 improves the workflow and operation
interface. The procedure of scousepy is discussed in detail by Henshaw et al.
(2016), but we highlight the key points here.

First, we identify both the spatial and spectral coverage over which it will
perform the fitting. We calculate the pixelwise noise from the baseline
without H13CO+ emission and mask the spectra whose peak flux is below
$5\sigma$. The reason to do this is that the primary-beam corrected ALMA data
has radial variation of noise. By doing so, the unmasked data has a uniform
$\mathrm{SNR}>5$ for further spectral line fitting. We also cut the velocity
range from -55 to -34 km s-1, both covering all the emission from H13CO+
$J=1-0$ and releasing the burden of loading data. Second, we set the Spectral
Averaging Areas (SAAs) to be $8\times 8$ pixel2 and scousepy then breaks up
the map (with 8873 pixels) into 555 small areas as shown in Figure 15. The
spatially averaged spectra extracted from each SAA are then manually fitted.
The spectra is assumed to be composed of a series of Gaussian profiles, as
justified in Section 4.1.1. Best-fitting solutions to the SAAs are then
supplied to the fully automated fitting procedure that targets all of the
individual spectra contained within each region. This process is controlled by
a number of tolerance levels. Last, we check the fitting result and the
quality due to the reduced $\chi^{2}$. Some minor manual adjustment can be
done here.

Figure 15: The breaking-up stage of scousepy where the HCO+ $J=1-0$ map is
broken into 555 SAAs.

##  Appendix H Clustering of decomposed Gaussian components by ACORNS

Henshaw et al. (2019) developed a Python written algorithm, acorns
(Agglomerative Clustering for ORganising Nested Structures)101010ACORNS is
publicly available at: https://github.com/jdhenshaw/acorns.. acorns is based
on a technique known as hierarchical agglomerative clustering whose primary
function is to generate a hierarchical system of clusters within discrete
data. Refer to Henshaw et al. (2019) for the principle and details of the
algorithm. In the following sections, we briefly describe how acorns
characterizes the velocity structure of “the Heart”.

We perform the acorns decomposition only on the most robust spectral velocity
components extracted by scousepy. We define “robust” as all velocity
components whose peak flux density is greater than 5<$\sigma_{\mathrm{rms}}$>
where <$\sigma_{\mathrm{rms}}$> is the averaged rms noise among the all
pixels, which is $\sim$9.5 mJy beam-1. The selected data (15879 components out
of 17143) constitute 92.6% of the total data set extracted by scousepy.

For the clustering, we set the minimum radius of a cluster to be 2.4 arcsec,
which is equal to the ALMA synthesized beam. This is to ensure that all
identified clusters are spatially resolved. The consequence is that 545 data
points are not assigned to clusters and isolated. Considering clustering in
p–p–v space, for two data points to be classified as “linked”, we specify that
the Euclidean distance between the points and the absolute difference in both
their measured centroid velocity and velocity dispersion can be no greater
than 2.4 arcsec and 0.21 km s-1, respectively111111The velocity difference
0.21 km/s is equal to the velocity resolution of the H13CO+ $J=1-0$
observations..

In the case of “the Heart”, acorns builds a forest with a total of 12 trees,
using 98 per cent of the selected “robust” data. As seen from Figure 16. the
forest consist of four hierarchical and eight non-hierarchical
trees121212Henshaw et al. (2019) expands on the nomenclature used in
dendrograms. The whole tree system is called “forest”, itself containing
numerous trees. Each tree may or may not then be further subdivided into
branches or leaves. Trees with no substructure are also classified as leaves..
The forest is dominated by four hierarchical trees: #0, #1, #3, and #7, which
in total contain 75.7% of all robust data. Given the enormity of the data set,
we call them as the major streams A, B, C, and D, respectively. “Stream”
refers to the gas flowing (see Section 4.2 for its physical nature). The
remaining data are classified as isolated points and are not included in the
further analyses.

Figure 17 shows the integrated flux maps of the four major streams.

Figure 16: The visualization of acorns outputs of “the Heart”. Top: acorns
forest, the hierarchical structure of the data points from acorns. In total,
12 trees are found and each is plotted in different color. Four are
hierarchical (#0, #1, #3, and #7) and eight are non-hierarchical. The four
hierarchical trees are named major streams A, B, C, and D, respectively.
Bottom: p–p–v 3D plot of the forest. The colors are the same as the top panel.
Figure 17: The integrated flux maps of the four major streams are normalized
to linear span from 0.0 to 0.4 Jy beam-1 km s-1 (0.0 to 0.8 Jy beam-1 km s-1
for the stream D) in colorscale. The integrated flux at each pixel is given by
the integration of decomposed Gaussian component(s). The black solid contours
depict the 5$\sigma$ level of H13CO+ $J=1-0$. The six groups of dashed gray
lines are the filaments identified from the Spitzer 8 $\mu$m extinction map
(the same as in the middle panel of Figure 1). The black bold dashed lines
mark the filaments which are assumed to be responsible for the streams in each
panel. Streams including two major ones (A and B) and four minor ones (C1, C2,
D1, and D2) are marked as red bold dashed vectors. The massive dense cores MM1
and MM2 together with their boundaries ($F_{\mathrm{cont,3mm}}=1.0$ mJy
beam-1) are shown in purple. The central circle marks how we visualize the 3D
data cube: blue axis is the collapsing axis, where integrated fluxes of pixels
are projected on the plane spanned by two black axes.