replicas with the same quenched disorder. Provided that the sample is kept in
the same experimental conditions for all the data acquisition time, then real
replicas can be associated to the different shots, each one thermalized into
different equilibrium states characterized by a specific spectral profile of
activated modes and sharp peaks. Thermalization is guaranteed by the fact that
during a single pulse of the external pumping, several stimulated emission
phenomena take place for each mode frequency ensuring a long enough mode
dynamics. To be precise, this would only be a partial thermalization, since
for a disordered system in the ergodicity broken phase, a complete
thermalization would require the system to visit all possible states. Given
the development of the theory and the current interpretation of experiments,
we cannot say how many equilibrium states a random laser actually visits for
each spectral shot, but we can quite safely say that the system has reached
the static transition predicted by mean-field theory. This is not only a
consequence of fast light mode dynamics, but also of a number of modes which
is small compared to usual thermodynamics degrees of freedom (e.g. $\sim
10^{23}$) and of the fact that these modes are highly connected as in the
dense interaction network of a mean-field model.

Let us first define the fluctuation of the intensity $I_{k}^{\alpha}$ of the
resonance at the frequency $\omega_{k}$ in a single spectrum $\alpha$ with
respect to the spectral intensity at that frequency averaged over all $N_{s}$
acquired spectra as

$\Delta_{k}^{\alpha}=I_{k}^{\alpha}-\frac{1}{N_{s}}\sum_{\gamma=1}^{N_{s}}I_{k}^{\gamma}.$
(2.72)

Each spectrum represents the realization of a replica. The experimental IFO
measured between two real replicas can be represented by the following matrix

$\mathcal{C}_{\alpha\beta}^{\text{exp}}=\frac{\sum_{k}\Delta_{k}^{\alpha}\Delta_{k}^{\beta}}{\sqrt{\vphantom{\big{|}}\sum_{k}(\Delta_{k}^{\alpha})^{2}}\sqrt{\sum_{k}(\Delta_{k}^{\beta})^{2}}}$
(2.73)

defined in the interval $[-1,1]$, where $I_{k}^{\alpha}$ denotes the intensity
of the mode $k$ in the spectrum corresponding to the replica (shot) $\alpha$,
with $\alpha=1,\dots,N_{s}$. Since thermalization is assumed, the experimental
value $I_{k}^{\alpha}$ can be thought as the equilibrium average of the
intensity, i.e.
$I_{k}^{\alpha}\equiv\frac{1}{\mathcal{T}}\int_{t_{0}}^{t_{0}+\mathcal{T}}\text{d}t|a_{k}^{\alpha}(t)|^{2}$,
where $\mathcal{T}-t_{0}$ is the time interval corresponding to random laser
lifetime, slightly longer than the pumping pulse. The overlap is defined
between intensity fluctuations rather than directly between intensities, in
order to exclude the effect of amplified spontaneous emission on the
measurements. Fluctuations are taken with respect to the intensity averaged
over many different replicas. From the $N_{s}$ measured spectra one can
extract $N_{s}(N_{s}-1)/2$ values of the IFO and determining their
distribution by building the histogram
$P(\mathcal{C})=\sum_{\alpha<\beta}\delta(\mathcal{C}-\mathcal{C}_{\alpha\beta})$.

Figure 2.3: Six histograms representing the probability distribution function
of the experimental intensity fluctuation overlap (2.73) built from the data
obtained by many single shot spectra collected at six increasing values of the
pumping rate (Reprinted from [Gho+14]). Notice that $P(\mathcal{C})$ gradually
changes from a Gaussian-like distribution centered in zero for low values of
the $\mathcal{P}$ to a bimodal distribution. The symbols represent the region
of the phase diagram in Fig. 2.2 to which these experimental $P(\mathcal{C})$
might correspond. The correspondence is clarified by the analytical study of
the IFO distribution function, whose results are presented in Fig. 2.4.

We report in Fig. 2.3 the results obtained in Ref. [Gho+14], concerning the
measurement of the IFO distribution. At low pumping rate, $P(\mathcal{C})$
appears as a Gaussian-like distribution centered in $\mathcal{C}=0$. Then, for
increasing values of $\mathcal{P}$, the distribution develops a nontrivial
structure with three distinguished peaks, one in $\mathcal{C}=0$ and two
symmetric side-peaks, and a continuous part between them. Eventually,
$P(\mathcal{C})$ reduces to a double-peaked distribution for high values of
$\mathcal{P}$: in this last case, as the pumping is varied, $\mathcal{C}$ can
in principle take all possible values in the interval $[-1,1]$, while for a
given value of $\mathcal{P}$ the position of the peaks is fixed. This kind of
behavior resembles the one of the Parisi overlap distribution function in the
replica symmetry broken phase of the model.

Figure 2.4: Probability distribution of the intensity fluctuation overlap
$P(\mathcal{C})$ for competing linear and nonlinear couplings
($\alpha=\alpha_{0}=0.4$), and in the fieldless case $R_{J}=1.1$ (reprinted
from [Ant16]). The pumping rate $\mathcal{P}$ increases from left to right and
from top to bottom, along the dashed line in Fig. 2.2. At $\mathcal{P}=1.3$
the distribution is a Dirac delta centered in zero. As soon as the pumping
rate exceeds the lasing transition threshold, all the RSB regime are
displayed: FRSB ($\mathcal{P}=2.07$), 1-FRSB ($\mathcal{P}=2.23$) and 1RSB
($\mathcal{P}=3.12,4.45,7.03$).

In order to build a precise correspondence at least in the mean-field fully-
connected model, one has to define a quantity, depending only on intensity
fluctuations to be analytically related to the Parisi overlap. In Ref.
[ACL15a], the IFO is expressed (in absolute value) by the following matrix

$\displaystyle\mathcal{C}_{\alpha\beta}=\frac{1}{8\epsilon^{2}N}\sum_{k=1}^{N}\left[\langle|a_{k}^{\alpha}|^{2}|a_{k}^{\beta}|^{2}\rangle-\langle|a_{k}^{\alpha}|^{2}\rangle\langle|a_{k}^{\beta}|^{2}\rangle\right]$
(2.74)

defined in $[0,1]$, where the average is taken with respect to the Gibbs-
Boltzmann measure of the spherical (2+4)-phasor model. Clearly this
distribution depends on the realization of the disordered couplings $J$ and
has to be averaged $P(\mathcal{C})=\overline{P_{J}(\mathcal{C})}$.

It is worth stressing that the definition of emission spectra in statistical
mechanical models is only possible in a model which takes into account both
phases and intensities, and hence the introduction of the IFO distribution
could not be possible within the phase-only approach originally developed in
Refs. [Ang+06, Ang+06a, Leu+09, CL11]. While the role of the phases is
essential for reproducing the phase transition phenomenology of mode-locking,
the role of the intensities is of crucial importance for bridging with the
experiments, where we do not have access to the phases.

The crucial result obtained in Ref. [ACL15a] is that the IFO matrix defined in
Eq. (2.74) can be expressed in terms of the overlap matrices $\mathcal{Q}$ and
$\mathcal{R}$ as

$\displaystyle\mathcal{C}_{\alpha\beta}=\mathcal{Q}_{\alpha\beta}^{2}-\frac{m^{4}}{4}~{}~{}~{}~{}~{}~{}~{}~{}~{}a\neq
b$ (2.75)
$\displaystyle\mathcal{C}_{\alpha\alpha}=\frac{1+\mathcal{R}_{\alpha\alpha}^{2}}{2}-\frac{m^{4}}{4},$
(2.76)

element by element, whatever the structure of $\mathcal{Q}$ and $\mathcal{R}$.
This result reveals that if a RSB structure is present at the level of the
configuration overlap, the same holds also for the IFO: in other terms the
structure of the state organization is the same whether we look at the
configurations of the modes or at their intensity spectra.

In Fig. 2.4 six different plots of the analytical IFO probability distribution
are displayed for increasing values of the pumping rate $\mathcal{P}$ along
the dashed line at $\alpha=0.4$ in Fig. 2.2, which goes through all kinds of
RSB phases. At low pumping rate the distribution is a Dirac delta centered in
zero, meaning that no correlations are present among the intensities and the
modes are independent and non-interacting (CW phase); increasing the pumping
the mode coupling becomes relevant and, accordingly, the overlap distribution
function is nontrivial, since the system enters a phase where the modes are
highly frustrated by disorder. First, $P(\mathcal{C})$ develops a small
continuous part around the central peak in $\mathcal{C}=0$, denoting the
typical continuous FRSB shape (second panel); then, besides the continuous
part, also symmetric side Dirac deltas emerge, corresponding to the 1-FRSB
phase (third panel). Eventually, the analytical $P(\mathcal{C})$ looses the
continuous part and becomes a linear combination of Dirac deltas, which is the
usual 1RSB structure. The resemblance with the experiments, though only
qualitative, is quite remarkable.

Despite this similarity, it should be noted that the correspondence between
theory and experiments is still under construction and many criticisms, both
on the theoretical and experimental sides, can be raised with the perspective
of improving it. The most obvious criticism is that, for now, no experimental
data corresponding to different samples are yet available for averaging over
the disorder. Due to the non-self-averageness of overlap probability
distribution functions, the average over disorder is essential to observe the
typical behavior of these observables and compare it with the theoretical
predictions. Furthermore, a major problem lies in the fact that the emission
of a random laser occurs in every direction, while the acquisition of spectra
is not performed in the whole solid angle (at least not in Ref. [Gho+14]).
This, coupled with spectral resolution issues, casts doubt on whether the
detected spectra correspond to all the laser modes oscillating in the sample.
Moreover, experimental data inevitably contain a part of dynamical relaxation
to equilibrium, which should be taken into account by the theory in order to
improve the comparison. Clearly, the analytical results obtained for the IFO
distribution in Ref. [ACL15a] and reported in Fig. 2.4 are purely at
equilibrium.

An additional issue that has to be addressed on the theoretical side is going
beyond the narrow bandwidth limit, which includes the FMC in a trivial way.
This is precisely the big goal of this Thesis work. In particular, we aim to
develop a numerical tool to simulate the model diluted with the FMC both out-
of and at equilibrium. Hopefully, this will also provide useful insight on the
diluted model, in view of its analytical solution.

##  Part I Numerical Simulations

###  Chapter 3 Mixed-Order Glass Transition in Random Lasers

The main goal of this work is to go beyond the fully-connected solution of the
spin-glass model for random laser presented in the previous chapter. By this,
we mean to release the narrow bandwidth approximation and include in a
nontrivial way the mode coupling selection induced by the FMC. Out of the
narrow bandwidth approximation, the interaction network of the glassy random
laser can no longer be considered fully-connected, but has to be diluted by
removing all the bonds that do not match the condition on the frequencies. The
FMC is of key importance for the study of mode-locking and for the
reproduction of real random laser spectra: for this reason its inclusion is
essential to bridge with the experiments. However, when the solution of the
mode-locked model is approached analytically, subtle technical difficulties
emerge, which require the development of new techniques with respect to
standard mean-field methods for disordered systems. The analytical approach
will be developed in the second part of this work. Here, we resort to
numerical simulations to get useful insights on the mode-locked model.

A first step towards the inclusion of the FMC has been taken on the ordered
version of the model [AIL15a, AIL15]. In this case, strong deviations from the
fully-connected behavior have been put in evidence. In particular, the mode-
locked low temperature (high pumping) phase exhibits lack of global order due
to the onset of _phase waves_ 111Phase waves is an evocative term reminding of
_spin waves_ in pairwise spin models with $O(2)$ global symmetry, such as the
XY model. In these models, when the symmetry is spontaneously broken, the
global magnetization is reduced by the onset of collective excitations
analogous to the Goldstone bosons in quantum field theory. In dimensions
$d\leq 2$ this phenomenon leads to lack of global magnetic order, as stated by
the Mermin-Wagner theorem [MW66]; however, in the special case of $d=2$,
topological transitions of the Kosterlitz–Thouless type [KT73] may be allowed.
In the case of multimode lasers, however, we are dealing with dense models, so
the comparison can not be pushed too far. (see also Ref. [ML15]) produced by
the tendency to align of modes which are close in frequency, due to the FMC.
In particular the phases $\phi_{k}$ of the modes are not all equal as in the
fully-connected case, but satisfy a linear relation with their frequencies
$\omega_{k}$, which for a linear comb can be written as
$\phi_{k}=\phi_{0}+k\Delta$, where $\Delta$, approximately independent from
$k$, is the slope of the phase wave and is a configuration-dependent quantity
[Ant16]. As a result, the coherency of the laser is not trivial as in the
narrow bandwidth limit, where the output consists of a train of almost
perfectly delta-like (unchirped) pulses, but there is a phase delay in the
emitted pulses, which depends on the slope $\Delta$, and hence on the
configuration. Therefore, in the magnetic analogy, when considering the
thermal average, the model has vanishing magnetization. This effect has not
been found in a model with ordered couplings and a random dilution of the same
order of the FMC: in this case, the physical properties of the system are
coherent with the fully-connected solution. Similarly, phase waves are not
expected in the presence of quenched disordered couplings.

After that, equilibrium numerical simulations of the spin-glass mode-locked
model have been performed in Ref. [GAL20], where evidence of a mixed-order
phase transition has been found and put in connection with an equipartition-
breaking transition at the same critical temperature. The common root of the
two transitions can be traced back to the same underlying phenomenon: the
breaking of ergodicity. In the present chapter and in the following one, we
focus on the mixed-order phase transition, while the analysis of the
equipartition-breaking transition will be deepened in Chap. 5. With respect to
the mean-field picture described in the previous chapter, we are particularly
interested in checking what of the RFOT scenario remains in the diluted mode-
locked model, which is much closer to real random lasers than the fully-
connected one.

In the following, first, we present the simulated spin-glass model, which is a
slightly simplified version of the mode-locked spherical (2+4)-phasor model.
Particular attention is devoted to the role played by the FMC in affecting the
topology of the model interaction graph. The numerical technique is explained
in detail, by presenting the Exchange Monte Carlo algorithm implemented to
shorten the thermalization time to equilibrium at low temperature. Moreover,
the simulated model presents the additional problem of being defined on an
interaction graph, which, though diluted, is still very dense. To address this
problem, parallel computing on graphic processing units has been adopted. The
results pointing towards the presence of a static glass transition are
collected and their problematic nature is discussed. In particular, the
unexpected scaling of the critical region found in [GAL20] motivates the need
for a new campaign of numerical simulations of the model aimed at collecting
data less affected by finite-size effects.

####  3.1 The Mode-Locked 4-Phasor Model

The simulated model is described by the following Hamiltonian

$\displaystyle\mathcal{H}[\bm{a}]$
$\displaystyle=-\sum_{\bm{k}|\text{FMC}(\bm{k})}J_{k_{1}k_{2}k_{3}k_{4}}\overline{a}_{k_{1}}a_{k_{2}}\overline{a}_{k_{3}}a_{k_{4}}+\text{c.c.}$
$\displaystyle=-\sum_{\bm{k}|\text{FMC}(\bm{k})}J_{k_{1}k_{2}k_{3}k_{4}}A_{k_{1}}A_{k_{2}}A_{k_{3}}A_{k_{4}}\cos(\phi_{k_{1}}-\phi_{k_{2}}+\phi_{k_{3}}-\phi_{k_{4}}),$
(3.1)

where $\bm{a}=\\{a_{1},...,a_{N}\\}$ is a $N$-dimensional complex vector of
electromagnetic field mode amplitudes. In the second expression, $A_{k}$ and
$\phi_{k}$ represents respectively the modulus and the phase of the mode
amplitude $a_{k}=A_{k}e^{i\phi_{k}}$ and a factor 2 has been absorbed in the
definition of the random couplings. Configurations are constrained to the
complex hypersphere of radius $\sqrt{\epsilon N}$, where
$\epsilon=\mathcal{E}/N$ measures the average optical power per mode available
in the system. The quenched disordered coupling constants
$J_{\bm{k}}=J_{k_{1}k_{2}k_{3}k_{4}}$ are independently drawn from a zero-mean
Gaussian distribution

