Glassiness and lack of equipartition in random lasers: The common roots of ergodicity breaking in disordered and nonlinear systems

We present here a unifying perspective for the lack of equipartition in nonlinear ordered systems and the low-temperature phase-space fragmentation in disordered systems. We demonstrate that they are just two manifestation of the same underlying phenomenon: ergodicity breaking. Inspired by recent experiments suggesting that lasing in optically active disordered media is related to an ergodicity-breaking transition, we studied numerically a statistical mechanics model for the nonlinearly coupled light modes in a disordered medium under external pumping. Their collective behavior appears to be akin to that displayed around the ergodicity-breaking transition in glasses, as we show measuring the glass order parameter of the replicasymmetry-breaking theory. Most remarkably, we also find that at the same critical point a breakdown of energy equipartition among light modes occurs, the typical signature of ergodicity breaking in nonlinear systems as the celebrated Fermi-Pasta-Ulam model. The crucial ingredient of our system that allows us to find equipartition breakdown together with replica symmetry breaking is that the amplitudes of light modes are locally unbounded, i.e., they are only subject to a global constraint. The physics of random lasers appears thus as a unique test-bed to develop under a unifying perspective the study of ergodicity breaking in statistical disordered systems and nonlinear ordered ones.

The study of ergodicity-breaking transitions is among the most challenging open problems in statistical mechanics. Ergodicity breaking is typical of a wide variety of models, spanning from the physics of viscous liquids to the problem of phase transitions in inference problems. The systematic study of ergodicity breaking in statistical systems with quenched disorder started with the seminal work of Parisi in 1979 [1], where he put forward his proposal of an "infinite number of order parameters" for the low-temperature phase of spin glasses. Since then, the framework to describe the low-temperature phase of spin glasses and glasses became a paradigmatic one: several strong evidences, gathered across the years, showed that the "many-states scenario" [2] is the correct one to understand the properties of a wide class of complex systems. This notwithstanding, the exact breaking of the symmetry between replicas, which is the formal way to represent ergodicity breaking in the presence of quenched disorder [2], is proven only for infinite-dimensional systems. The latter are models where the interaction network between the microscopic degrees of freedom has an infinite-dimensional topology, i.e., any degree of freedom is connected to any other. Therefore, in order to finally establish the relevance of exact results on the glass transition, it worth it to address the study of physical systems truly characterized by an infinitedimensional interaction network among their degrees of freedom. This is indeed the motivation to study the nonlinear interactions of electromagnetic field waves inside an optically active random medium: the complex amplitudes of a normalmode expansion of these waves form an infinite-dimensional interaction network where the scenario depicted by the meanfield theory of spin glasses should be robust. Furthermore, as we are going to show, in these systems ergodicity breaking can be detected in both ways it is historically known to take place: not only as the fragmentation of phase space tracked by nontrivial values of the Parisi's order parameter, but also as the lack of equipartition between the fundamental degrees of freedom. This unprecedented evidence is the most remarkable finding of the present study on random lasers: it is at the same time a consequence of the dilution of nonlinear couplings, which are reduced in number from O(N 4 ) to O(N 3 ) for the effect of a selection rule characteristic of random lasers-but what is relevant is dilution per se-and a consequence of the fact that the complex amplitudes of light modes are not locally bounded, unlike Ising, XY, or Heisenberg spins.
The breaking of ergodicity in the form of lack of equipartition was observed for the first time in the famous Fermi-Pasta-Ulam (FPU) numerical experiment [3], done at Los Alamos Laboratories in 1954 using one of the first computers worldwide, MANIAC I. The FPU "experiment" consisted of the numerical study of the relaxation to equilibrium for a chain of nonlinear oscillators with initial energy put only on a few modes-equipartition was never observed up to the maximum computer time available (at that time). The results were buried at Los Alamos Laboratories until when, in 1965, they were made public by Stanley Ulam. Immediately after that, the famous paper of Zabusky and Kruskal [4] on solitons came out, highlighting the deep connection between the evolution of far-from-equilibrium initial conditions in the FPU model and solitary wave solutions of the Korteweg-de Vries equation. Since then, the investigation of ergodicitybreaking localized solution in nonlinear systems became a very important research topic on its own [3][4][5][6][7][8][9][10], unfortunately without any interplay with the study of the ergodicity breaking "à la Parisi" in disordered systems.
The aim of the present work is to move the first step towards the filling of this gap. We present a system where ergodicity breaking in disordered system phase transitions and the lack of equipartition typical of nonlinear systems are manifest and independently measurable as concomitant phenomena.
In order to fully appreciate the large scope of our results, let us then take a step backward and provide a more specific introduction to random lasers in the context of nonlinear optics. Several diverse complex compounds of optically active random scattering media display lasing under external power pumping above a certain threshold. In their generality, they are called random lasers and are currently investigated in a large number of settings because of their peculiar properties relevant for applications to biosensing, medical diagnostics, on-chip spectroscopy, and optical imaging. They are easy to make, chip and robust, can be very small in size, down to the micron scale, and exhibit unique features such as low degree of spatial coherence, lack of directionality, and biocompatibility [11][12][13][14][15][16][17][18][19][20][21].
Recent experiments on random lasers provided evidence of particularly nontrivial correlations between the so-called shot-to-shot fluctuations of the emission spectrum [22][23][24][25][26][27][28]. The latter are compatible with an organization of mode configurations in clusters of states, similar to that occurring for the multitude of thermodynamic states composing the glassy phase in glass formers [2]. Such a correspondence has been theoretically explained proving the equivalence between the distribution of the intensity fluctuation overlaps (IFO) and the distribution of the overlap between states, the so-called Parisi overlap, the order parameter of the glass transition [29]. The analytic proof has been derived, though, assuming very narrow-band spectra, such that all modes can be considered at the same frequency [30][31][32]. This is not the case, however, for many realistic multimode lasers, both ordered and random. There, the four-wave nonlinear mixing between electromagnetic field modes is controlled by a deterministic selection rule depending on mode frequencies, and interactions are possible only for the quadruplets of modes whose frequencies satisfy the condition with γ being the typical linewidth of a mode. The importance of such selection rules in multimode random laser has been recently experimentally demonstrated in Ref. [33]. Hereafter we call an interaction network built on the mode-locking selection rule in Eq. (1) a mode-locked graph. Study of the repercussions of such a selection rule on nonlinear models is, therefore, a necessary step to understand the fundamental mechanisms at the ground level of the fascinating phenomenon of lasing in random media. As a first step of this work, we reproduce in numerical simulations a narrowing of the emission spectrum across the lasing threshold found in experiments, see Ref. [34] and Fig. 1. We furthermore demonstrate that the narrowing of the spectrum takes place concomitantly with an ergodicitybreaking transition, the glass transition being known to occur in random lasers [22,29,31,32]. We put in evidence two salient features of the glassy phase of light: (1) the breaking of an equipartition between the fundamental degrees of freedom of the system (Fig. 1); (2) the rise of nontrivial glassy correlations among different emission events, landmark of the phase-space shrinking typical of an ergodicity-breaking transition. The discovery that the breaking of equipartition in our model is just another facet of this ergodicity breaking not only provides a strong connection between the physics of nonlinear and of disordered systems but it is important also on its own as, from our analysis, it acquires the status of a new tool to detect glassy phases. We will explain below how breaking of equipartition is revealed by a spectral analysis, which is technically much more simple to perform than the study of the distribution of the state overlap, i.e., the order parameter for many-state disordered systems.
For the sake of clarity we emphasize that the P(q) distribution that we find is not the full replica symmetry-breaking one, describing, e.g., the prototypical spin-glass Sherrington-Kirkpatrick model. There the probability distribution of the overlap P(q) has a nontrivial shape with a continuous q support in the thermodynamic limit N → ∞. Here the P(q) is expected to be bimodal in the thermodynamic limit, though we will see that at finite sizes it appears to be a smooth function of q due to strong finite-size effects.
One last point is important to discuss in the introduction: the lack of equipartition has never been observed so far in disordered glassy systems. The reader might, in fact, be quite suspicious that after 40 years of the discovery of replica symmetry breaking such a fundamental phenomenon was never observed. The crucial point relies on the kind of variables and on the kind of network used in our model. Here we consider a model with continuous and locally unbounded variablesonly one global constraint keeps the energy bounded-and the structure of the interaction network is dense but not complete, therefore inhomogeneous. A specific combination of ingredients that for diverse reasons has never been studied in detail before. Let us briefly explain why. In the statistical physics framework complete graphs with continuous variables, pertaining to the set of p-spin models [35,36], are usually considered as analytically solvable approximations to models with discrete Ising spins. It is the completeness of the graph, i.e., the fact that in the presence of p-body interactions all the independent sets of p variables of degrees of freedom are considered, and hence its homogeneity, which makes an exact solution available [2,35,36]. While always considered for analytic computations, models on complete graphs are vice versa usually not studied in numerical simulations, because of their high computational cost. On the contrary,  Fig. 9 appearing in Ref. [34]. (b) Intensity spectrum I λ = A 2 λ as a function of the wavelength λ averaged over many instances of the quenched randomness, numerical simulations with N = 64 degrees of freedom. The wavelength λ = 2π/k is expressed here in arbitrary units (a.u.). Notice the narrowing of the spectrum at larger values of the pumping rate P. (c) Intensity spectrum I λ as a function of the wavelength λ for a single instance of the quenched randomness, numerical simulations with N = 64 degrees of freedom. Color code of the curves: P increases from bright (yellow) to dark (purple). Notice the crossover from a smooth, almost equipartitioned spectrum at low P to a disordered pattern of isolated peaks at high P. numerical simulations are the main investigation tool for finite-dimensional interaction networks. But in the case of a finite-dimensional topology, continuous and unbounded variables were never studied due to the problem of the power condensation catastrophe [37,38]. Models studied extensively on finite-dimensional topologies are the XY model or Heisenberg spins, where, due to the fixed modulus of local variables, is not possible by construction to observe breaking of equipartition. Hence, also from this respect, our study is completely original.
In conclusion, the dense yet noncomplete interaction network studied in the present work, i.e., the mode-locked graph, has all the ingredients to see some new physics: the structure of the graph is not completely homogeneous and the local variables can change their magnitude, a lucky combination which allowed us to probe the nontrivial local fluctuations which give rise to the so-far-unobserved breaking of equipartition at the glass transition. Most remarkably, the structure of the interaction network studied here is not a build-on-purpose one but is the one produced by the matching condition among mode frequencies in an optical system with nonlinear polarization in a random medium under external energy pumping [39].