$\displaystyle
P(J_{\bm{k}})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left\\{\frac{J_{\bm{k}}^{2}}{2\sigma^{2}}\right\\},$
(3.2)

with variance $\sigma^{2}=\overline{J_{\bm{k}}^{2}}=1/N^{2}$, ensuring the
extensivity of the energy. The scaling of the coupling distribution variance
takes into account the dilution order of the interaction graph, which is
determined by the condition

$\displaystyle\text{FMC}(\bm{k}):|\omega_{k_{1}}-\omega_{k_{2}}+\omega_{k_{3}}-\omega_{k_{4}}|\lesssim\gamma,$
(3.3)

where $\omega_{k}$ are the frequencies of the modes and $\gamma$ denotes their
typical linewidth. We refer to this simplified version of the general
(2+4)-phasor model discussed in the previous chapter as mode-locked (ML)
4-phasor model. With respect to the general model and in terms of the photonic
parameters introduced in the previous chapter, see eqs. (2.44) and (2.45), we
are here working in the limits $R_{J}\rightarrow\infty$ and $\alpha=1$. The
motivation for considering this model is that we are mainly interested in the
study of the non-linear term of the Hamiltonian defined in Eq. (2.40), which
is the most relevant one for reproducing the phenomenology of optical waves in
disordered media near the lasing transition. Indeed, the behaviour of
multimode optical systems in this regime is generally believed to be dominated
by non-linear mode interactions, see e.g. Refs. [Hau00, And+11].

The simplest choice for the frequency distribution is to consider a linear
comb as in the case of standard lasers, see Eq. (2.2). In this case the FMC
(3.3) can be mapped into a relation among the indices of the interaction
graph:

$\displaystyle|k_{1}-k_{2}+k_{3}-k_{4}|=0.$ (3.4)

More realistic dilution rules based on random frequency distributions will be
considered in future works in order to improve the modeling of real random
lasers. Besides being the simplest possible choice, the frequency comb
distribution is compatible with the strong-cavity approximation [CL11], which
amounts to neglect the off-diagonal elements of the linear interaction term in
the Hamiltonian of the 2+4 model. In fact, the 2-body FMC, which generally
looks like $|\omega_{k_{1}}-\omega_{k_{2}}|\lesssim\gamma$, admits off-
diagonal terms only in the case of modes whose frequencies differ less than
the threshold fixed by $\gamma$. In principle, modes of this kind exist in
random lasers [And+11]. However, in the case of a high-finesse linear comb,
these modes are excluded: the condition for mode selection reduces to
$|k_{1}-k_{2}|=0$, which leaves only the diagonal terms.

Furthermore, by assuming a flat gain curve, i.e. by taking $J_{kk}=g$ with
constant $g$, the linear part of the interactions becomes an additive
constant, which is only responsible for a shift of the energy and can be
neglect it. This follows as a consequence of the spherical constraint. The
assumption of a flat gain curve is also compatible with the regime we aim to
explore through numerical simulations: as shown in [AIL15] for the case of
standard multimode lasers, the inclusion of a more complex gain profile only
affects the fluorescence regime, while the transition and the lasing regime
are stable under perturbations of the gain.

The effective distribution of the phasor configurations which will be sampled
in numerical simulations is given by

$\displaystyle\mathcal{P}[\bm{a}]\propto
e^{-\beta\mathcal{H}[\bm{a}]}\delta\left(\epsilon
N-\sum_{k=1}^{N}|a_{k}|^{2}\right),$ (3.5)

where $\beta$ is the inverse of the spectral power of noise $T$. We notice
that, by rescaling the variables as $\tilde{a}_{k}=a_{k}/\sqrt{\epsilon}$, the
new variables are constrained on a fixed hypersphere at the cost of
introducing the effective inverse temperature $\beta_{\rm
photonic}=\beta\epsilon^{2}=\mathcal{P}^{2}$, which corresponds to the
photonic temperature introduced in Eq. (2.7). In these rescaled variables the
probability distribution of configurations reads as

$\displaystyle\mathcal{P}[\tilde{\bm{a}}]\propto
e^{-\mathcal{P}^{2}\mathcal{H}[\tilde{\bm{a}}]}\delta\left(N-\sum_{k=1}^{N}|\tilde{a}_{k}|^{2}\right),$
(3.6)

making explicit the role of the pumping rate (i.e. a parameter accounting for
both noise and external pumping) as the true control parameter of the system.
Now, $\mathcal{P}$ can be tuned in numerical simulations either by varying the
effective temperature $T=\beta^{-1}$ and fixing the optical power $\epsilon$,
or by working at fixed temperature and varying the value of $\epsilon$.
Simulations are performed at $\epsilon=1$ varying the temperature $T$ in order
to have a clear correspondence with the literature on glassy systems, but
results are often described in terms of pumping rate $\mathcal{P}$. One simply
needs to remember that, since $\epsilon=1$, the photonic temperature reduced
to the spectral power of noise $T$, and so $\mathcal{P}=1/\sqrt{T}$.

#####  3.1.1 Topological Properties

In this section, we aim to provide some details about the topology of the
interaction network of the ML 4-phasor model. A mode-locked graph can be
defined in full generality as a hypergraph whose hyperedges are selected
according to the FMC (2.36). Equivalently a mode-locked graph can be also
defined on a factor graph, with fixed connectivity of the function nodes and
connectivity of the variable nodes determined by the FMC. The following
treatment is restricted to the case of interest, which is characterized by
comb-like frequencies and 4-body interactions, but it can be extended to more
general situations.

An interesting quantity to compute is the total number of hyperedges that are
left by the FMC in the interaction graph of the ML 4-phasor model, with
respect to the fully-connected case. Let us denote by $N^{(f)}_{4}$ the number
of tetrads in the fully connected graph, which is given by

$\displaystyle N^{(f)}_{4}=\frac{N(N-1)(N-2)(N-3)}{4!}\sim\frac{N^{4}}{4!},$
(3.7)

with $N\gg 1$. We notice that even if the adjacency tensor defined by the FMC
with equispaced frequencies is not completely symmetric under permutations of
the indices, each term entering the Hamiltonian (3.1) has some symmetry. The
condition (3.4) can be satisfied by 24 permutations of the indices, which can
be grouped into 3 independent orderings with 8 equivalent permutations each.

Given a tetrad of indices $\bm{k}=\\{k_{1},k_{2},k_{3},k_{4}\\}$, the 3 non
equivalent orderings in Eq. (3.4) can be chosen to be

  1. 1.

FMC1. $\mathcal{P}_{1}=(k_{1},k_{2},k_{3},k_{4})$ identifying the combination
$k_{1}-k_{2}+k_{3}-k_{4}=0$ and all indices permutations;

  2. 2.

FMC2. $\mathcal{P}_{2}=(k_{1},k_{3},k_{2},k_{4})$ identifying
$k_{1}-k_{3}+k_{2}-k_{4}=0$ and all indices permutations;

  3. 3.

FMC3. $\mathcal{P}_{3}=(k_{2},k_{1},k_{3},k_{4})$ identifying
$k_{2}-k_{1}+k_{3}-k_{4}=0$ and all indices permutations;

All the permutations inside each of the 3 groups correspond to terms inside
the Hamiltonian (3.1) which have the same value of the energy. Consider for
example the first ordering:

$\displaystyle\text{FMC}_{1}:\omega_{k_{1}}+\omega_{k_{3}}=\omega_{k_{2}}+\omega_{k_{4}}\rightarrow
k_{1}+k_{3}=k_{2}+k_{4},$ (3.8)

where we have used Eq. (2.2). Following Ref. [Mar+18], we consider uniformly
distributed indices, i.e. $P(k)=1/N$ for $k\in[1,N]$. In this case, the
probability for the sum of two indices $k_{ij}^{+}=k_{i}+k_{j}$ to take the
value $k^{+}\in[2,2N]$ can be easily determined:

$\displaystyle P_{+}(k^{+})=\left\\{\begin{array}[]{rcl}\begin{aligned}
&\frac{k^{+}-1}{N^{2}}&\mbox{if}&~{}~{}~{}k^{+}\in[2,N+1]\\\
&\frac{2N-(k^{+}-1)}{N^{2}}&\mbox{if}&~{}~{}~{}k^{+}\in[N+2,2N].\end{aligned}\end{array}\right.$
(3.10)

We can now evaluate the probability that the quadruplet $\bm{k}$ satisfies the
condition $\mbox{FMC}_{1}$ as the probability that the left and right hand
side of Eq. (3.8) take the same value $k^{+}$:

$\displaystyle P(\mbox{FMC}_{1})$
$\displaystyle=\sum_{k^{+}=2}^{2N}P_{+}(k^{+})^{2}$
$\displaystyle=\frac{1+2N^{2}}{3N^{3}}\sim\frac{2}{3N},$

where the last relation holds in the large-$N$ limit. The same occurs for the
other independent orderings, say $\mbox{FMC}_{2,3}$. This leads to the removal
of the factor $1/3$ from the probability that the FMC is satisfied by any
ordering. Eventually, the number of couplings in the interaction network of
the ML 4-phasor model is

$\displaystyle
N_{4}^{*}=\frac{2}{N}\left[1+\mathcal{O}\left(\frac{1}{N}\right)\right]N_{4}^{(f)}.$
(3.11)

Hence, the FMC tends to cut $\mathcal{O}(N)$ couplings from the fully-
connected graph, reducing the total number of couplings in the system to
$\mathcal{O}(N^{3})$. This prevision has been verified numerically with great
accuracy. It is worth stressing that, though diluted, the graph is still
dense.

The FMC condition is a deterministic selection rule, which induces non-trivial
correlations among the interacting modes. To gain an insight into the kind of
correlations, we consider an analogy with random networks, following Ref.
[AIL15]. A way to build random but correlated networks is to introduce a
distance among the nodes: in the case of a random graph, e.g. Erdös-Rényi
graph, the distance can be chosen as the absolute value of the difference of
the node indices $d_{k_{1}k_{2}}=|k_{1}-k_{2}|$. Then, one can select bonds
accordingly to a probability that depends on that distance. In the case of the
ML 4-phasor model, which has factor nodes of connectivity 4, one needs a
4-indexed metrics in order to select the interacting quadruplets. This metrics
can be taken as $d_{\bm{k}}=|k_{1}-k_{2}+k_{3}-k_{4}|$: therefore, the FMC
with equispaced frequencies is equivalent to including only quadruplets
presenting the minimum value, $d_{\bm{k}}=0$. In this way, the mode
frequencies are not degrees of freedom, but coordinate driving correlations
playing the role of a distance on a graph. It should be stressed again that in
the present case the mode coupling is deterministic and not random: no
probability is associated with the distance. As a result, modes that are in
the center of the spectrum are preferred for combinatorial reasons. Indeed,
the central modes have a higher probability of having close frequencies in the
sense of the distance $d_{\bm{k}}$. This is the reason for the narrowing of
the intensity spectrum observed through the lasing transition, see e.g. Fig.
1.1, when the external pumping is increased (or equivalently the temperature
is reduced).

#####  3.1.2 Generation of a Mode-Locked Graph

Let us here describe in detail how the FMC is implemented in our code in order
to build the interaction network of the ML 4-phasor model.

First, a virtual complete graph with
$N_{4}^{(f)}=\binom{N}{4}\sim\mathcal{O}(N^{4})$ interactions is generated
with ordered quadruplets of indices $k_{1}<k_{2}<k_{3}<k_{4}$. Then the FMC is
applied to the complete graph. Notice that for each ordered quadruplet of
indices, the FMC can be satisfied only by the permutation class
$\mathcal{P}_{3}(k_{1},k_{2},k_{3},k_{4})$. Each time a quadruplet of indices
matches the previous condition, the corresponding interaction is added to the
real graph and a random value extracted from the Gaussian distribution Eq.
(3.2) is assigned to it. This procedure is repeated by randomly picking a
quadruplet from the complete graph until a preassigned number
$N_{4}\sim\mathcal{O}(N^{3})$ of interactions for the ML graph is reached.

In order to be able to perform a neat finite-size scaling analysis, the number
$N_{4}$ is chosen to be the largest power of $2$ below the total number of
couplings satisfying the FMC, which is given by the quantity $N_{4}^{*}$
computed in Eq. (3.11). In practice, the number $N_{4}$ is chosen first and
then the corresponding size $N$ to be simulated is selected in order to
minimize the difference $\Delta N=N_{4}-N_{4}^{*}$.

We notice that this way of building the mode-locked interaction network
introduces an artificial source of disorder in the model, besides the original
one. Each one of the $N_{\text{s}}$ simulated disordered samples is
characterized by a realization of the couplings $\\{J_{\bm{k}}\\}$ that
differs from the others both in the quadruplet network and in the numerical
values. However, the fluctuations of the observables with respect to the
randomness of the quadruplet topology turn out to be much smaller than the
fluctuations due to the numerical values of the couplings, already for the
smallest simulated sizes. In fact, as it has been observed also for the
ordered mode-locked graph [Ant16], when compared on the log-scale the energy
fluctuations occurring during the equilibrium dynamics (i.e. the ensemble
fluctuations) are at least two orders of magnitude larger than the graph-to-
graph fluctuations, which are therefore negligible for all practical purposes.

####  3.2 Numerical Analysis

This section is devoted to present the details of the Monte Carlo algorithm
implemented for the simulation of the model.

#####  3.2.1 Exchange Monte Carlo Algorithm

The numerical simulations of the ML 4-phasor model have been performed by
means of an Exchange Monte Carlo algorithm parallelized on GPU’s 222Graphic
Processing Units. The code, written in CUDA, has been running on three types
of GPU: Nvidia GTX680 (1536 cores), Nvidia Tesla K20 (2496 cores) and Nvidia
Tesla V100 (5120 cores). in order to sample the equilibrium probability
distribution Eq. (3.5). In this section, we briefly describe the most salient
features of the method, by following Ref. [NB99]; then, we provide a few
details on our implementation for the simulation of the ML 4-phasor model.

The Exchange Monte Carlo method, else called Parallel Tempering (PT), was
introduced by Hukushima and Nemoto in Ref. [HN96] as a variation of simulated
tempering [MP92]. In fact, PT is the most simple and general form of simulated
tempering, which has been proposed as a finite-temperature generalization of
the famous simulated annealing [SV83]. All these algorithms have been
developed in order to cope with complex optimization problems, characterized
by the presence of many minima of the cost function, where usual Monte Carlo
algorithms are not feasible. In particular, PT has been widely used for
equilibrium simulations of finite-dimensional spin glasses, see e.g. Refs.
[Bal+00, HPV08, PM14], which are known to be “hardly-relaxing” systems.
Indeed, for glassy systems standard iterative algorithms tend to get stuck in
small regions of the state space from which they cannot escape, facing the so
called critical slowing down.

One explanation for this phenomenon is that each state in the Markov chain of
a Monte Carlo algorithm is chosen from the previous one and is in some sense
close to it. Therefore, starting from a certain initial configuration, there
are states that can be reached with a small number of moves, while there are
other states, farther in the configuration space, which can only be reached in
a large number of moves. The state space of glassy models contains many stable
and metastable states, as it can be revealed by several techniques
333Metastable states can be revealed by directly studying the model dynamics,
typically in the Martin-Siggia-Rose functional integral formalism (developed
by De Dominicis and Janssen for disorderd models, see e.g. Ref. [De ̵78]), but
also by the study of the complexity or by an analysis based on the Thouless-
Anderson-Palmer (TAP) equations [TAP77]. The prototype case is that of the
spherical p-spin model, for which the dynamical equations have been closed and
solved in Ref. [CHS93, CK93] and the TAP approach has been developed in Refs.
[KPV93, AH95]. Moreover, the structure of the metastable states has been
carefully analyzed by means of the Franz-Parisi potential in Refs. [FP95,
CGP97], while their basins of attraction have been studied in Ref. [BF98].. It
turns out that the metastable states have a relatively low energy with respect
to the states by which they are surrounded. If a simulated system is
initialized in a configuration close to a metastable state (or directly in a
stable state), to escape its basin of attraction, the algorithm must pass
through one of the surrounding states with higher energy, an occurrence that
has an exponentially low probability, since configurations are sampled with
Boltzmann weights, which depend on the energy difference.

It has to be pointed out that, of course, multiple states are also present in
much simpler systems, such as the standard Ising model below the critical
temperature. In this case there are only two low temperature states in which
the Gibbs measure breaks down when the system size is sent to infinity: a
positively magnetized state and a negatively magnetized one. In fact, to pass
from one state to the other, the system configuration has to jump an energy
barrier whose height scales exponentially with the size of the system. This
event has a probability that is exponentially small, similarly to the
corresponding case in spin-glass models. However, in this case the low
temperature states are symmetric under spin reversal transformation and one
gets the same information on the measured properties of the system, no matter
whether the initial configuration is chosen close to a state or to the other.
Conversely, in glassy models the states are usually not related by any
symmetry, so it happens that for different simulations the algorithm gets
stuck in different basins each time depending on the initial condition, giving
completely different answers for the observables. This can be regarded as a
finite-size evidence of ergodicity breaking.

To avoid this situation, PT has been developed based on the idea that system
thermalization may be facilitated by a reversible Markovian dynamics of
configurations among heat baths at close temperatures. In particular,
configurations belonging to copies of the system at higher temperature may
help the copies at lower temperature to jump out of the local minima of the
rugged free energy landscape. While the dynamics is carried out in parallel
for all the heat baths simulated, once after a fixed number of steps an
exchange of configurations between baths at neighboring temperatures is
proposed. We refer to this kind of move as swap, to distinguish it from the
usual Monte Carlo step. A swap is proposed sequentially for all pairs of
neighboring inverse temperatures $\beta_{i}$ and $\beta_{i+1}$, with the
following acceptance probability implementing detailed balance with the
equilibrium Boltzmann distribution for each thermal bath:

$\displaystyle p_{\text{swap}}=\min\ [1\ ,\
e^{(\beta_{i}-\beta_{i+1})(\mathcal{H}[\bm{a}_{i}]-\mathcal{H}[\bm{a}_{i+1}])}].$
(3.12)

Thus, each state in each of the simulations is sampled with exactly its
Boltzmann weight, so that in PT simulations measurements can be performed in
the same way as in a usual Monte Carlo simulation. We will come back on the
measurement process in the following.

It should be noted that a key role is played by the time after which a move is
proposed. This time has to be large enough in order to avoid exchanges among
configurations that are very similar to those just exchanged, and, on the
other hand, not too large, otherwise thermalization will require a very long
time. In one word, one needs the time between two subsequent swaps to be the
smallest possible in order to make the most of a PT algorithm.

#####  3.2.2 The choice of temperatures

The reason why the PT method overcomes energy barriers is strictly related to
the choice of the simulation temperatures. To give a more intuitive
understanding of the situation, we go through the following argument, taken
from Ref. [NB99].

Let us focus of two copies of the system at temperatures $T_{1}<T_{2}$, one
below and the other above the glass transition of the system. Suppose that the
two copies start from configurations which belong to the same energy basin.
Since the high temperature simulation do not show ergodicity breaking, it will
freely explore the phase space on a time scale similar to that of a simulation
of a normal, non-glassy system. On the contrary, the low temperature copy of
the system will remain stuck in the initial energy basin. If one attempts to
swap the states of the two simulations, Eq. (3.12) says that, unless the swap
does not increase the energy of the low temperature simulation by a great
deal, then it is unlikely to be accepted. However, from time to time, it will
happen that the system at $T_{2}$ finds its way into a region of low energy,
that is another basin with respect to the initial one, where the simulation at
$T_{1}$ is stuck. When this happens a swap will quite likely be accepted. Thus
the low temperature copy is transported in one move to another energy basin,
and the high temperature one finds itself back in the basin that it started
in. By repeating the process over a long time, the low temperature simulation
is moved repeatedly to new energy basins. Thus, the PT algorithm effectively
overcomes the problem of barrier crossing, which makes simulation of glassy
systems so hard, and allows us to sample a significant fraction of the state
space, while still sampling with the correct Boltzmann weights for a
temperature below the glass transition [NB99].

In view of this argument, before running the simulations one must have an
approximate knowledge of the critical temperature of the model, in order to
establish properly the temperature interval, i.e. define a
$\beta_{\text{min}}$ and a $\beta_{\max}$ such that the critical temperature
falls inside the interval. Moreover, one has to bear in mind that temperatures
should be close enough, so that the typical configuration domains at nearby
temperatures overlap. If this does not occur, the energy distributions at some
nearby heat-baths might display no sensitive overlap, thus yielding an
extremely low probability of a swap between them. If a critical point is
there, this is likely to occur when one heat-bath is at a temperature above
the critical one and the other one at a temperature below it. If this is the
case there will be a drastic drop in the exchange frequencies (swapping rate)
between these two temperatures, above and below the phase transition, making
the algorithm extremely inefficient.

One essential criterion to decide if the algorithm is working efficiently,
both regarding on how often we propose a swap move and the choice of
temperatures, is to compute the swapping rate between adjacent temperatures,
that is the fraction of accepted swaps. The algorithm works efficiently only
if, for all couples of temperatures is not too small. An optimal value lies
between 0.6 and 0.8 and this interval has been taken as a reference in this
work.

#####  3.2.3 Details of the Algorithm

In what follows, the specific PT algorithm designed for the simulation of the
ML 4-phasor model is described. The first step is to build the mode-locked
graph of interaction is built according to the procedure described before.
Then a PT dynamics is run for copies of the system with the same quenched
disorder at different temperatures. Each of the system copies follows its own
dynamics, except when a swap is accepted among neighboring heat baths. In the
following, we focus on the main features of the algorithm.

###### Local update

The algorithm uses local Metropolis updates for the dynamics of each PT copy
of the system. Then, a configuration update has to be proposed with the
requirement to keep

$\displaystyle\sum_{k=1}^{N}|a_{k}|^{2}=\text{const}.$

In order fulfill the constraint, each update of the configurations is carried
out by choosing at random two variables $a_{i}=A_{i}e^{i\phi_{i}}$ and
$a_{j}=A_{j}e^{i\phi_{j}}$ and extracting three random numbers [Ant16, GAL20]:
$x,y\in[0,2\pi]$ and $z\in[0,\pi/2]$. The first two correspond to the new
(attempted) phases $x=\phi_{i}^{\prime}$ and $y=\phi_{j}^{\prime}$, the third
one mixes the intensities of the modes selected by preserving the spherical
constraint

$\displaystyle A_{j}^{\prime}=\sqrt{A_{i}^{2}+A_{j}^{2}}\cos
z~{}~{}~{}~{}A_{j}^{\prime}=\sqrt{A_{i}^{2}+A_{j}^{2}}\sin z$

Then the attempted update is accepted according to the usual Metropolis
formula, in order to implement detail balance.

###### Parallel computation

It is worth noticing that for the case of the ML 4-phasor model the updates
have to performed sequentially, rather than in parallel. In fact, the parallel
Monte Carlo algorithm of a system of interacting variables needs a sparse
network of interaction, such as nearest neighbors, e.g., to be implemented, in
order to split the system in smaller non-interacting sub-systems, that can be
updated in parallel. This procedure clearly speeds up the computation of the
update. In our system, however, this is not possible due to the density of the
mode-locked interaction network, in which each variable participate in
$O(N^{2})$ interacting quadruplets.

However another kind of code optimization can be implemented for the ML
4-phasor model, by exploiting the computing capability of GPU’s: the
parallelization of the energy computation in the local Metroplis update. In
order to accept or reject the update of two spins $a_{i}$ and $a_{j}$, one has
to compute the energy difference between the attempted configuration and the
current one. This operation has a computational complexity which scales like
the number of quadruplets involved in the computation, i.e.
$N_{4}^{(i,j)}=O(N^{2})$:

$\displaystyle\Delta E=\sum_{k=1}^{N_{4}^{(i,j)}}\Delta E_{k},$

where $\Delta E_{k}$ denotes the energy difference between each quadruplet.
The computation of each $\Delta E_{k}$ is realized in parallel on a distinct
kernel on GPU.

Last, but not least, also the PT dynamics at different temperatures has been
parallelized on GPU’s. The two kinds of parallelization considered together
reduce the execution time of the entire simulation by a factor of 8 [GAL20].

#####  3.2.4 Observables of Interest

Before listing the observables considered, let us clarify the measurement
procedure. In order to properly estimate statistical errors, time correlations
have been taken into account. A correlation time $\tau_{\rm corr}$ can be
identified as the maximum among all the correlation times of each thermal bath
dynamics. Consequently, the observables can be measured every $\tau_{\rm
corr}$ Monte Carlo steps. If $N_{\rm MCS}$ is the total amount of Monte Carlo
steps of the simulation, for each disordered sample the number of thermalized,
uncorrelated configurations is given by

$\displaystyle\mathcal{N}\equiv\frac{N_{\rm MCS}-\tau_{\rm eq}}{\tau_{\rm
corr}},$ (3.13)

where $\tau_{\text{eq}}$ denotes the thermalization time of the replica with
the lowest temperature. We will come back in a while on how the equilibrium
time is defined. Then, for a given observable $O$ function of the
configurations $\bm{a}$, the ensemble average is estimated by the following
time average

$\displaystyle\langle O[\bm{a}]\rangle=\frac{1}{\mathcal{N}}\sum_{t=\tau_{\rm
eq}/\tau_{\rm corr}}^{N_{\rm MCS}/\tau_{\rm corr}}O[\bm{a}_{t}].$ (3.14)

On top of that, the average over the quenched randomness of the couplings is
defined as follows. For each $\\{J_{\bm{k}}\\}$ realization we have a thermal
average $\langle O[\bm{A}]\rangle_{J}$. Averaging over the random samples
yields the least fluctuating finite-$N$ proxy for the average in the
thermodynamic limit:

$\displaystyle\overline{O[\bm{a}]}=\frac{1}{N_{\text{s}}}\sum_{j=1}^{N_{\text{s}}}\langle
O[\bm{a}]\rangle_{J}^{(j)}.$ (3.15)

The observables that will be considered in the present chapter and in the next
one are the following:

  * •

_intensity spectrum_ : normalized to the square root of the temperature, in
order to connect with the physical intensities

$\displaystyle I_{k}=\frac{A_{k}^{2}}{\sqrt{T}},$ (3.16)

  * •

_specific heat_ : measured as the equilibrium energy fluctuations as

$\displaystyle
c_{V_{N}}=\frac{1}{N}\frac{\overline{\langle\mathcal{H}^{2}\rangle-\langle\mathcal{H}\rangle^{2}}}{T^{2}},$
(3.17)

  * •

_Parisi overlap distribution_ , $P(q)$, where for the ML 4-phasor model the
overlap among configurations (see Eq. (2.61)) is given by

$\displaystyle q_{\alpha\beta}$ $\displaystyle=\frac{1}{N}\mbox{Re
}\sum_{k=1}^{N}\overline{a}_{k}^{\alpha}a_{k}^{\beta}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}A_{k}^{\alpha}A_{k}^{\beta}\cos(\phi_{k}^{\alpha}-\phi_{k}^{\beta}),$
(3.18)

  * •

_plaquette overlap distribution_ , $P(\mathcal{Q})$, where the overlap among
the plaquettes of two replicas is defined as [GAL20]

$\displaystyle\mathcal{Q}_{\alpha\beta}=\frac{1}{N_{4}}\sum_{\bm{k}}\mathcal{E}_{\bm{k}}^{\alpha}\mathcal{E}_{\bm{k}}^{\beta}$
(3.19)

with the plaquette overlap given by

$\mathcal{E}_{\bm{k}}^{\alpha}=a_{k_{1}}\overline{a}_{k_{2}}a_{k_{3}}\overline{a}_{k_{4}}+\text{c.c.}$
(3.20)

and $N_{4}$ the number of quadruples satisfying the FMC, as discussed in Sec.
3.1.2.

  * •

_intensity fluctuations overlap distribution_ , $P(\mathcal{C})$, where the
intensity fluctuation overlap (IFO) among two replicas of the system (see Eq.
(2.74)) is given by

$\displaystyle\mathcal{C}_{\alpha\beta}=\frac{1}{N}\sum_{k=1}^{N}\Delta_{k}^{\alpha}\Delta_{k}^{\beta},$
(3.21)

where the intensity fluctuations are defined as

$\displaystyle\Delta_{k}^{\alpha}=\frac{I_{k}^{\alpha}-\langle
I_{k}^{\alpha}\rangle}{2\sqrt{2}\epsilon}.$ (3.22)

###### Thermalization

In order to guarantee that the data used to compute the displayed observables
are taken from correctly equilibrated samples thermalization can be tested in
several ways.

First, one can look at energy relaxation on sequential time windows whose
length is each time twice the length of the previous one. For each simulated
heat bath dynamics, a minimal requirement is that the time average of the
energy $\langle\mathcal{H}\rangle$ takes the same value at least on the last
and second but last windows. A similar test is performed on the specific heat,
by computing energy fluctuations over the last and the second-last
“logarithmic” window and checking that the values obtained for each
temperature match inside their statistical errors.

However, one can not rely solely on the trend of the energy and of its
fluctuations over time to assess thermalization. In fact, as mentioned before,
due to the presence of energy barriers that scale exponentially with the
system size, one may mistake a single local minimum for a good equilibrium
state. Therefore, a further and stronger requirement for thermalization is the
symmetry of the Parisi overlap distribution $P_{J}(q)$ for each disordered
sample, which can be tested by checking that its skweness is approximately
zero in the low temperature phase.

Once dynamical thermalization to equilibrium has been assessed and a
thermalization time $\tau_{\rm eq}$ identified, the time average Eq. (3.14)
coincides with the canonical ensemble average. The number of Monte Carlo steps
necessary to reach thermalization for each simulated size are reported in
Tables 4.1 and 4.2, for the simulations performed to obtain the results of
Ref. [Nie+22].

####  3.3 Evidence of a Random First Order Transition

In this section we present the results of numerical simulations, where
evidence of a glass transition in the ML 4-phasor model has been first
reported [GAL20]. Here, we are just interested in discussing the physical
picture drawn from these simulations: more technical details on measurement
and data analysis will be provided in the next chapter, where the results of
new simulations are presented.

As already mentioned in the previous chapter, a Random First-Order Transition
(RFOT), the paradigm of the glass transition [KW87, LW07, LN08], is a mixed-
order phase transition, characterized by the divergence of the thermodynamic
susceptibilities and, at the same time, by the discontinuity of the order
parameter at the static transition point. The former is the signature of a
critical phenomenon, which is determined by a continuous second-order phase
transition; the latter is, instead, a feature which is typical of first-order
transitions, where the new dominant thermodynamic state is already present
before the transition, differently from the continuous case, where it arises
at the transition.