I. MODEL AND OBSERVABLES
The disordered model studied here is the mode-locked four-phasor model, characterized by a deterministic selection rule on the interacting quadruplets of modes which emerges naturally from the study of mode dynamics in the stationary regime [30][31][32]40,41]. The dynamic variables of the system are the complex slow amplitudes a k (t ) = A k (t ) e iφ k (t ) of the electromagnetic field expansion in normal modes: Being A k = |a k | ∈ R + and φ k = arg a k ∈ [0, 2π ], the dynamics of the stationary regime can be shown to be a stochastic potential dynamics whose Hamiltonian reads [31,42] An important ingredient of the model is then the implementation of gain saturation, formally rephrased into a global constraint on the intensity [30,41]: The diagonal part of the g ii effective linear coupling is the net gain profile, which plays a role mainly below the lasing threshold. The off-diagonal terms can interpreted as effective interactions between radiant modes [43,44], though only modes of very similar frequency can display a coupling. They are all zero when the frequencies are very distinct, so that the frequency-matching condition of the linear term is never satisfied unless different modes have overlapping frequency. While this is generally true for standard high-quality-factor multimode mode-locked lasers [45], in random media a significant frequency overlap between modes can occur and off-diagonal linear contributions are actually important in the fluorescence regime.
To express the interactions in the slow-amplitude mode basis actually used in the dynamics is actually a complicated problem; the nature and behavior of modes in random media is still the subject of active research [19]. In some cases the solution can be found using some self-consistent procedures, starting from the solution obtained without the nonlinear coupling [46][47][48][49]. In particular, when the nonlinear term is entirely neglected, a possible (though not unique) solution is the one that diagonalizes the linear interaction. When the lasing threshold is overcome, however, the nonlinear term becomes nonperturbatively relevant, and the diagonalization of the linear term does not correspond to a slow amplitude basis anymore, in the most general case of lasing in random media.
Our theory works for any basis, as long as it is well defined and complete, and it focuses on the onset to the lasing phase and on the properties of the modes in this regime. For simplicity, according to this focus we make the working choice of a constant gain profile g ii = g, ∀ i, in the whole wavelength band of the random laser and zero off-diagonal g contributions. Together with the global constraint Eq. (4), this implies the following form for the Hamiltonian: (5) The sum termed FMC is generated by choosing the nonlinear interactions according to a selection rule which depends on mode frequencies, the so-called frequency-matching condition (FMC) [38]: Here γ is the typical linewidth, and we explicitly wrote the permutation that can satisfy the constraint in the ordered sum over The FMC constraint introduces in the topology of the interaction network inhomogeneities such that standard fully connected mean-field approximations used to solve the thermodynamics of disordered systems [2] cannot be applied. Let us consider, for instance, the simple case of a linear dispersion relation with equispaced angular frequencies, ω j = ω 0 + j δ, with δ γ and i = 1, . . . , N. Equation (6) very simply reads as the constraint on the summation indices in Eq. (5), diluting by an order N the number of interactions from the fully connected case, [50]. We call the resulting interaction network the mode-locked graph. The values of the coefficients J ı={i 1 ,i 2 ,i 3 ,i 4 } in Eq. (5) are linked to the electromagnetic normal modes and the nonlinear optical susceptibility as In comparison to the slow-amplitude dynamics these couplings are even slower, so that they can be taken as quenched.
The spatial distribution of the normal modes E(r) and the susceptibility of χ (3) (r) are not exactly known in realistic random lasers; nevertheless, they are expected to be largely inhomogeneous. Here we then consider the interaction coefficients J i 1 ,i 2 ,i 3 ,i 4 to be Gaussian distributed random variables, as the statistical behavior of the system would not depend on the specific form for large N.
Despite the fact that energy is continuously injected and dissipated within a random laser, according to Refs. [30,31,41,42] one can assume an effective equilibrium distribution for the amplitudes: where β is some effective inverse temperature (related to the rate of spontaneous emission of photons), and measures the optical power per mode available to the system. Rescaling (9), the new variables are constrained on the same hypersphere at the cost of a rescaling of the effective temperature as where P is the so-called pumping rate parameter. In these rescaled variables the Boltzmann weight reads making explicit the role of pumping as an effective heat bath for the stationary regimes of the lasing medium. 1