The observables which help to investigate the presence of a RFOT444Remember
that it is the static transition to be relevant for experiments on random
lasers, and not the dynamic one, as it usually is in structural glasses. This
is why our attention is devoted to the simulation of the equilibrium
properties of the model at the static glass transition. in the ML 4-phasor
model are the specific heat (3.17) and the overlap distribution function (• ‣
3.2.4). A singularity in the specific heat puts in evidence the second-order
nature of the transition, whereas a jump in the overlap probability
distribution $P(q)$ is a signature of its first-order nature. In models with
continuous variables, the $P(q)$ is expected to be a distribution with a
single peak in $q=0$ in the high temperature phase and to develop side peaks,
as well, in the low temperature glassy phase. At finite $N$, of course, exact
Dirac delta peaks in the $P(q)$ appear as smoothed functions of $q$, due to
finite-size effects. One does not have to confuse this finite-size behavior of
the $P(q)$ in a RFOT, with the behavior of the $P(q)$ in the spin-glass
transition of the SK model [SK75], where the overlap distribution is expected
to take a non-trivial shape, different from a bimodal one, even in the
$N\rightarrow\infty$ limit.

Figure 3.1: (a) Specific heat
$c_{V_{N}}=\overline{\langle\mathcal{H}^{2}\rangle-\langle\mathcal{H}\rangle^{2}}/(NT^{2})$
as a function of $T$; Different curves represent different system size.
(Inset) Specific heat as a function of $\tau N^{3/2}$ where
$\tau=(T/T_{c}(N)-1)$: curve collapse in the scaling region. The four sizes
are $N=32,48,64,96$. (b) Configuration overlap probability distribution for
$N=64$ at temperatures $T/T_{c}=1.71,1.15,0.97,0.64$. Reprinted from [GAL20].

In Fig. 3.1, we display the behaviour of the specific heat and of the
configuration overlap distribution function. The specific heat diverges as the
size increases: in the inset panel data are collapsed in the critical region
with an exponent $3/2$, which turns out not to be compatible with a mean-field
theory of second-order phase transitions. The reason why this comes about will
be clarified in the next chapter: there, we will see how this unexpected
exponent $3/2$ turns out to be a preasymptotic finite size effect.

The configuration overlap distribution function in Fig. 3.1 turns out to be
Gaussian in the low-$\mathcal{P}$ phase. Then, for
$\mathcal{P}>\mathcal{P}_{c}$, the distribution shows a clear deviation from
Gaussianity, but only “shoulders” are displayed at the simulated sizes, rather
than proper side peaks.

These results can be compared with those obtained through numerical
simulations of a 4-phasor model with random dilution of the same order of that
induced by the FMC, see Ref. [GAL20]. The comparison reveals two important
differences: first, the scaling of the critical region in the case of random
dilution yields an exponent of $1/2$ which perfectly matches the expectations
of standard $\phi^{4}$ mean-field theory; secondly, in the case of random
dilution the overlap distribution function exhibits clear secondary peaks at a
finite distance from the origin in the high-$\mathcal{P}$ phase, signaling a
glassy RSB phase. It is evident, then, that the finite-size effects are
stronger in the mode-locked model than in the randomly diluted one. This is
also quite intuitive: at the small simulated sizes the correlations induced by
a deterministic selection rule are not negligible. In fact, the diluted mode-
locked graph of interactions, though comprised by an extensive number of
couplings less than the fully-connected graph, is still dense, suggesting
compatibility with mean-field theory. However, up to the precision of the
study reported here, the question whether the ML 4-phasor model is a mean-
field theory or not remains open and needs a more refined analysis to be
answered.

Figure 3.2: (a) Plaquette overlap distribution $P(\mathcal{Q})$ for $N=64$,
computed over $2^{14}$ quadruplets and plotted for four different values of
the pumping rates: $(\mathcal{P}_{c}/\mathcal{P})^{2}=1.71,1.15,0.97,0.64$
(Inset) Multimodality parameter $b$ measured for $P(\mathcal{Q})$ as a
function of $T=\mathcal{P}^{-2}$: values above the threshold $b^{*}=5/9$ (full
black line) indicate a bimodal distribution. (b) Intensity fluctuation overlap
(IFO) probability distribution $P(\mathcal{C})$ for system size $N=64$,
plotted for the same values of the pumping rate as panel (a). (Inset) Binder
parameter $\mathcal{B}$ measured for $P(\mathcal{C})$ as a function of $T$ ;
the behavior is typical of first-order transitions, with the transition at the
minimum of $\mathcal{B}$. Reprinted from [GAL20].

Regarding the first-order nature of the transition a stronger indication comes
from the study of the plaquette overlap distribution (3.19). At variance with
the configurational overlap, which is computed over $N$ variables, the
plaquette overlap is computed over $\mathcal{O}(N^{3})$ quadruplets, hence it
is less plagued by finite-size effects. In Fig. 3.2 we display both the
plaquette overlap and the IFO (3.21) probability distribution functions. In
the low-$\mathcal{P}$ phase the plaquette overlap has a very peaked
distribution in $\mathcal{Q}\simeq 0$, while for
$\mathcal{P}\approx\mathcal{P}_{c}$ there is clear evidence of a secondary
peak at $\mathcal{Q}>0$. A similar behavior is shown by the IFO distribution
function: in this case, although at first sight there is no clear evidence of
secondary peaks at high pumping rates, for the same reason why this happens in
the $P(q)$, non-Gaussian tails appear in the vicinity of the transition.

To deepen the analysis one can study the moments of these distributions. In
particular, it is useful to consider the third and fourth moments, which are
related to the symmetry of the distribution and to the swelling of its tails.
By denoting with $q$ a general overlap (be it the Parisi, the plaquette or the
intensity fluctuation overlap) the skewness and the kurtosis of its
distribution function are given by

$\gamma=\frac{\overline{\Delta q^{3}}}{(\overline{\Delta
q^{2}})^{3/2}}~{}~{}~{}~{}~{}~{}~{}~{}k=\frac{\overline{\Delta
q^{4}}}{(\overline{\Delta q^{2}})^{2}}.$ (3.23)

In the high temperature phase, where the order parameter is zero, the
distribution at finite size will be a Gaussian centered in zero, with $k=3$.
At the phase transition and in the low temperature phase the distribution will
be different from a Gaussian. A very useful quantity to measure the deviation
from Gaussianity is the Binder parameter [Bin81, BL84], which is defined as

$\mathcal{B}=\frac{1}{2}\left(3-k\right).$ (3.24)

If we are dealing with a first-order phase transition, then the Binder
parameter displays a nonmonotonic reversed bell behavior with a maximal
deviation from Gaussianity in the coexistence region, where the distribution
is bimodal. This is precisely what is observed for the IFO distribution
function and has been reported in the inset of panel (b) in Fig. 3.2, where
$\mathcal{B}$ is plotted as a function of the effective temperature
$T=\mathcal{P}^{-2}$.

As a term of comparison, we refer to Ref. [PRS01], where the Binder parameter
is computed in numerical simulations of both the SK model and the Ising
$p$-spin model (with $p=3$). In the former case, characterized by a continuous
transition to a FRSB phase, the Binder parameter is zero at high temperature
and then increases monotonically. Moreover, the Binder parameter for different
sizes exhibits a crossing of the curves for different sizes, from which the
critical temperature of the model can be estimated (cf. also Ref. [BY88] for
the case of the Edwards-Anderson model in two or three dimensions). On the
other hand, in the case of the Ising $p$-spin model, characterized by a
discontinuous transition to a 1RSB phase (regarding the distribution of the
order parameter) the Binder parameter exhibits the reverse bell behavior that
has been here reported in the case of the IFO probability distribution
function.

For the study of the plaquette overlap distribution, one can not use the
Binder parameter as a good indicator because, although clearly bimodal at the
transition, the distribution is not Gaussian far from the transition. However,
in order to study the bimodal nature of the distribution at the transition,
one can introduce another parameter, the so-called _bimodality_ parameter,
which is defined as

$b=\frac{\gamma^{2}+1}{k+\frac{3(n-1)^{2}}{(n-2)(n-3)}}$ (3.25)

where $n$ is the number of data composing the histogram of the probability
distribution. In the inset of panel (b) in Fig. 3.2 the behavior of $b$ as a
function of temperature is reported: the region where the parameter $b$
signals a bimodal distribution of the overlap is precisely the interval of
pumping rates around $\mathcal{P}_{c}$.

###  Chapter 4 Universality Class of the Transition

In this chapter, the study of the static glass transition in the ML 4-phasor
model is improved with respect to the previous analysis (Chap. 3). The main
goal is to determine the universality class of the model, through a refined
finite-size scaling analysis of the transition. In particular, we aim to
assess whether the unexpected scaling exponent of the critical region found in
Ref. [GAL20] is a genuine non-mean-field feature of the mode-locked diluted
model or if it is only a preasymptotic effect, due to the very strong finite-
size effects, which derive from the difficulty to simulate large enough sizes
in dense disordered models.

First, a preliminary argument for the scaling of the critical region in mean-
field theory is presented in a very general way, which does not require any
specific knowledge about the glass transition, besides its second-order nature
concerning the divergence of the susceptibilities at the critical temperature.
The outcome of the analysis is that the exponent for the scaling of the
critical region in a generic mean-field model must fall in a compact interval
and its specific value is related to the order of the leading nonlinear term
of interaction. The result obtained is coherent with the scaling exponent of
the Random Energy Model (REM), the simplest model exhibiting a glass
transition [Der80, GM84].

In the second part of the chapter, results obtained from new numerical
simulations are presented [Nie+22]. Besides performing simulations of the
original ML 4-phasor model increasing the number of simulated sizes,
temperatures and samples, a slightly different version of the model is
introduced, in order to reduce the specific kind of finite-size effects
induced by the FMC. In particular, the bond filtering action of the FMC is
considered on a frequency space with periodic boundary conditions, leading to
simulations of the bulk spectrum of the original model. The study of the
critical region in both cases leads to a reduction of the value of the
exponent for the scaling of the critical region, which turns out to fall into
the interval defined by the mean-field argument, inside the experimental
uncertainty. Finally, the study of the glass transition is completed by
analyzing data in order to obtain the overlap distribution functions
introduced in the previous chapter.

####  4.1 A Mean-Field Argument for the Scaling of the Critical Region

Second-order phase transitions are characterized by scale invariance at
criticality: in finite-dimensional systems (e.g. lattices in $d$ dimensions),
fluctuations extend over regions of all possible sizes and a characteristic
scale of length no longer exists. As a consequence, the behavior of physical
quantities near the critical point is described with respect to some control
parameter, e.g. temperature, by power laws, whose exponents are universal, in
the sense that they depend only on very general properties of the system (such
as the dimensionality of space $d$, the dimensionality of the order parameter
and the symmetries of the Hamiltonian), but not on the details of the
microscopic interactions [Zin02]. Out of the mean-field approximation, the
critical exponents can be computed with renormalization group techniques, such
as the $\epsilon$-expansion [Par88], or can be extrapolated from the study of
the system at finite size $L$ (i.e. the linear size of a lattice), through the
finite-size scaling analysis [FB72].

However, here we are interested in the case of infinite-dimensional models
defined on graphs not embeddable into any space with finite dimensions $d$,
for which mean-field theory is exact. Even in this case, the scaling regime
amplitude has a dependence on the system size $N$ (i.e. the number of nodes in
the graph of interactions), governed by an exponent, which we will denote as
$\nu_{\rm eff}$, through the following relation

$\displaystyle|T-T_{c}|\sim 1/N^{1/\nu_{\rm eff}},$ (4.1)

where $T$ is an effective control parameter for the transition - from now on
we will refere to it as temperature - and $T_{c}$ denotes the critical point.
We will make sense of $\nu_{\rm eff}$ in terms of standard critical exponents
in the following. The prototype of a mean-field theory of continuous phase
transitions is represented by the Landau effective potential111Let us stress
that $\phi$ is a global quantity, not a local magnetization field $\phi(x)$ as
in the Landau-Ginzburg $\lambda\phi^{4}$ theory. The potential $V(\phi)$ in
Eq. (4.2) is the result of the Landau approximation of the $\lambda\phi^{4}$
field theory, which consist in taking $\phi(x)=\phi$ for all points in space,
hence neglecting the Laplacian term: this is nothing but a mean-field
approximation.

$\displaystyle V(\phi)=\frac{\tau}{2}\phi^{2}+\frac{g}{4!}\phi^{4},$ (4.2)

where $\phi$ is the global order parameter of the transition, $\tau=T/T_{c}-1$
is the reduced temperature and $g$ is the coupling constant. The probability
distribution of the order parameter $\phi$ is

$\displaystyle p(\phi)=\frac{e^{-NV(\phi)}}{Z},$ (4.3)

where the partition function $Z$ is given by

$\displaystyle Z=\int\text{d}\phi~{}e^{-NV(\phi)}.$ (4.4)

In a second-order transition, the relevant quantities that have to be
considered are the fluctuations of the order parameter around the minimum of
the effective potential

$\displaystyle\delta\phi^{2}=\langle\phi^{2}\rangle-\langle\phi\rangle^{2},$
(4.5)

where the brackets denote the average with respect to the distribution
$p(\phi)$. The critical behaviour of the susceptibilities, which include the
specific heat, is related to the fluctuations of the field near the critical
point. In order to estimate the value of the exponent $\nu_{\rm eff}$ we aim
to match the fluctuations above and below $T_{c}$.

When $\tau\gtrsim 0$ the effective potential is well approximated by
$V(\phi)\approx\frac{1}{2}\tau\phi^{2}$ for values of the order parameter
close enough to the minimum $\phi^{*}=0$. The probability distribution of
$\phi$ is a zero-mean Gaussian distribution and its normalization is given by

$\displaystyle
Z\approx\int\text{d}\phi~{}e^{-\frac{N\tau}{2}\phi^{2}}\sim\frac{1}{\sqrt{N\tau}}.$
(4.6)

In this regime, the fluctuations of the order parameter centred around the
minimum $\phi^{*}=0$ are given by the variance of the Gaussian distribution:

$\delta\phi^{2}_{{}_{T>T_{c}}}\sim 1/N\tau.$ (4.7)

We notice that the previous equation is the usual scaling law for the linear
susceptibility above the critical point in the standard mean-field $\phi^{4}$
theory:

$\chi\sim\frac{1}{N}\tau^{-\gamma},$ (4.8)

with $\gamma=1$, see, e. g. Ref. [Car96].

On the other hand, when $\tau\lesssim 0$ the quartic term of the potential
becomes relevant and can not be neglected. In this regime, the fluctuations of
the order parameter are centered around one of the two symmetric minima of the
potential, namely $\phi_{\pm}=\pm\phi^{*}$, depending on the initial
conditions. Since we are interested in matching the fluctuations above and
below the critical temperature, we assume the temperature to be sufficiently
close to $T_{c}$ in order for the amplitude of the fluctuations to be of the
order of the distance from the origin:

$\delta\phi^{2}_{{}_{T<T_{c}}}\sim(\phi^{*})^{2}.$ (4.9)

Clearly, this is no longer valid well below $T_{c}$, where the curvature of
the potential has to be considered. The minima $\phi_{\pm}$ can be easily
determined according to the saddle-point approximation of the partition
function

$\displaystyle Z=\int\text{d}\phi~{}e^{-NV(\phi)}\approx e^{-NV(\phi^{*})},$
(4.10)

where $\phi^{*}$ such that
$\frac{\text{d}V(\phi)}{\text{d}\phi}\Big{|}_{\phi^{*}}=0$ is

$\phi^{*}=\sqrt{-\frac{6\tau}{g}}.$ (4.11)

Hence, $\delta\phi^{2}_{{}_{T<T_{c}}}\sim-\tau/g$. Therefore, we have an
estimate of the fluctuations on the two sides of the critical point,
respectively

$\displaystyle\delta\phi^{2}_{{}_{T>T_{c}}}$
$\displaystyle\sim\frac{1}{N\tau}$ (4.12)
$\displaystyle\delta\phi^{2}_{{}_{T<T_{c}}}$ $\displaystyle\sim-\tau/g.$
(4.13)

By matching the previous expressions we find

$\displaystyle\delta\phi^{2}_{{}_{T>T_{c}}}\sim\delta\phi^{2}_{{}_{T<T_{c}}}~{}~{}~{}\Longrightarrow~{}~{}~{}|\tau|\sim\frac{1}{N^{1/2}},$
(4.14)

yielding the estimate $\nu_{\rm eff}=2$.

Notice, that Eq. (4.11) corresponds to the usual scaling law for the order
parameter below $T_{c}$ [Car96]

$\phi\sim-\tau^{\beta}$ (4.15)

with exponent $\beta=1/2$ for the mean-field $\phi^{4}$ theory. In terms of
the exponents $\beta$ and $\gamma$ defined by the usual scaling laws, Eq.
(4.14) can be written as

$\tau^{2\beta}\sim\frac{1}{N}\tau^{-\gamma}~{}~{}~{}\Longrightarrow~{}~{}~{}|\tau|\sim\frac{1}{N^{\frac{1}{2\beta+\gamma}}}.$

This allow us to identify the scaling exponent $\nu_{\rm eff}$ for the
critical region of infinite-dimensional models with the expression

$\nu_{\rm eff}=2\beta+\gamma.$ (4.16)

Therefore, $\nu_{\rm eff}$ can be taken just as a shorthand symbol for the
expression $2\beta+\gamma$.

This argument can be straightforwardly extended to more general mean-field
potentials, in order to obtain a range of values for the exponent $\nu_{\rm
eff}$. Let us consider the family of potentials

$\displaystyle V(\phi)=\frac{1}{2}\tau\phi^{2}+\frac{g}{(2n)!}\phi^{2n},$
(4.17)

where the choice of an even lowest order non-linearity is still compatible
with the phenomenology of a second-order phase transition. The fluctuations of
the order parameter above the critical temperature are the same as in the case
$n=2$, which means that the exponent $\gamma=1$ for every $n$. Below the
critical temperature, by using the saddle-point method we find

$\displaystyle\phi^{*}=\left[-\frac{(2n-1)!~{}\tau}{g}\right]^{\frac{1}{2(n-1)}}.$
(4.18)

Incidentally, this means that the general expression of the exponent $\beta$
for the scaling of the order parameter in a $\phi^{2n}$ mean-field theory is
$\beta=1/2(n-1)$. Therefore, the amplitude of the fluctuations above and below
$T_{c}$ are given by

$\displaystyle\delta\phi^{2}_{{}_{T>T_{c}}}$
$\displaystyle\sim\frac{1}{N\tau}$ (4.19)
$\displaystyle\delta\phi^{2}_{{}_{T<T_{c}}}$
$\displaystyle\sim(-\tau/g)^{\frac{1}{(n-1)}},$ (4.20)

and their matching leads to

$\displaystyle\delta\phi^{2}_{{}_{T>T_{c}}}\sim\delta\phi^{2}_{{}_{T<T_{c}}}~{}~{}~{}\Longrightarrow~{}~{}~{}|\tau|\sim\frac{1}{N^{\frac{n-1}{n}}}.$
(4.21)

The range of values that $\nu_{\rm eff}=2\beta+\gamma$ can take in a mean-
field model can be found by taking $n=2$ and $n\rightarrow\infty$ in the
previous expression. Thus, we have

$\displaystyle 1<\nu_{\rm eff}\leq
2~{}~{}~{}\Longleftrightarrow~{}~{}~{}\frac{1}{2}\leq\frac{1}{\nu_{\rm
eff}}<1.$ (4.22)

With respect to this argument, it is clear that the value
$1/\nu_{\text{eff}}=3/2$ found in Ref. [GAL20], is out of the interval found
for mean-field values of the exponent $\nu_{\text{eff}}$.

It is worth noting that the validity of this argument is restricted to the
large-$N$ limit, where the saddle-point approximation holds. This is something
that we have to keep in mind, when comparing results of numerical simulations
at finite $N$ with the estimate on the scaling exponents obtained from the
previous argument.

#####  4.1.1 The Random-Energy Model

In order to compare the previous argument with a well-known model we have
performed a numerical analysis of the Random-Energy Model (REM), which is the
reference mean-field model for disordered systems with non-linear
interactions. Let us briefly review the main features of the model, following
Ref. [MM09], before presenting our results.

The REM is the simplest statistical mechanics model of a disordered system
exhibiting a phase transition [GM84]. It was originally introduced by Derrida
in Refs. [Der80, Der81] as an exactly solvable model arising from the limit
$p\rightarrow\infty$ of the p-spin model. As a result of this procedure, the
model does not take into account any specific interaction among the variables:
the energy is a random process rather than a deterministic function. For this
reason, the REM has the remarkable advantage that results obtained through the
replica method can be compared with formal mathematical approaches, see e.g.
Refs. [GD89, Bov06, OK09, Gue13, Dot11].

Given $M=2^{N}$ energy levels (possibly corresponding to the configurations of
$N$ Ising spin variables) the corresponding energies
$\\{E_{\nu}\\}_{\nu\in\\{1,\dots,M\\}}$ are taken as independent random
Gaussian variables222A generalized version of the model, where correlations
among the energy levels are introduced in a hierarchical way, has been
formulated by Derrida in Ref. [Der85] and exactly solved by Derrida and
Gardner in Ref. [DG86]. Here, we are just interested in the original version
of the model, since the generalized version is in the same universality class.
extracted from the distribution function

$\displaystyle p(E)=\frac{1}{\sqrt{\pi
NJ^{2}}}\exp\left(-\frac{E^{2}}{NJ^{2}}\right),$ (4.23)

where the scaling of the variance with $N$ ensures the extensivity of the
thermodynamic potentials and $J$ is a parameter. An instance of the quenched
disorder corresponds to an extraction of the $M$ energy levels. A Boltzmann
weight $p_{\nu}=\exp(-\beta E_{\nu})/\mathcal{Z}$ is then assigned to each
configuration, where $\mathcal{Z}$ is the partition function of the model

$\displaystyle\mathcal{Z}=\sum_{\nu=1}^{M}e^{-\beta E_{\nu}}.$ (4.24)

We have developed a simple enumeration algorithm to study the REM, which works
as follows: for each disorder sample of a given system size $N$ the energy
levels $\\{E_{\nu}\\}$ are generated, by independently extracting a set of
$2^{N}$ random numbers from the Gaussian distribution Eq. (4.23) with $J=1$. A
set of equispaced temperatures $T$ is generated in the interval
$[T_{\text{min}},T_{\text{max}}]$. The internal energy of the model is
computed as a function of temperature by evaluating the following thermal
average

$\displaystyle\langle E\rangle=\frac{\sum_{\nu}E_{\nu}~{}e^{-\beta
E_{\nu}}}{\sum_{\nu}e^{-\beta E_{\nu}}},$ (4.25)

for each of the $\beta=1/T$ values extracted before. The specific heat is
computed from the fluctuations of the internal energy as

$\displaystyle c_{V_{N,i}}(T)=\frac{1}{N}\frac{\langle E^{2}\rangle-\langle
E\rangle^{2}}{T^{2}},$ (4.26)

where the index $i$ accounts for the sample. The procedure is repeated for
several independent extractions of the random energies $\\{E_{\nu}\\}$.
Eventually, when a sufficiently large number of samples $N_{\text{sam}}$ is
collected for a certain simulated size, the disorder average of the specific
heat is computed by averaging over the samples:

$\displaystyle
c_{V_{N}}(T)=\frac{1}{N_{\text{sam}}}\sum_{i=1}^{N_{\text{sam}}}c_{V_{N,i}}(T).$
(4.27)

The number of samples $N_{\text{sam}}$ is chosen in such a way that the
estimated value of the specific heat remains stable upon fluctuations of the
samples included in the average.

#####  4.1.2 Finite-Size Scaling Analysis

Let us consider momentarily a model defined on a $d$-dimensional lattice of
linear size $L$.

We consider an observable $Y_{L}$ that in the thermodynamic limit scales like

$Y_{\infty}(T)\approx At^{-\psi}$

near the critical point, where $t$ is the modulus of the reduced temperature,
$A$ is some constant and $\psi$ some critical exponent. The fundamental
assumption of the Finite-Size Scaling (FSS) Ansatz [FB72, Car88] is that the
behaviour of $Y_{L}$ near the critical temperature is controlled by the ratio
$\xi_{\infty}/L$. The parameter $\xi_{\infty}$ denotes the correlation length
of the infinite-size system that scales as

$\displaystyle\xi_{\infty}(T)\approx\xi_{0}t^{-\nu},$ (4.28)

where $\xi_{0}$ is a constant. The previous equation is the definition of the
critical exponent $\nu$ governing the scaling of the correlation length. The
scaling hypothesis for $Y_{L}$ can be then written as

$\displaystyle Y_{L}(T)=L^{\omega}f_{Y}(\xi_{\infty}/L),$ (4.29)

where $\omega$ is the critical exponent for the scaling of the peak of the
observable and $f_{Y}$ is a dimensionless function that depends on the
observable $Y$. The function $f_{Y}$ is such that in the limit
$L\rightarrow\infty$ one recovers the scaling law $Y_{\infty}(T)\approx
At^{-\psi}$, therefore $\omega=\psi/\nu$ [FB72]. Moreover, by exploiting Eq.
(4.28) the scaling relation (4.29) can be rewritten as

$\displaystyle Y_{L}(T)=L^{\psi/\nu}\hat{f}_{Y}(L^{1/\nu}t_{L})$ (4.30)

where $t_{L}=|T/T_{c}(L)-1|$, being $T_{c}(L)$ the finite-size critical
temperature, and $\hat{f}_{Y}$ is another scaling function that depend on $Y$.

It is important to notice that, since the lattice size $L$ is linked to the
number of degrees of freedom $N$ by $L=N^{1/d}$, the scaling relation for the
observable $Y$ can be written as

$\displaystyle Y_{N}(T)=N^{\psi/\nu d}\hat{f}_{Y}(N^{1/\nu d}t_{N})$ (4.31)

i.e. in terms of $N$, rather than $L$. By considering the well-known
hyperscaling relation $\nu d=2\beta+\gamma$ [Car96], we find

$\displaystyle
Y_{N}(T)=N^{\psi/(2\beta+\gamma)}\hat{f}_{Y}(N^{1/(2\beta+\gamma)}t_{N}).$
(4.32)

This is the only possible scaling relation that one can use when studying an
infinite-dimensional model, where there is no lattice size $L$, but the only
finite-size parameter is $N$.

In the case of the specific heat, the previous finite-size scaling law takes
the following form (see e.g. [ID89])

$\displaystyle c_{V_{N}}(T)=N^{\alpha/\nu_{\rm
eff}}\hat{f}_{C_{V_{N}}}(N^{1/\nu_{\rm eff}}t_{N}),$ (4.33)

where $\alpha$ denotes the critical exponent of the specific heat peak
divergence, and we have used fact that $\nu_{\rm eff}=2\beta+\gamma$. Since
the dimensionless function $\hat{f}$ is scaling invariant, if one uses the
correct values of the exponents $\alpha$ and $\nu_{\rm eff}$, the curves
$c_{V_{N}}(T)/N^{\alpha/\nu_{\rm eff}}$ for different values of $N$ should
collapse on the same curve.

Figure 4.1: Specific heat of the Random Energy Model. Different colors
represent different simulated sizes. The simulated sizes are $N=16,20,24,28$.
(Inset) Specific heat scaled by $N^{\alpha/\nu_{\rm eff}}$ as a function of
$\tau N^{1/\nu_{\rm eff}}$. The data have been collapsed with exponents
$\alpha=0.52\pm 0.07$ and $\nu_{\rm eff}=1.94\pm 0.22$.

In Fig. 4.1 we show the finite-size behaviour of the specific heat around the
critical temperature. In the main panel, the specific heat is plotted as a
function of temperature, for different system sizes. The simulated sizes are
$N=16,20,24,28$, and for each size $N_{\text{sam}}=100$. In the inset data
belonging to different sizes are collapsed near the critical temperature with
exponents $\alpha=0.52\pm 0.07$ and $\nu_{\rm eff}=2\beta+\gamma=1.94\pm
0.22$. As expected, due to the fact that the scaling hypothesis holds near the
critical point, the collapse of the data succeeds around the peaks of the
specific-heat curves. Moving away from the critical temperature the curves
begin to separate one from another. In the inset we show only the interesting
part of the collapse. At variance with the result found in Ref. [GAL20] for
the ML 4-phasor, the result $1/\nu_{\rm eff}=1/2$ for the scaling exponent of
the critical region is in perfect agreement with the mean-field argument we
have presented at the beginning of the chapter, in the case of a Landau
potential. Moreover, we notice that the same behavior of REM specific heat
divergence, has been found in a 4-phasor model with a random diluted topology
[GAL20]. Both these model belong to the same universality class, the one of
the standard $\phi^{4}$ Landau theory.

####  4.2 Simulation with Periodic Boundary Conditions in the Frequency Space

There are many sources of finite-size effects when simulating a dense model,
most of which can be coped with by using powerful algorithms, such as the PT
Monte Carlo algorithm. However, the ML 4-phasor model has another specific
source of finite-size effects induced by the FMC. As discussed in the previous
chapter, because of this condition, modes near the boundaries of the spectrum
($k\gtrsim 1$, $k\lesssim N$) interact much less than modes whose frequency
lays in the middle of the spectrum ($k\sim N/2$). Though their dynamic
evolution is taken into account in the simulations, the edge modes are less
and less important as the external pumping increases. To circumvent this
problem we introduce a slightly different model network, imposing periodic
boundary conditions on the frequencies [Nie+22].

In the perspective in which the FMC can be thought as a distance on the graph,
we can introduce a distance with periodic boundary conditions. This has the
effect of eliminating band-edge modes, or, equivalently, it is like
considering only modes at the center of the spectrum, as if pertaining to a
larger system. The periodic boundary conditions on the frequencies are
obtained in practice by representing the frequency indices as variables on a
ring, see Fig. 4.3, and taking their distance as the smallest one between any
two of them