II. NUMERICAL ALGORITHM
We have studied systems with N complex variables a k = |A k |e iφ k interacting with the Hamiltonian in Eq. (5).

A. Exchange Monte Carlo
The sampling of the probability distribution in Eq. (9) was done by means of a parallel tempering Monte Carlo algorithm (PT) [51]. In the PT algorithm one runs K simulations with local Metropolis dynamics for identical replicas of the same system, i.e., with the same quenched disorder J ı . Then exchanges of configurations between heat baths at nearby temperature are proposed, provided that eventually in the thermalized systems for each one of the K copies the configurations are sampled according to the equilibrium Gibbs distribution e −β i E .
We have run K independent simulations at temperatures At each of the 64 steps of the local Metropolis algorithm an exchange of configurations among simulations running at neighboring temperatures is proposed. That is, if the Metropolis algorithm for a i runs with β i and that of a j with β j , one makes the following attempt: {(a i , β i ); (a j , β j )} ⇒ (a j , β i ); (a i , β j ), which is accepted with probability p swap : The replica exchange update is proposed sequentially for all pairs of neighboring temperatures β i and β i+1 . For all simulations the N PT temperatures where taken with a linear spacing in β, i.e., β i+1 = β i + β.

B. The Monte Carlo update under spherical constraint
In the local Metropolis algorithm the configuration of complex "spins" (a 1 , . . . , a N ) is updated with the requirement of 1 We have investigated how the system behaves by varying the pumping rate P. According to Eqs. (9) and (10), in numerical simulations it is identical to fix the constraint and change the effective temperature T = β −1 or work at fixed temperature varying the value of . We have done our simulations varying the temperature T in order to leave a clear term of comparison with the literature on glassy systems (see Methods for the numerical algorithm), but we often discuss our results in terms of pumping rate P. The reader simply needs to always remember that P ∼ 1/ √ T .
keeping k |a k | 2 = const. In order to fulfill such a constraint and the detailed balance, each update is realized choosing at random two spins a i = A i e iφ i and a j = A j e iφ j and proposing an update to randomly chosen values of a i and a j such that This is simply achieved by extracting three random numbers φ i , φ j , and θ with uniform probability in the interval [0, 2π ] and proposing the four simultaneous updates , and A j → A j = C sin(θ ). We have simulated systems at four different sizes: N = 32, 48, 64, 102.
For the parallel tempering we used N PT = 32 temperatures for all sizes.