$\displaystyle|k_{i}-k_{j}|=\left\\{\begin{array}[]{c c}|k_{i}-k_{j}|&\mbox{if
}~{}~{}~{}|k_{i}-k_{j}|\leq\frac{N}{2}\\\ &\\\ \frac{N}{2}-|k_{i}-k_{j}|\
\mbox{mod}\left(\frac{N}{2}\right)&\mbox{if
}~{}~{}~{}|k_{i}-k_{j}|>\frac{N}{2}\end{array}\right.$ (4.37)

The generalization to the case of quadruplets is straightforward and will be
denoted as $d_{\bm{k}}^{\text{PB}}$. In Fig. 4.3 we plot the distribution of
tetrads per mode $k$. Data in green pertain to the original mode-locked
network built by selecting modes according to $d_{\bm{k}}=0$; data in purple
pertain to the mode-locked network with periodic boundary conditions on the
frequencies, where modes are selected according to
$d_{\bm{k}}^{\text{PBC}}=0$. The distribution is approximately flat in the
case of modes selected by the FMC with periodic boundary conditions, whereas
it is peaked on the center of the spectrum in the original case.

From now on, we will refer to the version of the model with periodic boundary
conditions on the frequencies as PBC, whereas the original one, with free
boundary conditions will be termed FBC.

An important remark is that, for a certain number of modes $N$, the total
number of quadruplets $N_{4}^{*}$ satisfying the FMC with periodic boundary
conditions is slightly greater than the corresponding number in the case of
free boundary conditions, since all the modes participate in approximately the
same number of interactions. However, the order of the dilution with respect
to the fully-connected graph is the same in the two cases.

It is appropriate to anticipate that the fact that with periodic boundary
conditions the central modes are no more preferred by the interactions clearly
affects the intensity spectrum of the model, by eliminating the global
narrowing at high pumping, see e.g. Fig. 4.4. Actually, one can think that
there is no more distinction among central modes and edge modes, though the
interacting quadruplets are still selected through a deterministic condition.
We are loosing the ability of qualitatively reproducing the central narrowing
of random laser spectra, in favor of a significant reduction of finite-size
effects.

Figure 4.2: Periodic boundary conditions on the mode frequency indexes for the
frequency matching condition. The modes can be thought as distributed on a
ring, instead of a linear chain.

Figure 4.3: Number of tetrads $M$ to which each mode $k$ of the size $N=64$
participate. Green: standard FMC; purple: FMC with periodic boundary
conditions. Notice the absence of a central peak in the purple distribution.

$N$ | $N_{4}$ | $T_{\text{min}}$ | $T_{\text{max}}$ | $N_{\text{PT}}$ | $N_{\text{MCS}}$ | $N_{\text{rep}}$ | $N_{\text{s}}$
---|---|---|---|---|---|---|---
18 | $2^{8}$ | 0.35 | 1.6 | 50 | $2^{19}$ | 4 | 400
32 | $2^{11}$ | 0.45 | 1.6 | 46 | $2^{19}$ | 4 | 350
48 | $2^{13}$ | 0.5 | 1.6 | 44 | $2^{20}$ | 4 | 300
62 | $2^{14}$ | 0.55 | 1.6 | 42 | $2^{20}$ | 4 | 250
62 | $2^{14}$ | 0.3 | 1.6 | 52 | $2^{20}$ | 4 | 100
96 | $2^{16}$ | 0.65 | 1.6 | 38 | $2^{20}$ | 2 | 100
120 | $2^{17}$ | 0.7 | 1.6 | 36 | $2^{20}$ | 2 | 75
Table 4.1: Details for the simulations of the ML 4-phasor with free boundary conditions on the frequencies. Notice that for the size $N=62$ a second series of simulation has been performed with lower $T_{\text{min}}$ and $T_{\text{max}}$, in order to better explore the low temperature phase. $N$ | $N_{4}$ | $T_{\text{min}}$ | $T_{\text{max}}$ | $N_{\text{PT}}$ | $N_{\text{MCS}}$ | $N_{\text{rep}}$ | $N_{\text{s}}$
---|---|---|---|---|---|---|---
18 | $2^{9}$ | 0.05 | 1.2 | 46 | $2^{19}$ | 2 | 200
28 | $2^{11}$ | 0.1 | 1.2 | 44 | $2^{19}$ | 2 | 200
42 | $2^{13}$ | 0.2 | 1.2 | 40 | $2^{20}$ | 2 | 150
54 | $2^{14}$ | 0.25 | 1.25 | 40 | $2^{20}$ | 4 | 100
66 | $2^{15}$ | 0.25 | 1.25 | 40 | $2^{20}$ | 2 | 100
82 | $2^{16}$ | 0.3 | 1.3 | 40 | $2^{20}$ | 2 | 100
104 | $2^{17}$ | 0.35 | 1.3 | 38 | $2^{21}$ | 2 | 80
Table 4.2: Details for the simulations of the ML 4-phasor model with periodic
boundary conditions on the frequencies.

#####  4.2.1 Details of the Simulation

We have performed numerical simulations of the ML 4-phasor model, both with
free and with periodic boundary conditions on the frequencies.

First, the number of interacting quadruplets $N_{4}$ is chosen as a power of
$2$. Then, the corresponding size $N$ is selected in such a way that the
difference between $N_{4}$ and the true number of tetrads $N_{4}^{*}$ is
minimum, as explained in the previous chapter. Since $N_{4}^{*}$ is greater
with PBC rather than with FBC, we managed to perform simulations in the PBC
case with networks which are larger from the point of view of the number of
couplings, though with smaller sizes. In both cases, however, the simulated
network with the highest number of quadruplets has $N_{4}=2^{17}$.

For each size $N$ of the simulated systems, we have run PT simulations with
$N_{\text{PT}}$ replicas of the system at temperatures $T_{i}\in[T_{\rm
min},T_{\rm max}]$, with $i=1,\dots,N_{\rm PT}$. On top of that, we have
performed simulations on $N_{\mbox{rep}}$ identical replicas of the PT
simulations in order to compute overlap distributions. The temperature
interval has been chosen properly in order to have a sufficient number of
replicas above the size-dependent critical temperature $T_{c}(N)$ and ensure
faster thermalization of the replicas in the low temperature glassy phase. We
have chosen equispaced temperatures with spacing $\Delta T=0.025$ and a number
of Monte Carlo steps $N_{\text{swap}}=64$ after which a swap of configurations
between adjacent heat baths is proposed. Both these choices have been checked
to be compatible with high acceptance ratios for the swaps for the whole
duration of the simulation. Each copy of the system at each temperature share
the same realization of quenched disordered couplings $\\{J_{\bm{k}}\\}$. The
number of simulated disordered samples has been chosen in relation with the
size of the system. The values of the simulation parameters are reported in
Table 4.1 for the model with FBC and Table 4.2 for the model with PBC.

####  4.3 Numerical Results

We devote this section to the results of the Monte Carlo analysis of the ML
4-phasor model Eq. (3.1) with FBC and PBC on the frequencies.

#####  4.3.1 Spectra

Figure 4.4: (a,c) Intensity spectrum for a single realization of quenched
disorder of the ML $4$-phasor model with free (periodic) boundary conditions
on the frequencies and $N=120$ ($N=104$) modes. Color map: temperature
$T\in[0.7,1.45]$ ($T\in[0.35,1.1]$). (b,d) Intensity spectrum with free
(periodic) boundary conditions on the frequencies and $N=120$ ($N=104$)
averaged over $N_{\text{s}}=75$ ($N_{\text{s}}=80$) instances of quenched
disorder. All spectra are normalized to their area and modes $k$ are divided
by $N$.

The first observable we display is the intensity spectrum Eq. (3.16). Let us,
first, briefly comment on the relationship between the physical intensities
$I_{k}$ and the complex amplitude variables of the simulated model (3.1). In
real experiments the heat bath temperature $T$ is typically kept fixed (there
are exceptions like, e.g., in Ref. [WC01]) and the overall system energy
${\mathcal{E}=\epsilon N}$ is varied by tuning the pumping power. As already
discussed (Chap. 3), in our simulations, $\epsilon$ is fixed and kept equal to
one in the spherical constraint, $\sum_{k}|a_{k}|^{2}=\sum_{k}A_{k}^{2}=N$,
whereas $T$ is varied. Therefore, according to
$\mathcal{P}=\epsilon/\sqrt{T}$, a change in the pumping rate $\mathcal{P}$
because of a shift in the energy $\epsilon$ pumped into the system corresponds
to a shift of $1/\sqrt{T}$. If we rescale the intensity of the mode $k$ as in
Eq. (3.16), i.e. $I_{k}=A_{k}^{2}/\sqrt{T}$, we have
$\sum_{k}I_{k}=N/\sqrt{T}=N\epsilon$, as in Eq. (2.41).

In Fig. 4.4, we show the emission spectra at equilibrium for the ML $4$-phasor
model with FBC and PBC, both for a single instance of disorder (panels (a,c))
and averaged over roughly a hundred instances of disorder (panels (b,d)). The
most relevant difference between the model with FBC and PBC is the complete
absence of global narrowing in the spectrum of the PBC case, which corresponds
to the absence of band-edge modes: all modes interact with identical
probability with the rest of the system.

On the other hand in the FBC case, one can observe the typical central
narrowing occurring in random lasers [Cao+99, Cao+00, Wie08, Gho+14] as the
pumping energy increases. This phenomenon becomes particularly evident in the
averaged spectrum (panel (b)), which is smoother than the single sample one.
Finally, we notice that the averaged spectrum in the PBC case looks like the
central part of the FBC spectrum.

One of the most relevant features of all the intensity spectra shown in Fig.
4.4 is that they become more and more structured and heterogeneous upon
decreasing the temperature. The pattern of the peaks is disordered and
strongly depends on the random sample and on the single dynamic history. A
first analysis of this phenomenon in terms of intensity equipartition breaking
among the different modes has been performed in [GAL20], and a deepening of
the collective inhomogeneous behavior of the modes will be presented in Chap.
5.

#####  4.3.2 Specific Heat

In Figs. 4.6 and 4.6 we show the specific heat respectively for the ML
4-phasor model with FBC and PBC. In the main panel data are plotted with
respect to temperature: each point corresponds to the equilibrium fluctuations
of energy within a certain heat bath averaged over the disordered samples. In
the inset, data are collapsed according to the scaling hypothesis previously
discussed. The finite-size scaling analysis of the specific heat has been
performed in a more refined way with respect to the simple case of the REM. In
order to get the two exponents $\alpha$ and $\nu_{\rm eff}=2\beta+\gamma$ of
Eq. (4.33) from our numerical data we follow the method proposed in Refs.
[al13, Leu+15].

Figure 4.5: Specific heat $c_{V_{N}}$ for the ML 4-phasor model with free
boundary conditions on the frequencies as a function of $T$. Different curves
represent different simulated sizes of the system. The simulated sizes are
$N=18,32,48,62,96,120$. (Inset) Specific heat scaled by $N^{\alpha/\nu_{\rm
eff}}$ as a function of $\tau N^{1/\nu_{\rm eff}}$, with $\alpha=0.48$,
$\nu_{\rm eff}=0.9$.

Figure 4.6: Specific heat $c_{V_{N}}$ for the ML $4$-phasor model with
periodic boundary conditions on the frequencies as a function of $T$.
Different curves represent different simulated sizes of the system. The
simulated sizes are $N=18,28,42,54,66,104$. (Inset) Specific heat scaled by
$N^{\alpha/\nu_{\rm eff}}$ as a function of $\tau N^{1/\nu_{\rm eff}}$, with
$\alpha=0.27$, $\nu_{\rm eff}=1.16$.

First, the size-dependent critical temperatures $T_{c}(N)$ are more precisely
assessed by fitting the points around the peak of each curve in the main
panels of Figs. 4.6 and 4.6 with a quadratic function of the temperature
$f_{N}(T)=a_{N}+b_{N}T+c_{N}T^{2}$. The critical temperatures are identified
with the maximum of each of the fitting functions $T_{c}(N)=-b_{N}/(2c_{N})$,
with a statistical error estimated accordingly. The results of this procedure
are reported in Table 4.3. The critical temperature $T_{c}(\infty)$ of the
models can be extrapolated from the fit of the finite-size critical
temperatures with the following function: $T_{c}(N)=T_{c}(\infty)+aN^{-b}$,
where the exponent $b$ gives a first rough estimate of the critical exponent
$1/\nu_{\rm eff}$. The results of the fit are:

FBC: $\displaystyle T_{c}(\infty)=0.86\pm 0.03,\quad b=1.6\pm 0.5,$ (4.38)
PBC: $\displaystyle T_{c}(\infty)=0.61\pm 0.03,\quad b=0.98\pm 0.3.$ (4.39)

As one can see, already from this rough estimate the model with PBC has an
exponent $1/\nu_{\rm eff}$ falling inside the interval derived for generic
mean-field models (4.14) at the beginning of this chapter. In comparison to
the estimate provided from FBC data, indeed, there appears to be a drastic
reduction of finite-size effects.

| FBC | PBC
---|---|---
$N_{4}$ | $N$ | $T_{c}$ | $\Delta T_{c}$ | $N$ | $T_{c}$ | $\Delta T_{c}$
$2^{8}$ | 18 | 0.55 | 0.04 | - | - | -
$2^{9}$ | - | - | - | 18 | 0.42 | 0.02
$2^{11}$ | 32 | 0.63 | 0.025 | 28 | 0.49 | 0.02
$2^{13}$ | 48 | 0.69 | 0.02 | 42 | 0.52 | 0.02
$2^{14}$ | 62 | 0.75 | 0.03 | 54 | 0.55 | 0.03
$2^{15}$ | - | - | - | 66 | 0.56 | 0.04
$2^{16}$ | 96 | 0.8 | 0.07 | 82 | 0.56 | 0.05
$2^{17}$ | 120 | 0.83 | 0.09 | - | - | -
Table 4.3: Values of the critical temperatures for each simulated size of the
ML 4-phasor model with fixed and periodic boundary conditions on the
frequencies.

Then we take the following Ansatz on the form of the scaling function
$\hat{f}$ in Eq. (4.33)

$\displaystyle\hat{f}(x)=A+Cx^{2},$ (4.40)

where $x=N^{1/\nu_{\rm eff}}t_{N}$, with $t_{N}$ computed by using the
$T_{c}(N)$ reported in Table 4.3. In the previous Ansatz we have not included
the linear term, since the points are translated in order for the peak of each
curve to be in the origin and we expect the linear term not to matter. With
this Ansatz the scaling hypothesis for the specific heat Eq. (4.33) reads as

$\displaystyle c_{V_{N}}(T)=\tilde{A}_{N}+\tilde{C}_{N}t_{N}^{2},$ (4.41)

where $\tilde{X}_{N}=X_{N}N^{\frac{\alpha+x}{\nu_{\rm eff}}}$, with
$X_{N}=\\{A_{N},C_{N}\\}$ and $x=\\{0,2\\}$. For each size we select a
temperature interval centered in $T_{c}(N)$ corresponding to the points
plotted in the insets of Figs. 4.6 and 4.6. We fit the points in the selected
interval with the previous function and determine the values of the
coefficients. We notice that the behaviour of the logarithm of the
coefficients absolute value, i.e.

$\displaystyle\log|\tilde{X}_{N}|=\log|X_{N}|+\frac{\alpha+x}{\nu_{\rm
eff}}\log N$ (4.42)

is linear in $\ln N$; hence, the estimates of $\alpha$ and $\nu_{\rm eff}$ can
be obtained by linear interpolation.

Figure 4.7: Overlap distributions for five instances of disorder and for the
average over all instances at $T=0.25\simeq 0.45~{}T_{c}$. Simulation size
$N=54$ of the ML 4-phasor model with periodic boundary conditions on the
frequencies.

For the systems with FBC, this finite-size analysis provides the following
values for the critical exponents

$\displaystyle\mbox{FBC:}~{}~{}\alpha=0.48\pm 0.05,\quad 1/\nu_{\rm
eff}=1.1\pm 0.16.$ (4.43)

With respect to the estimate $1/\nu_{\rm eff}\simeq 3/2$ found in [GAL20], a
much larger statistics allows to find an estimate of $1/\nu_{\rm eff}$ closer
to, and compatible with, the mean-field threshold, suggesting that deviations
from mean-field theory might be due to pre-asymptotic effects in $N$. The
confirmation that this is, indeed, the origin of the anomalous value
previously found for $1/\nu_{\rm eff}$ comes from the analysis with PBC. In
this case, the critical exponents turn out to be

$\displaystyle\mbox{PBC:}~{}~{}\alpha=0.27\pm 0.05,\quad 1/\nu_{\rm
eff}=0.86\pm 0.14.$ (4.44)

With PBC we find an estimate of $1/\nu_{\rm eff}$ well below the threshold for
a mean-field universality class. Therefore, up to the precision of our
analysis, despite being possibly still of a different universality class with
respect to the REM, for which $1/\nu_{\rm eff}=1/2$, we can assess the mean-
field nature of the glass transition in the ML 4-phasor model.

#####  4.3.3 Overlap Distribution Functions

In this section we complete the study of the glass transition in the ML
4-phasor model, by presenting results of our numerical simulations about the
first-order nature of the transition.

###### Parisi overlap

Figure 4.8: Parisi overlap distribution for the size $N=62$ of the ML 4-phasor
model with free boundary conditions on the frequencies. The distribution is
averaged over $N_{\text{s}}=100$ instances of disorder. Color map: same as
Fig. 4.4. The blue curve corresponding to the lowest temperature is at
$T\simeq 0.4T_{c}$, with $T_{c}=0.86(3)$.

Figure 4.9: Parisi overlap distribution for the size $N=54$ of the ML 4-phasor
model with with periodic boundary conditions on the frequencies. The
distribution is averaged over $N_{\text{s}}=100$ instances of disorder. Color
map: same as Fig. 4.4. The blue curve corresponding to the lowest temperature
is at $T\simeq 0.45T_{c}$, with $T_{c}=0.61(3)$.

Let us first discuss the case of the Parisi overlap distribution (• ‣ 3.2.4).
The protocol used in numerical simulations to measure overlaps corresponds to
the definition of replicas as independent copies of the system with the same
quenched disorder. For each sample, i.e. each realization of disorder, we run
dynamics independently for $N_{\text{rep}}$ replicas of the system, starting
from randomly chosen initial phasor configurations. In this way, replicas
explore different regions of the same phase space, passing through
configurations typically belonging to separated equilibrium states, if there
are many of them, and, sometimes, to the same state.

To study the behavior of the $P_{J}(q)$ we choose $N_{\text{rep}}=4$, so that
at any measurement time six values of the overlap are available
$q_{\alpha\beta}=\\{q_{01},q_{02},q_{03},q_{12},q_{13},q_{23}\\}$. Hence,
passing from two replicas to four increases the statistics by a factor six. In
order to accumulate statistics, we measure the value of $q_{\alpha\beta}$
using $\mathcal{N}$ equilibrium, time uncorrelated, configurations of replicas
at the same iteration of the simulated dynamics, see Eq. (3.13). Hence, for
each disordered sample the $P_{J}(q)$ histograms are built with
$\mathcal{N}\times N_{\text{rep}}(N_{\text{rep}}-1)/2$ values of the overlap.

The overlap distribution functions $P_{J}(q)$ are computed as the normalized
histograms of the overlaps for each one of the samples. This has been done for
each simulated size of the ML 4-phasor model with both FBC and PBC. In Fig.
4.7 we present the overlap distributions for five samples at the temperature
$T=0.25\simeq 0.45T_{c}$ of the size $N=54$ of the ML 4-phasor model with PBC,
together with the overlap distribution averaged over $100$ samples. Given the
fluctuations of $P_{J}(q)$ among the different samples, it is clear that the
only physical quantity to be considered in order to assess the glass
transition is the averaged $P(q)\equiv{\overline{P_{J}(q)}}$.

This is particularly important in the case of the overlap distribution
function, since, it is not a self-averaging quantity [MPV87], i.e., the
average $P(q)$ cannot be reached simply by increasing the size of the system
over which a single sample $P_{J}(q)$ is built, as, for instance in the case
of the free energy or of the specific heat, but only by averaging over
disorder.

In Fig. 4.9 and Fig. 4.9 the average overlap distribution function of the ML
4-phasor model with FBC and PBC are, respectively, reported for the whole
simulated temperature range in systems whose size correspond to
$N_{4}=2^{14}$, i.e. respectively $N=62$ and $N=54$ spins. The reduction of
the finite-size effects obtained by using periodic boundary conditions in the
choice of interacting modes leads to display $P(q)$ with more distinct
secondary peaks in the case of the ML 4-phasor model with PBC.

We have showed the overlap distribution functions for the $N=62$ (for the
model with FBC) and $N=54$ (for the model with PBC), because these are the
largest sizes for which we performed simulations with $N_{\text{rep}}=4$.
However, it can be useful to show also the behavior of the $P(q)$ for a higher
size. In Fig. 4.10 we show the overlap distribution function for the size
$N=82$ of the model with periodic boundary conditions. One can see that the
side peaks of the distribution at low temperatures are not appreciably more
pronounced than than those of the distribution for $N=54$. However, a
reduction of the finite size effects with respect to $N=54$ can be appreciated
at high temperature, where the distribution is still a Gaussian, but it
slightly narrows around the peak in $q=0$. We expect that in the large $N$
limit the overlap distribution function at high temperature becomes a delta
function peaked in $q=0$.

Figure 4.10: Parisi overlap distribution for the size $N=54$ of the ML
4-phasor model with with periodic boundary conditions on the frequencies. The
distribution is averaged over $N_{\text{s}}=100$ instances of disorder. Color
map: same as Fig. 4.4.

###### Plaquette overlap and IFO

The study of the plaquette overlap and IFO distribution functions has been
performed both for the model with FBC and PBC. In the following we only
display results for the PBC case, which is less affected by finite-size
effects.

In Fig. 4.11 the plaquette overlap distribution is plotted for all the
simulated temperatures and for the system size $N=54$. For each of the
$\mathcal{N}$ uncorrelated configurations at equilibrium, the plaquette
overlaps are computed over $N_{4}\sim O(N^{3})$ quadruplets, leading to a
reduction of the statistical error on the overlap values. Then, for each
sample, the distribution is computed with $\mathcal{N}\times
N_{\text{rep}}(N_{\text{rep}}-1)/2$ values of the plaquette overlap. Data in
figure are averaged over $N_{\rm s}=100$ disordered samples.

In the low temperature region, the distribution clearly develops a nontrivial
shape, with a heavy tail for values of $\mathcal{Q}>0$. The presence of the
three visible peaks in the distribution tail at the lowest temperatures is a
consequence of the single sample behavior of the $P(\mathcal{Q})$, which we
present in Fig. 4.12 for the simulated temperature $T=0.25$. As one can see,
the shape of the distribution dramatically changes from sample to sample: in
particular, for some instances of disorder a single peak, for others two peaks
are displayed in the $P(\mathcal{Q})$, at sample dependent positions. The
effect of taking the average over all the simulated samples is reported in the
sixth panel in Fig. 4.12, which corresponds to the lowest $T$ curve in Fig.
4.11. We notice that, even with $N_{s}=100$ disordered samples, the average is
far from a smooth function, displaying the strong lack of self-averaging of
the plaquette distribution.

In Fig. 4.14 the averaged IFO distribution is displayed for the system size
$N=54$, at all simulated temperatures. Even with the significant reduction of
finite-size effects obtained through the PBC, the distribution obtained with
our data at equilibrium does not show clear side peaks in the low temperature
region. Actually, side peaks are more evident very close to the transition
temperature (gray curves in figure), rather than for the lowest simulated
temperatures. Once again, we get a clear understanding of the averaged
distribution by looking at single instances of disorder. In Fig. 4.14, one can
see that, while many samples have a distribution which is still peaked in
$\mathcal{C}=0$, a signature that they may not have entered in the glassy
phase yet, in many others the $P(\mathcal{C})$ has developed nice side-peaks
at the temperature $T=0.25$. However, when taking the average, the combined
effect of samples of the first kind and of the variable position of the peaks
in samples of the second kind, leads to the smoothed and almost flat shape of
the IFO distribution at low temperature.

Figure 4.11: Plaquette overlap distribution for the size $N=54$ of the ML
4-phasor model with periodic boundary conditions on the frequencies. The
distribution is averaged over $N_{\text{s}}=100$ instances of disorder. Color
map: same as Fig. 4.4. The blue curve corresponding to the lowest temperature
is at $T\simeq 0.45T_{c}$, with $T_{c}=0.61(3)$. Figure 4.12: Plaquette
overlap distributions for five instances of disorder and for the average over
all instances at $T=0.25\simeq 0.45~{}T_{c}$. Simulation size $N=54$ of the ML
4-phasor model with periodic boundary conditions on the frequencies. As for
the Parisi overlap distribution, each sample exhibits a particular
distribution shape, with one or two peaks appearing at different positions:
the relevant quantity in the thermodynamic limit is the averaged distribution
in the sixth panel. Figure 4.13: IFO distribution for the size $N=54$ of the
ML 4-phasor model with with periodic boundary conditions on the frequencies.
The distribution is averaged over $N_{\text{s}}=100$ instances of disorder.
Color map: same as Fig. 4.4. The blue curve corresponding to the lowest
temperature is at $T\simeq 0.45T_{c}$, with $T_{c}=0.61(3)$. Figure 4.14: IFO
distributions for five instances of disorder and for the average over all
instances at $T=0.25\simeq 0.45~{}T_{c}$. Simulation size $N=54$ of the ML
4-phasor model with periodic boundary conditions on the frequencies.

###  Chapter 5 Breaking of Equipartition and Pseudo-Localization at the
Transition

The present chapter is devoted to the study of a particular phenomenon taking
place in the ML 4-phasor model, which has been mentioned previously in this
work, when describing the intensity emission spectra of random lasers:
intensity localization of light, else termed _power condensation_ [ACL15a,
AIL15a], and its relationship to the high pumping replica symmetry breaking
(RSB) random lasing phase.

Localization is a widespread phenomenon in physics, which has always been
related to the breaking of ergodicity. In the pioneering work of Anderson on
semiconductors, the spatial localization of the electron wavefunction induced
by a large degree of disorder was identified as the underlying mechanism for a
semiconductor/insulator transition [And58]. The absence of diffusion in the
insulating-localized phase was interpreted as a clear manifestation of
dynamical ergodicity breaking. Besides the seminal work of Anderson,
localization was related to ergodicity breaking also in classical systems,
since the famous numerical study of Fermi-Pasta-Ulam-Tsingou (FPUT) on the
anharmonic chain [Fer+55]. In this case, localization was observed in the
Fourier space of the chain modes. Starting from an atypical initial condition,
e.g. only the lowest harmonic excited, one would expect that a slight
anharmonicity is sufficient to cause the system to relax on a state where the
energy is equally divided among all the modes111If one initializes a harmonic
chain on a eigenmode of the Hamiltonian, the time evolution will leave the
system on that eigenmode, resulting in a breaking of ergoditicy. Nonlinearity
is believed to facilitate the recovery of ergodicity, because it introduces a
coupling which makes the eigenmodes of the unperturbed Hamiltonian more
connected.. Contrary to expectations, the system showed a recurrent dynamics
for all the duration of the experiment, with no sign of relaxation to
equipartition. Therefore, localization was coming along with dynamical
ergodicity breaking also in the case of energy localization in the Fourier
power spectrum.

Ergodicity breaking was understood as a purely dynamical phenomenon until the
’80s, when the theory of replica symmetry breaking was established as a new
thermodynamic paradigm and statistical ensembles formalism for ergodicity
breaking transitions in complex disordered systems. Quite interestingly, while
localization phenomena in quantum many-body systems have been widely
investigated during the last decades [Alt+97, BAA06, NH15, RMS15, VR16, AP17,
AL18], they have seldom been probed in disordered systems, apart from the
attempts of Refs. [MS18, GAL20]. The lack of a broad analysis of localization
phenomena in the context of disordered glassy systems is mainly due to the
nature of the variables which are customary for these systems: Ising, XY or
Heisenberg spins, all locally bounded, $|\vec{s}|=1$. In those models, such as
the spherical $p$-spin model, where variables are used with continuous locally
unbounded magnitude as a proxy for magnetic spins in spin-glasses or density
fluctuations in structural glasses, in order to be able to perform analytical
computations, the interaction network is usually fully connected, thus
hindering any sort of magnitude localization. Indeed, this kind of mean-field
representation on a complete graph, together with a local potential (soft
spins) or a global constraint (spherical spins) guarantees magnitude
equipartition.

Therefore a careful investigation of how magnitude localization coexists with
replica symmetry breaking in disordered systems is a gap that needs to be
filled. Once again, the spin-glass theory of random lasers represents a very
fortunate research field for this kind of analysis. The fundamental variables,
i.e. the amplitudes of the light modes, are naturally continuous and locally
unbounded; the global constraint which they are subjected to is a quite
natural requirement for the stationary regime of a lasing system (see Chap.
2); the presence of dilution is a direct consequence of the specific selection
rule in light mode coupling, i.e. the FMC. Hence, these systems have all the
ingredients necessary to exhibit a power condensation transition.

As a disclaimer, it is worth stressing that this is not the generalization to
light waves of the spatial wavefunction localization occurring in Anderson
theory, that is known to be inhibited in 3D random lasers because of the
vectorial nature of light waves [SS14, Spe+16]. Rather, it is a condensation
of the overall magnitude on a few variables in a set of locally unbounded
variables subjected to a global constraint, i.e. a global conservation law.
This kind of condensation has already been observed and studied analytically
in non-interacting systems [MEZ05, EMZ06, TS01]. When the value of this
globally conserved quantity exceeds a given (non-universal) threshold the
system undergoes a transition where a macroscopic fraction of the conserved
quantity concentrates on a finite portion of the system. Given the
clarification, from now on we will use the terms ‘localization’ and
‘condensation’ interchangeably to refer to the phenomenon of our interest.

This chapter is organized as follows: first, a few details are provided
regarding the condensation transition in non-interacting systems, pointing out
the emergence of a pseudo-localized phase, which has been discovered in Refs.
[Gra+21a, Gra+21]; then we clarify to what extent, in the presence of an
interaction network, the order of dilution is key to understand possible
regimes of equipartition or localization of a globally conserved quantity;
finally, results from the numerical simulations of the ML 4-phasor model are
presented, revealing the presence of a _hybrid_ phase analogous to the pseudo-
localized phase of non-interacting systems, where the intensity of light modes
is neither equipartitioned among all modes nor really localized on a few of
them.

####  5.1 Pseudo-Localization in Non-Interacting Systems

Condensation of a global quantity on a finite number of degrees of freedom has
been found and very precisely described in the framework of large deviation
calculations and ensemble inequivalence in the case of mass-transport models
[MEZ05, EMZ06] or for bosonic condensates in optical lattices described by the
Discrete Non-Linear Schrödinger Equation (DNLSE) [TS01]. The fundamental
ingredient involved in the localization phenomenology for non-interacting
systems is the presence of a global constraint.

Let us consider, for instance, the DNLSE. We focus on the infinite temperature
limit in which the hopping term, i.e. the kinetic two-body term in the
Hamilotnian of the DNLSE, can be neglected. In this case, the partition
function of the model reads

$\displaystyle\Omega(\mu,E)=\int\prod_{i=1}^{N}\text{d}\psi_{i}\text{d}\bar{\psi}_{i}\
e^{-\mu\sum_{i=1}^{N}|\psi_{i}|^{2}}~{}\delta\left(E-\sum_{i=1}^{N}|\psi_{i}|^{4}\right),$
(5.1)

where $\mu$ is the chemical potential. Clearly the quantity

$\displaystyle A=\sum_{i=1}^{N}|\psi_{i}|^{2},$ (5.2)

represents the mass of the condensate, so that the partition function in Eq.
(5.1) corresponds to an ensemble where exact conservation of energy is
enforced by means of a Dirac delta function, whereas mass is conserved only on
average by means of a field $\mu$. In the following, we will refer to the
former way of imposing a global conservation law as “hard” constraint and to
the latter as “soft” constraint.

In the case of the DNLSE, as in all cases where it takes place, the physical
quantity that localizes is the one controlled by the hard constraint, hence in
this case it is the energy. It is only thanks to the global action of the
constraint on the total energy that configurations with a strongly
heterogeneous distribution of energy on lattice sites are allowed. The
analytical calculations of Refs. [Gra+21a, Gra+21] show that as soon as energy
is constrained above a certain critical value, $E>E_{c}$, these localized
configurations dominate the partition function. It can be shown analytically,
but it can be easily guessed by looking at Eq. (5.1), that localization cannot
take place for a quantity controlled “on average”. It is not possible to have
strongly inhomogeneous fluctuations and/or localization of something which is
controlled homogeously by means of a Lagrange multiplier like the chemical
potential in Eq. (5.1). This is precisely the same mechanism characterizing
the Bose-Einsten condensation, which is a form of localization in Fourier
space: the condensed phase cannot be reached by controlling density with a
chemical potential; density must be tuned directly, for instance, by
decreasing the volume for a given number of bosons [Hua87]. The fact that the
Bose-Einstein condensation cannot be implemented by tuning the chemical
potential of a reservoir in contact with the system is analogous to the fact
that energy localization cannot be achieved by tuning the temperature of a
thermostat, i.e. by studying the partition function where conservation of
energy is imposed “on average”, as $\exp(-\beta\mathcal{H})$, rather than
exactly, as $\delta(E-\mathcal{H})$. In both the above examples localization
entails lack of statistical ensemble equivalence: for the energy localization
in the DNLSE it is the lack of equivalence between fixed temperature
(canonical) and fixed energy (microcanonical) ensembles, for the Bose-Einstein
condensation it is the lack of equivalence between fixed chemical potential
(grancanonical) and fixed density (canonical) ensembles.

One important feature that has been predicted for the DNLSE [Gra+21a, Gra+21]
is the existence of a pseudo-localized phase similar to what we are going to
show in numerical simulations of the ML 4-phasor model. The peculiarity of the
localization phenomenon in the case of the DNLSE is that it takes place in two
steps. First, by increasing the energy, i.e. the quantity which controls
localization, one first encounters a second-order transition at a value of the
energy – $E_{\rm th}$ – where the equivalence of ensembles breaks down and
temperature becomes negative. Then, at a larger value of energy,
$E_{c}>E_{th}$, there is a first-order transition to a localized phase. For
$E\in[E_{\rm th},E_{c}]$ the system finds itself in a pseudo-localized phase,
where a thermodynamic anomaly is indicated by the lack of ensemble equivalence
and the presence of negative temperature, but localization has not been really
achieved yet.

The main difference between the ML $4$-phasor model of the glassy random laser
and DNLSE, is that in the random laser case the joint distribution of the
variables over which the global constraint is imposed is not factorized.
Therefore, the analytical results discussed in [Gra+21a, Gra+21] cannot be
straightforwardly extended to this case. In the following, first we try to
derive some understanding of localization in generic interacting systems;
then, we resort to numerical simulations in order to probe the localization
transition in the ML 4-phasor model.

####  5.2 Scaling Argument for the Occurrence of Intensity Localization in
Interacting Systems

In this section we go through a scaling argument on generic $p$-spin
interacting systems (both with ordered and disordered interactions) with
continuous variables and a global constraint, of which the ML 4-phasor model
with the spherical constraint (2.41) is a special case.

Let us consider the $p$-spin Hamiltonian

$\mathcal{H}[\sigma]=-\sum_{k_{1}\ldots k_{p}}^{\\#N^{A}}J_{k_{1}\ldots
k_{p}}\ \sigma_{k_{1}}\ldots\sigma_{k_{p}},$ (5.3)

whose $N$ continuous spin variables $\sigma$ are subjected to a generic
$\rho$-metrical constraint (e.g., $\rho=2$ is the case of the spherical
constraint)

$\sum_{i=1}^{N}\sigma_{i}^{\rho}=N,$ (5.4)

where $N^{A}$ on top of the sum, with $A\in[1,p]$, denotes the scaling with
the size of the number of $p$-uples contributing to the energy and the spin
indices $k_{i}$ run from $1$ to $N$. If $A=p$ we have a fully-connected
interaction graph, i.e. each spin contributes in $\mathcal{(}N^{p-1})$
$p$-uples. At the other extreme, if $A=1$ the graph is sparse, i.e. each spin
only interacts in a finite number of $p$-uples, not growing with the size of
the system. All dilutions in between will be considered hereafter.

For simplicity, we will take the variables as real valued here, yet keeping
the word intensity for the spin magnitude $|\sigma|$. We notice that the
glassy random laser is a model belonging to this family, with $p=4$ but with
complex spins. Though complex variables yield new physical features (see Chap.
2), we stress that these are not relevant for what concerns the influence of
connectivity on the possible onset of intensity localization.

#####  5.2.1 Ordered Couplings

First we consider the ordered case $J_{k_{1}\ldots k_{p}}=J$. The typical
ground state spin configurations in the canonical ensemble, where
equipartition is expected to hold (in the average, not strictly), are those
minimizing Eq. (5.3). The energy is extensive

$E=\mathcal{H}[\sigma_{\rm gs}]=\mathcal{O}(N)$

provided that the coupling constant scales as

$J\propto\frac{1}{N^{A-1}}.$ (5.5)

In the occurrence of intensity localization, that is, if only a few modes take
the overall intensity, equal to $N$ according to Eq. (5.4), whereas all the
other are zero, we are interested in the energy contribution of a localized
spin configuration. First of all, let us notice that in order to have a non-
zero contribution the intensity at least $p$ coupled spins must localize. If
we represent by $\square$ such a localizing $p$-uple, the intensity localized
configuration is

$\\{\sigma_{\rm loc}\\}:\qquad\sigma_{k\in\square}\propto
N^{1/\rho},\qquad\sigma_{k\notin\square}=0.$ (5.6)

According to Eqs. (5.3), (5.5) the energy of such a configuration of spins
scales with $N$ like

$E_{\rm loc}=\mathcal{H}[\sigma_{\rm
loc}]=\mathcal{O}\left(\frac{N^{p/\rho}}{N^{A-1}}\right).$

To figure out whether intensity condensation might occur and dominate, one
eventually has to compare the scaling behaviors of the energies of an
equipartitioned and a localized configuration:

$\mathcal{O}(N)\quad\mbox{vs}\quad\mathcal{O}(N^{\frac{p}{\rho}+1-A}).$

Hence, we notice that the kind of global constraint imposed is also key to
understand whether a localization transition may take place. The following
cases occur for $\rho=2$, depending on the interaction connectivity scaling
$N^{A}$:

  1. 1.

$A>p/2$. Any possible intensity localized configuration of spins would yield
subextensive contributions to the energy. The equipartition regime is,
therefore, dominant. The case $A=p$ is the fully connected interaction graph.

  2. 2.

$A=p/2$. Both kinds of spin configurations yield an $\mathcal{O}(N)$
contribution to the energy. In this case a pseudo-localized phase might occur.

  3. 3.

$A<p/2$. Intensity localization provides the most prominent contribution to
the energy, that is, $\mathcal{O}(N^{>1})$. The case $A=1$ is the sparse case.

Figure 5.1: Schematic representation of the possible regimes of low
temperature intensity equipartition, intensity localization or pseudo-
localization with respect to the coupling dilution in the $4$-spin spherical
model, both for the ordered unfrustrated version (left) and the quenched
disordered, frustrated one (right). The case of the mode-locked random laser
whose connectivity $c\sim N^{3}$ is pointed by the arrow. Here, we have used
the phasor notation $a$, instead of $\sigma$, used in the text for real
variables.

#####  5.2.2 Quenched Disordered Couplings

If the interaction couplings $J_{k_{1}\ldots k_{p}}$ are quenched disordered,
independently distributed and with zero mean ${\overline{J_{k_{1}\ldots
k_{p}}}}=0$, the typical ground state of the Hamiltonian (5.3) is extensive –
$E=\mathcal{H}[\sigma_{\rm gs}]=\mathcal{O}(N)$ – if the variance of the
distribution of the couplings scales like

${\overline{J^{2}_{k_{1}\ldots k_{p}}}}\propto\frac{1}{N^{A-1}}.$ (5.7)

If the total intensity of the system is localized in a single interacting
$p$-uple, as in (5.6), Eqs. (5.3), (5.7) imply that the energy scales with the
size like

$E_{\rm loc}=\mathcal{H}[\sigma_{\rm
loc}]=\mathcal{O}\left(\frac{N^{p/\rho}}{N^{(A-1)/2}}\right).$

Comparing the equipartitioned contribution $E=\mathcal{O}(N)$ and the
localized contribution $E_{\rm loc}$ we find for $\rho=2$ the following three
regimes as the exponent $A$ varies:

  1. 1.

$A=p$. Localized energy contibutions are subextensive. The equipartition
regime is dominant. This is the fully connected interaction graph case.

  2. 2.

$A=p-1$. Both kinds of spin configurations yield an $\mathcal{O}(N)$
contribution to the energy. This is the case, e.g., of the mode-locked glassy
random laser, where $p=4$, $A=3$. In this case one might conjecture the
occurrence of a pseudo-localized phase.

  3. 3.

$A<p-1$. Intensity localization provides superextensive contributions to the
energy.

In figure (5.1) we summarize the scaling argument predictions on a pictorial
diagram with the known cases of equipartition, localization and pseudo-
localization for the $4$-spin model.

####  5.3 Evidence of Pseudo-Localization in the ML 4-Phasor Model

Let us turn back to our system. Besides the presence of interactions, there is
another important difference between the glassy random laser and the case of
the DNLSE. In the partition function of the ML 4-phasor model, which reads as

$Z_{N}(\beta,\mathcal{E})=\int\prod_{k=1}^{N}\text{d}a_{k}\text{d}\bar{a}_{k}\
e^{-\beta\mathcal{H}[\bm{a}]}\delta\left(\mathcal{E}-\sum_{k=1}^{N}|a_{k}|^{2}\right),$
(5.8)

with $\mathcal{H}$ given by Eq. (3.1), the conservation of energy is realized
on average, by imposing homogeneously a temperature for all interacting
quadruplets, while the conservation of the total intensity is realized
exactly, by means of the hard global constraint (2.41). It is then possible to
guess that in the ML $4$-phasor model a localization transition might occur at
the level of intensity, rather than energy. If that were the case, though, the
algorithms used in our numerical simulations would not be suitable anymore. An
intensity localized system, indeed, would display a small, fixed number of
variables whose energy fluctuations are no longer dependent on temperature.

In a way analogous to the condensation in the DNLSE, intensity localization is
achieved by tuning the physical quantity $\mathcal{E}=\epsilon N$ controlled
by the overall hard constraint (2.41). More precisely, when a certain
threshold is overcome in the controlling parameter of the constraint,
$\mathcal{E}>\mathcal{E}_{c}$, configurations for which a finite amount of the
overall intensity is stored in $\mathcal{O}(1)$ quadruplets might become
thermodynamically dominant. However, as already discussed in Chap. 3,
numerical simulations are performed by keeping the optical power fixed
($\epsilon=1$) and varying the temperature $T$ or the pumping rate
$\mathcal{P}$. We notice that this is equivalent to sample configurations from
the equilibrium distribution

$P[\bm{\hat{a}}]\propto
e^{-\mathcal{H}[\bm{\hat{a}}]}\delta\left(\mathcal{P}N-\sum_{k=1}^{N}|\hat{a}_{k}|^{2}\right)$
(5.9)

as one can see by performing the change of variables
$\hat{a}_{k}=a_{k}\beta^{-1/4}$. A decrease (increase) in $T$ (in
$\mathcal{P}$) can be read off equivalently as an increase of the spherical
constraint value. Therefore, the occurrence of a localization transition will
be revealed in terms of a critical temperature $T_{c}$ or equivalently a
critical pumping rate $\mathcal{P}_{c}$.

Qualitative information about the presence of a localization transition can be
already traced in the behaviour of the emission spectra, when the temperature
is lowered (or equivalently the pumping rate is increased), see Fig. 1.1. It
can be very clearly seen that, as the pumping is increased the overall
intensity is heterogeneously distributed among the modes. This might hint that
a localization phenomenon in intensity occurs, but it is not enough to
establish it.

In the following sections the analysis is refined by introducing and studying
suitable observables for localization and equipartition breaking transitions.
Data are referred to the simulations of the ML 4-phasor model with PBC, whose
details are reported in Table 4.2.

#####  5.3.1 Participation Ratio

Whether the system truly localizes or not can be ascertained only from the
study of the participation ratio, i.e. the localization order parameter, which
for our system is defined as

$\displaystyle
Y_{2}=\left\langle\frac{\sum_{k=1}^{N}{I}^{2}_{k}}{\left(\sum_{k=1}^{N}{I}_{k}\right)^{2}}\right\rangle=\frac{1}{N^{2}}\left\langle\sum_{k=1}^{N}{I}^{2}_{k}\right\rangle,$
(5.10)

where ${I}_{k}=|a_{k}|^{2}$ and we have used the fact that

$\displaystyle\sum_{k=1}^{N}{I}_{k}=N$ (5.11)

because of equation (2.41), with $\epsilon=1$. In this case, it is irrelevant
to normalize the intensities to the square root of the temperature (as done
for the spectra in Eq. (3.16)), since this factor cancels out between
numerator and denominator. As in the previous chapters, $\langle\cdot\rangle$
denotes the thermal average, i.e. the average computed over all the
uncorrelated equilibrium configurations, which have been sampled.

The dependence on the number of degrees of freedom of $Y_{2}$ can be easily
rationalized in two extreme situations: equipartition and localization of the
mode intensity. Let us consider localization first: in this case a finite
fraction of the whole intensity is taken by a finite number of modes that does
not increase with $N$. That is, in the localized phase, a few modes $k$ have
intensity

$\displaystyle{I}_{k}\propto N,$ (5.12)

whereas all the others have $I_{\neq k}=0$. Then, we have

$\sum_{k=1}^{N}I^{2}_{k}\simeq\sum_{k=1}^{\\#{\rm loc\ modes}}I^{2}_{k}\propto
N^{2}.$

This implies that in a localized phase the participation ratio $Y_{2}$ in the
limit $N\rightarrow\infty$ is a constant that does not depend on $N$:

$\displaystyle\text{localization}~{}\Longleftrightarrow~{}\lim_{N\to\infty}Y_{2}=\mbox{const}.$
(5.13)

On the contrary, in the equipartition phase of nearly homogeneous spectral
intensities, any of the $N$ modes has intensity ${I}_{k}=\mathcal{O}(1)$, so
that in the thermodynamic limit

Figure 5.2: Main: Participation ratio $NY_{2}$ of the mode intensities versus
$T$ at different systems sizes $N$. The vertical line is the asymptotic value
of the critical temperature for the system $T_{c}\simeq 0.61$ (from [Nie+22],
see also Table 4.3). Inset: scaling of the participation ratio near the peak
for the three largest sizes, as $N^{\Psi}Y_{2}$ versus $T$, with
$\Psi=0.35>0$. The peak height scaling is less than $N$, thus Eq. (5.10) tends
to zero as $N\to\infty$.

$\displaystyle\text{equipartition}~{}\Longleftrightarrow~{}Y_{2}\sim\frac{1}{N}.$
(5.14)

Now we are ready to display, in Fig. 5.2, the first important quantitative
information obtained from the study of the equilibrium distribution of the
intensity among modes. In the figure we have plotted for convenience the
average over quenched disorder of $NY_{2}(T)$, which we expect of
$\mathcal{O}(1)$ in the equipartition phase and of $\mathcal{O}(N)$ in a
possible localized phase. We observe that in the high temperature phase
$NY_{2}\sim\mbox{const}$ and, therefore, the system is in the equipartition
regime. Below the critical point, indicated by the vertical line at
$T_{c}=0.61$, which is the glass transition temperature (see [Nie+22]), we
find, instead, an anomaly. Despite the main panel of Fig. 5.2 gives a clear
indication that below the glass temperature $NY_{2}(T)$ grows with the system
size, the collapse of data in the inset shows that the scaling of the growth
is definitely less than $N$. Therefore, the regime is not localized, though it
is not equipartitioned, either. In fact, for $T\lesssim T_{c}$ it occurs to be
$NY_{2}(T)\sim N^{1-\Psi}$ with a value close to $\Psi\simeq 1/3$. This means
that, since

$Y_{2}=\frac{1}{N^{2}}\sum_{k=1}^{N}\langle|a_{k}|^{4}\rangle\sim N^{-\Psi},$