C. Parallel computation on GPU
The update of the variables A i must be done sequentially, since the interaction network is dense and there is no way to partition the variables in subsets which can be updated independently in parallel. We have, instead, implemented parallelization on GPU graphic cards in two ways: (i) the update energy shift and (ii) the PT replicas dynamics between configuration exchanges.
In order to accept or reject the update of two spins a i and a j , one has to compute the energy update on N (i, j) The calculation of each E κ is realized in parallel on a distinct kernel on GPU. Moreover, the execution of the K evolutions at different temperatures of the parallel tempering is quite naturally implemented in parallel on the GPU. We have run simulations on two type of graphic cards: GTX 680 and Tesla K20, achieving an overall speedup with respect to the sequential code on CPU, up to a factor of 8.

III. BREAKING OF EQUIPARTITION
Let us first compare the dependence on the pumping rate P of the emission spectrum measured in experiments [panel (a) of Fig. 1, data taken from [52]] with the emission spectrum in our numerical simulations averaged over the random coupling distribution [panel (b)]; the increase of P is accompanied by the same kind of narrowing. We stress that the symmetry of the spectral distribution around the central wavelength is due to the fact that the gain is taken as constant in this case, cf. Eq. (5), whereas in reality a material-dependent, nonconstant gain curve g(λ) is present.
The equipartition breaking obtained by raising P is then even clearer if we consider the behavior for a single instance of disorder rather than its average over many realizations, see panel (c) of Fig. 1. From a smooth profile at low P a pattern of random sharp peaks emerges in the spectrum at high P. Let us point out that the spectrum shown in panel (c) of Fig. 1 represents an average over a large time window (see Appendix), in agreement with the experimental nature of the emission spectrum, which is a time-integrated signal. From the literature on nonlinear systems [5,7], we know that the degree of equipartition in the spectrum can be measured by the so-called spectral entropy: whereÎ k is the normalized thermal average of the intensity for given wave number k. From the spectral entropy one then defines the effective number of degrees of freedom [5,7]: where n eff = 1 for perfect equipartition and n eff = 1/N when the total energy is concentrated in a single mode [5,7]. For our system the behavior of its average, n eff , as a function of the pumping rate P is reported in Fig. 2(a) (data are plotted as a function of P −2 = T ). At low pumping rate we find n eff ≈ 1, which signals a good degree of equipartition for all the sizes N studied. On the contrary, by increasing P we find that n eff rapidly decreases and the decrease depends on N. We identify as the pumping value that at which n eff starts decreasing below 1 with a size-dependent lasing threshold P c (N ) = √ T c (N ). The figure clearly shows that the larger the N, the steeper the n eff decrease.
A stronger indication that we are dealing with a first-order transition with respect to n eff comes from the study of the Binder B and the bimodality b, parameters of the distribution P(n eff ), displayed in Figs. 2(b) and 2(c). The Binder parameter, which is simply defined as the ratio between the fourth and second moment of the distribution P(n eff ), with n eff = n eff − n eff as a measure of the deviation from Gaussianity. The order parameter of a first-order transition is typically characterized by a Gaussian distribution, which is centered around two different values, depending on whether its is probed well below or well above the transition temperature (if the parameter control is temperature). In the coexistence region where the two phases have comparable free energies, the distribution of the order parameter is bimodal. For this reason the Binder parameter as a function of the temperature also has a characteristic behavior, pointed out first by Binder [53]; it is a nonmonotonic reversed bell behavior with a maximal deviation from Gaussianity in the coexistence region, where the distribution is bimodal. Indeed, in order to be completely sure that the deviation from Gaussianity is due to bimodality and not to other reasons, we also measured the bimodality parameter b, defined as where n is the number of data composing the histogram of the probability distribution, κ is the curtosis of the distribution, i.e., the same ratio of fourth-and second-order momentum appearing in the definition of the Binder parameter-see Eq. (17)-and γ is the skewness, defined as Indeed, what we find is precisely what Binder presents as the phenomenology of a first-order transition: the value P c (N ) at which B is mostly far from a Gaussian corresponds to the peak of the bimodality indicator b, as well as to the breaking point of n eff . Eventually, in Fig. 2(d) we compare the spectral profiles obtained for pumping rates either well below or well above the lasing threshold. The evidence that n eff behaves as the order parameter is very important for the purpose of the whole discussion. The main point of presenting data here on equipartition breaking is to compare them with the standard way to detect ergodicity breaking in disordered systems, i.e., studying the Parisi's order parameter. As we are going to show, the glass transition that we detect for our system is a first-order one with respect to its order parameter. It is therefore very important that also the "order parameter" for equipartition breaking behaves as the parameter of a first-order transition; this puts on a solid basis the whole scenario we are proposing, namely, that we are looking at the same underlying phenomenon from different perspectives. The goal of next section is to show that the distribution of the equilibrium overlap between replicas confirms this picture.

IV. GLASS TRANSITION
The so-called random first-order transition (RFOT) [54][55][56] is a mixed-order ergodicity-breaking transition characterized by a specific heat anomaly and a discontinuity in the order parameter. The latter is identified as the probability distribution of the overlap q between equilibrium configurations at a given temperature [2]. In structural glasses, across the glass transition the distribution P(q) passes from a unimodal distribution at high temperatures, where there is a unique state, to a bimodal one with a secondary peak (actually two symmetric peaks) at low temperatures because of states fragmentation into nonequivalent clusters [2]. This is the typical behavior when phase space splits in disjoint ergodic components: configurations inside the same ergodic component have typical overlap q 1 , the overlap between configurations in disjoint ergodic components being q 0 < q 1 .
Let us stress that this kind of ergodicity breaking is not the so-called "spin-glass" transition of the Sherrington-Kirkpatrick model [57], where the probability distribution of the overlap P(q) is expected to take a nontrivial shape, different from a bimodal one, even for N → ∞. In the RFOT case, due to finite-size effects, P(q) appears continuous but it tends to a bimodal distribution in the thermodynamic limit, as is expected for the glass transition in the p-spin model with p > 2. The difference is the phase transition is well accounted for by the specific heat behavior. C v displays a peak, divergent with N, in our simulations, whereas in the properly defined spin glass no specific heat singularity is expected to occur at the spin-glass transition.
We show now that the four-phasor random laser model displays these RFOT features at the critical pumping power.
Let us start from the specific heat anomaly. In panels (a) and (b) of Fig. 3 the specific heat C V curves for different sizes of the system are shown, where C V per phasor is The and represent, respectively, thermal average and average over quenched disorder. In order to have a term of comparison, in Fig. 3 we also show the results for the model (called here random diluted model) that has the same number of modes, the same distribution of coupling values, and with the same total number of interactions, yet whose set of interacting modes are not chosen according to Eq. (7) but, rather, uniformly randomly selected from the set of all quadruplets (see Methods for details). At all sizes N, C V /N has a characteristic nonmonotonic behavior, with the position of the peak depending on N. We identify this point as a finite-size "critical temperature" T c (N ) = 1/P 2 c . The good collapse of the curves at different N, shown in the inset of Figs. 3(a) and 3(b), demonstrates critical scaling, a typical feature of second-order phase transitions. The width of the C V scaling region for the model with random dilution is τ ∼ N −1/2 , where τ = T /T c − 1: the scaling exponent 1/ν = 1/2 is the same predicted by the simplest mean-field theory for a second-order phase transition and is also in agreement with the reference mean-field model of the RFOT [54]. On the contrary, the scaling region for our four-phasor mode-locked model follows the behavior τ ∼ N −3/2 , with an exponent 1/ν = 3/2, not compatible with a mean-field theory of second-order transition. This unexpected exponent 3/2 is probably due to the correlations in the interaction network induced by the frequency-matching condition, Eq. (7). Though it might be a preasymptotic effect due to too small finite sizes of the simulated systems, with the present data we cannot rule out the possibility that it is a behavior persistent even in the thermodynamic limit.
In order to detect the breaking of ergodicity, we then studied the overlap in the numerical simulations, where it can be measured as the similarity between two mode configurations evolving at equilibrium at the same temperature and with the same realization of random couplings: two replicas [2]. Labeling with the Greek indices α and γ two replicas, we have initially measured (cf. Methods) the phasor overlap q αγ : The Greek indices in Eq. (20) denote different replicas, i.e., independent configurations at equilibrium at the same temperature T . The Parisi's overlap is characterized by a low-temperature nontrivial distribution in the presence of a glass phase [2]. In panels (c) and (d) of Fig. 3 the phasor overlap distribution P(q) is shown, respectively, for the randomly diluted fourphasor model [panel (c)] and the mode-locked model [panel (d)] at different values of P. In both cases, in the low-P ergodic phase, it turns out to be a symmetric Gaussian. Then, for P > P c , in the case of random bond dilution one finds secondary peaks at a finite distance from the origin, signaling a glassy broken-ergodicity phase. In the mode-locked graph, actually, this effect is not very pronounced, and shoulders are displayed at the simulated sizes and powers, rather than proper side peaks. As for the specific heat, finite-size effects turn out to be stronger in the mode-locking graph than in the randomly bond-diluted one.
We further considered another observable-the quadruplets overlap probability distribution P (Q) (see Methods). The quadruplet overlap Q αβ between the energy stored in the same set of four coupled modes, at equilibrium in state α and in state γ , reads where κ runs over the ordered list { ı} of all N 4 nonzero four-body couplings. While the phasor overlap is computed as the average over N correlated random variables, i.e., the local phasor overlapsā α k a γ k , the quadruplet overlap Q αγ μ is the average of N 4 = O(N 3 ) variables. We hence expect it to be less plagued by finite-size effects. This is, indeed, the case, as shown in Fig. 4. While at low pumping rate P the quadruplets overlap has a very peaked distribution at Q 0, for P ≈ P c we find a clear signature of a secondary peak at Q > 0. In the inset of Fig. 4(a) we see the corresponding multimodality parameter b.
The definition of b is the same given in Eq. (18) above, where skewness γ and curtosis κ are now defined as follows: with q = q − q . We find that the region where the parameter b signals a bimodal distribution of the overlap is precisely the interval of pumping rates around P c . For the study of the plaquette overlap Q distribution, we could 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. We have presented the numerical evidence that in a statistical mechanics model for random lasers there is a RFOT glass transition concomitant with the breaking of equipartition of energy among the modes. One question now is to which extent this picture can be assessed even in experiments. The available technology for the measurements of light mode phases [58,59] applies only to high-power directional impulses: unfortunately, this is not the operating regime of standard random lasers [17,60]. Therefore we cannot rely on observables which require the measurement of phases. The phasor and quadruplet overlaps defined in Eqs. (20) and (22) depend on phases, so they are not quantities that can be measured in current experiments. A further overlap C αβ between intensity fluctuations k ∝ I k − I k can then be introduced: where I k = A 2 k . At the mean-field narrow-band level, the intensity fluctuation overlap (IFO) C αβ is proved to be in a one-to-one correspondence with the standard overlap, i.e., C αγ ∝ q 2 αγ ∀ α, γ [29]. The advantage of C αβ with respect to q αβ is that it can be measured in real random lasers [22][23][24][25][26][27][28]. In panel (b) of Fig. 4, the distribution P (C) is shown for N = 64 at four different values of the pumping rate P, two above and two below the critical P c , as determined from the caloric curve (Fig. 2). At first sight there is no clear evidence of secondary peaks at high pumping rates for this system size, although non-Gaussian tails appear in the vicinity of the transition. It is the study of the Binder parameter B dependence on P which reveals how P (C) brings the signature of a first-order transition: data are shown in the inset of panel (b) of Fig. 4. The behavior of B as a function of P (plotted as B vs P −2 = T to have a clearer term of comparison with the literature) is the one characteristic of first-order transitions [53]. This behavior of the IFO order parameter is, therefore, also consistent with the RFOT scenario and with the study of the breaking of equipartition in the spectrum: in all the three cases the ergodicity-breaking parameter behaves as the order parameter of a first-order transition.

V. DISCUSSION
By means of Monte Carlo numerical simulations we have analyzed a statistical mechanical model of the nonlinear interactions of light in a random medium under external optical pumping. We have shown that the model reproduces the narrowing of the emission spectrum characterizing the onset of lasing in experiments, see Fig. 1, and how it signals an ergodicity-breaking transition. Indeed, the onset of the random lasing regime from an incoherent regime of fluorescence displays the properties of a glassy phase transition, characterized by a diverging specific heat (typical of continuous transitions) and a discontinuous order parameter, the overlap, a feature of first-order transitions. This is termed random first-order transition in the literature [55].
Further on, thanks to the quite general properties of the model analyzed, we were able to show unprecedented evidence of a deep connection between lack of equipartitiontypical of ergodicity breaking in ordered systems with nonlinear interactions [3,5,9]-and the breaking of ergodicity as described within the paradigm of replica symmetry breaking in disordered systems [2].
The mode-locked four-phasor model is the first example of a system where ergodicity breaking manifests itself at the same time as a breaking of the symmetries between replicas and as a lack of equipartition. We stress that the possibility to reveal this phenomenon is due to the following two salient features of the model: the O(N ) dilution of the nonlinear couplings with respect to the corresponding meanfield model of the p-spin class, induced by the frequencymatching condition, and the local unboundedness of the light mode amplitudes. While the rule for the dilution of couplings is specific of random lasers, it is dilution itself which allowed us to see equipartition breaking. The model of the random laser presented here is thus, in our opinion, only a specific instance of a broader class of systems characterized by the concomitance of equipartition and replica symmetry breaking: we expect a similar phenomenon to take place for all p-spin models on inhomogeneous networks and with locally unbounded variables. This important result suggests a possible way to overcome the intrinsic difficulty usually encountered in the measure of q, which is not a singleexperiment observable. The standard protocol is that one needs to compare the results of several experiments done on the same sample or of several numerical simulations with the same realization of quenched disorder to obtain a measure P(q). Using the jargon of disordered systems, one needs more replicas of the same system. As we have shown, the occurrence of a nontrivial P(q) is simultaneous with the loss of spectral equipartition, and one can thus simply take advantage of the latter to detect the ergodicity-breaking glass transition.
In conclusion, as weird as it may sound, the results presented in this work reinforce the notion that light in random media really looks like a glassy system and offers, among all physical systems, a good benchmark to test the existence and the hidden nature of a glass transition and its connection to nonlinearity. The first step of the numerical study is the generation of the mode-locked graph with disordered couplings. Our goal is to study what happens beyond the narrow-band mean-field approximation of [29,31,42,61], in particular, when the nonlinear interaction is that of Eq. (5). Operatively, it is easier to describe the structure of the mode-locked interaction network as a bipartite graph where interaction nodes J μ , labeled by Greek letters, are connected to variable nodes A k , labeled with Latin letters. Since we have a four-body interaction, each interaction node is always attached to four variable nodes and is defined by the ordered list of their indices, J μ (i, j, k, l ). The order is relevant, because from the point of view of the energy stored in the interaction [see Eq. (5)] there are nonequivalent permutations. The steps to generate the mode-locked graph are as follows: (1) A virtual complete graph with N 4 interaction nodes is generated.
(2) For each interaction node the three nonequivalent index permutations in Eq. (5) should be considered: P μ (i 1 , i 2 , i 3 , i 4 ), P μ (i 2 , i 1 , i 3 , i 4 ), and P μ (i 1 , i 2 , i 4 , i 3 ), even though for the index order i 1 < i 2 < i 3 < i 4 only P μ (i 2 , i 1 , i 3 , i 4 ) can be satisfied at most. Each time that the mode frequencies (or indices) satisfy condition (6), the corresponding interaction of the virtual graph is added to the real graph.
(3) The procedure at point 2 is repeated until a preassigned number of interactions in the complete graph is reached. For computational reasons this number is the largest power of 2 below the total number of possible couplings satisfying (6).
For large N, the above procedure tends to cut O(N ) out of all interacting quadruplets [50]. Operatively, for a system with N complex variables we have drawn a bipartite graph with a number of interactions scaling as O(N 3 ) and equal to the power of 2, soon smaller than the number of all possible interactions fulfilling the FMC constraint. The number of interaction nodes N 4 corresponding to each N is listed in Table I.
Concerning the structure of the topology of the interaction network, it is important to stress that the dimensionality of the disordered optical medium in real space is scarcely important for the thermodynamics of the problem: the interaction among modes remains in any case highly nonlocal in the basis of normal modes. The only quantity which depends on the realspace dimensionality is the spatial overlap between the normal modes, the information about which is stored in the disordered coefficients J ı of the Hamiltonian in Eq. (5) as [31,61] where χ (3) ν is the nonlinear susceptibility of the system. We do not consider contributions from the first nonlinear term χ (2) in the polarization expansion because of the relatively limited bandwidth of the emission spectra in random lasers and the consequent lack of second harmonic generation, but that term can be easily inserted, leading to no qualitative difference in the results. Furthermore, we do not explicitly consider linear polarization contributions, which would read g k = g(ω k ) being the net gain profile (gain taken away loses due to the open cavity) that can be expressed as a function of χ (1) (ω, r). In our model, though, where gain saturation is expressed by the global constraint on the overall intensity distributed in the system, and considering for simplicity a constant net gain profile g k g, this contribution amounts to a constant H 2 = −g, with no role in the dynamics. The random laser couplings as expressed in Eq. (A1) are, in general, disordered because modes display different spatial shape and extension [18,62]. The constituents of the integrals in Eq. (A1) are very difficult to calculate from first principles. The only specific form of the nonlinear susceptibility has been computed by Lamb [39,63] for few-modes ordered lasers, and no analog study for random lasers has been performed so far, to our knowledge. Integrals like Eq. (A1) in a random medium can be regarded as a sum over many random variables. Different couplings involving a given mode might, in general, be correlated [64]. Since, however, we are interested in the critical behavior, thus in the large size limit of our simulated systems and since correlations decay with the size of the system, we adopt as working hypothesis a Gaussian distribution for each J ı : The scaling of variance with N, J 2 ∼ N −2 , guarantees energy extensivity. There are very many couplings, but each one is vanishingly small. Macroscopic phenomena such as lasing occur because the system undergoes a transition to a collective behavior. Even in the presence of randomness, cooperativity is the leading mechanism. We further stress that, from the perspective of probing RFOT, considering correlated J's leads to a qualitatively analog phase diagram, as it is well known in spin-glass systems such us, e.g., the random orthogonal model [65,66].
If we look at interactions in the space of normal modes the system is infinite dimensional: any degree of freedom participates to an infinite [O(N 2 )] number of interactions in the thermodynamic limit N → ∞. That is why we expect the mean-field glass transition scenario drawn in [29,31] to be quite robust for the mode-locked p phasor, even if the narrow-band hypothesis [30] is removed. Last but not least, the nonlocality of interactions between light modes also guarantees that phenomena like energy localization, a pathology of sparse networks [37,38], are avoided.

Replicas in numerical simulations
Replicas are independent equilibrium configurations sampled with the same quenched disorder. This definition corresponds to the protocol used in numerical simulation. One "replica" of the system is represented by the swarm of N PT configurations used for a given instance of the parallel tempering Monte Carlo dynamics. Different instances of the PT dynamics characterized by the same set of quenched couplings J ı and the same interaction network between modes are different replicas.
For N = 32, N = 48, and N = 64, for each disorder instance we simulated four replicas, which gave us the availability of six independent values of q αβ : q 12 , q 13 , q 14 , q 23 , q 24 , q 34 . For N = 102 we simulated two replicas for each instance of the disorder. To accumulate statistics for P(q), we measured values of q αβ comparing replicas at the same iteration of the PT dynamics, each of 640 iterations. Since the distribution P(q) is not self-averaging [2], for each size of the system we have sampled the equilibrium measure for N sample ≈ 100 instances of the disorder.
It is useful to clarify also how we measured in practice all thermal averages indicated with angular brackets, O [A] . Once that the system reaches equilibrium at a given temperature/power, at a time in Monte Carlo sweep units of N term < N sweep /4, thermal averages were measured as time averages along the dynamics, along the second half of each run: (A5)