Glassiness and Lack of Equipartition in Random Lasers: the common roots of Ergodicity Breaking in Disordered and Non-linear Systems

TL;DR

This paper unifies the understanding of ergodicity breaking in disordered and non-linear systems, showing that random laser behavior and glassy phase-space fragmentation share a common root, with equipartition breakdown occurring at a critical transition.

Contribution

It introduces a unifying perspective linking ergodicity breaking in disordered and non-linear systems through a statistical mechanics model of random lasers, highlighting equipartition breakdown at the transition.

Findings

Ergodicity breaking occurs in a model of random lasers at a critical point.

Equipartition among light modes breaks down at the same transition.

The system exhibits replica symmetry breaking similar to glassy systems.

Abstract

We present here for the first time a unifying perspective for the lack of equipartition in non-linear 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 the one displayed around the ergodicity-breaking transition in glasses, as we show measuring the glass order parameter of the replica-symmetry-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 non-linear systems as the celebrated Fermi-Pasta-Ulam model. The crucial ingredient of our system which 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 non-linear ordered ones.

Figures

Tables

Table 1. Table 1: Details for the simulations for different system sizes N 𝑁 N . The number of quadruplets involved in the Hamiltonian, N 4 subscript 𝑁 4 N_{4} , is always a power of 2 2 2 .

$N$	$N_4$	$T_{\text{min}}$	$T_{\text{max}}$	$N_{\text{PT}}$	$N_{\text{sweep}}$	$N_{\text{sample}}$
32	$2^{11}$	0.2	2.0	32	$2^{20}$	100
48	$2^{13}$	0.3	2.0	32	$2^{20}$	100
64	$2^{14}$	0.5	1.4	32	$2^{20}$	100
102	$2^{16}$	0.7	1.2	32	$2^{21}$	100

Equations

$$|{\omega }_{i_1}-{\omega }_{i_2}+{\omega }_{i_3}-{\omega }_{i_4}|<\gamma ,$$

$$\mathbf{E}(\mathbf{r},t)=\sum _{k=1}^N{a}_k(t)e^{i{\omega }_{k}t}{\mathbf{E}}_k(\mathbf{r})+\text{c.c.}.$$

$$\mathcal{H}[\mathbf{a}]=-\sum _{i_1<i_2}^{\mathrm{FMC}}{g}_{i_1{i}_2}{\bar{a}}_{i_1}{a}_{i_2}$$

$$+\sum _{i_1<i_2<i_3<i_4}^{\mathrm{FMC}}{J}_{\vec{\text{i}}}{\bar{a}}_{i_1}{a}_{i_2}{a}_{i_3}{\bar{a}}_{i_4}+\text{mbox}c.c.$$

$$=-\sum _{i_1<i_2}^{\mathrm{FMC}}{g}_{i_1{i}_2}{A}_{i_1}{A}_{i_2}\cos ({\varphi }_{i_1}-{\varphi }_{i_2})$$

$$-\sum _{i_1<i_2<i_3<i_4}^{\mathrm{FMC}}{J}_{\vec{\text{i}}}{A}_{i_1}{A}_{i_2}{A}_{i_3}{A}_{i_4}\cos ({\varphi }_{i_1}-{\varphi }_{i_2}+{\varphi }_{i_3}-{\varphi }_{i_4}),$$

$$\vec{\text{i}}\equiv i_1,i_2,i_3,i_4$$

$$ϵN=\sum _{k=1}^N{A}_k^2$$

$$\mathcal{H}[\mathbf{a}]=-gϵN$$

$$-\sum _{i_1<i_2<i_3<i_4}^{\mathrm{FMC}}{J}_{\vec{\text{i}}}{A}_{i_1}{A}_{i_2}{A}_{i_3}{A}_{i_4}\cos ({\varphi }_{i_1}-{\varphi }_{i_2}+{\varphi }_{i_3}-{\varphi }_{i_4}),$$

$$|{\omega }_{i_2}-{\omega }_{i_1}+{\omega }_{i_3}-{\omega }_{i_4}|\lesssim \gamma .$$

$$|i_1-i_2+i_3-i_4|=0$$

$$J_{i_1,i_2,i_3,i_4}\propto \int d\mathbf{r}\sum _{\vec{\nu }}^{x,y,z}{E}_{i_1}^{({\nu }_1)}(\mathbf{r})E_{i_3}^{({\nu }_2)}(\mathbf{r})E_{i_3}^{({\nu }_3)}(\mathbf{r})E_{i_4}^{({\nu }_4)}(\mathbf{r}){\chi }_{(\vec{\nu })}^{(3)}(\mathbf{r}|{\omega }_{i_1},{\omega }_{i_2},{\omega }_{i_3},{\omega }_{i_4})$$

$$P(a_1,\dots ,a_k)=e^{-\beta \mathcal{H}[\mathbf{a}]}\delta \left(ϵN-\sum _{k=1}^N{A}_k^2\right),$$

$$\beta \to \beta {ϵ}^2={\mathcal{P}}^2$$

$$\rho [\{\mathbf{a}\}]=\frac{1}{Z}\exp \{-{\mathcal{P}}^2\mathcal{H}[\mathbf{a}]\},$$

$$p_{\text{swap}}=\text{min}\left[1,e^{-({\beta }_i-{\beta }_j)(\mathcal{H}[{\mathbf{a}}_j]-{\beta }_j\mathcal{H}[{\mathbf{a}}_i])}\right]$$

$$|a_i{|}^2+|a_j{|}^2=|a_i^{\prime }{|}^2+|a_j^{\prime }{|}^2=C^2.$$

$$\Delta E=\sum _{\kappa =1}^{N_4^{(i,j)}}\Delta E_{\kappa }$$

$$S_{\text{sp}}$$

$${\hat{\mathcal{I}}}_k$$

$$n_{\text{eff}}=\frac{\exp \left(S_{\text{sp}}\right)}{N},$$

$$\mathcal{B}=\frac{1}{2}\left(3-\frac{\overline{(\Delta n_{\text{eff}}{)}^4}}{{\left[\overline{(\Delta n_{\text{eff}}{)}^2}\right]}^2}\right),$$

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

$$\gamma =\frac{\overline{(\Delta n_{\text{eff}}{)}^3}}{{\left[\overline{(\Delta n_{\text{eff}}{)}^2}\right]}^{3/2}}.$$

$$\frac{C_V}{N}=\frac{\overline{⟨E^2⟩-⟨E{⟩}^2}}{N{T}^2}$$

$$q^{\alpha \gamma }=\frac{1}{N}\sum _{k=1}^N{\bar{a}}_k^{\alpha }{a}_k^{\gamma }=\frac{1}{N}\sum _{k=1}^N{A}_k^{\alpha }{A}_k^{\gamma }\cos ({\varphi }_k^{\alpha }-{\varphi }_k^{\gamma }).$$

$${\mathcal{E}}_{\vec{i}}$$

$${\mathcal{Q}}^{\alpha \gamma }=\frac{1}{N_4}\sum _{\kappa =1}^{N_4}{\mathcal{E}}_{\kappa }^{\alpha }{\mathcal{E}}_{\kappa }^{\gamma },,$$

$$\kappa$$

$$\gamma$$

$${\mathcal{C}}^{\alpha \beta }$$

$${\Delta }_k^{\alpha }$$

$$J_{i_1{i}_2{i}_3{i}_4}=\frac{ı}{2}{\omega }_{i_1}{\omega }_{i_2}{\omega }_{i_3}{\omega }_{i_4}{\int }_{V}d\mathbf{r}\sum _{\vec{\nu }}^{x,y,z}{\chi }_{\vec{\nu }}^{(3)}({\omega }_{i_1},{\omega }_{i_2},{\omega }_{i_3},{\omega }_{i_4};\mathbf{r})E_{i_1}^{{\nu }_1}(\mathbf{r})E_{i_2}^{{\nu }_2}(\mathbf{r})E_{i_3}^{{\nu }_3}(\mathbf{r})E_{i_4}^{{\nu }_4}(\mathbf{r}).$$

Full text

The study of ergodicity-breaking transitions is one 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 Parisi79 , 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” Mezard87 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 presence of quenched disorder Mezard87 , 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.

It is therefore worth, in order to finally establish the relevance of exact results on the glass transition, to address the study of physical systems truly characterized by an infinite-dimensional interaction network among their degrees of freedom. This is indeed the motivation to study the non-linear interactions of electromagnetic field waves inside an optically-active random medium: the complex amplitudes of a normal mode expansion of these waves form an infinite-dimensional interaction network where the scenario depicted by the mean-field theory of spin-glasses should be robust. Furthermore, as we are going to show, in these systems ergodicity breaking can be detected in both the ways it is historically known to take place: not only as the fragmentation of phase space tracked by non-trivial 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 non-linear couplings, which are reduced in number from $\mathcal{O}(N^{4})$ to $\mathcal{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 Fermi55 , done in Los Alamos laboratories in 1954 using one of the first computers worldwide, MANIAC I. The FPU ’experiment’ consisted in the numerical study of the relaxation to equilibrium for a chain of non-linear oscillators with initial energy put only on few modes: equipartition was never observed up to the maximum computer time available (at that time). The results were buried in Los Alamos laboratories until when, in 1965, they were made public by Stanley Ulam. Immediately after that, the famous paper of Zabusky and Kruskal Zabusky65 on solitons came out, highlighting the deep connection between the evolution of far from equilibrium initial conditions in the FPU model and solitary waves solutions of the Kortweg De-Vries equation. Since then, the investigation of ergodicity-breaking localized solution in non-linear systems became a very important research topic on its own Fermi55 ; Zabusky65 ; Livi85 ; Livi87 ; Cretegny98 ; KLR04 ; Chaos05 ; Benettin13 , 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, for the first time, ergodicity breaking in disordered systems phase transitions and the lack of equipartition typical of non-linear 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 non-linear 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 bio-sensing, 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 bio-compatibility Cao03 ; Florescu04b ; Anni04 ; Cao05 ; LCOW07 ; vanderMolen07a ; Wiersma08 ; Fallert09 ; Andreasen11 ; Leonetti13c ; Bachelard14 .

Recent experiments on random lasers provided evidence of particularly non-trivial correlations between the so-called shot-to-shot fluctuations of the emission spectrum Ghofraniha15 ; Gomes16 ; Pincheira16 ; Supratim16 ; Tommasi16 ; Lopez18 ; Tommasi18 . The latter are compatible with an organization of mode configurations in clusters of states, similar to the one occurring for the multitude of thermodynamic states composing the glassy phase in glass formers Mezard87 . 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 Antenucci15f . The analytic proof has been, though, derived assuming very narrow-band spectra, such that all modes can be considered at the same frequency Gordon02 ; Antenucci15a ; Antenucci15b . This is not the case, however, for many realistic multimode lasers, both ordered and random. There, the four-waves non-linear mixing between electromagnetic field modes are controlled by a deterministic selection rule depending on modes frequencies: interactions are possible only for the quadruplets of modes whose frequencies satisfy the condition

[TABLE]

with $\gamma$ being the typical line-width of a mode. The importance of such selection rule in multimode random laser has been recently experimentally demonstrated in Antenucci19 . We will call, hereafter, an interaction network built on the mode-locking selection rule in Eq. (1) a mode-locked graph. The study of the repercussions of such selection rule on nonlinear models is, therefore, a necessary step to understand the fundamental mechanisms at the ground of the fascinating phenomenon of lasing in random media.

As a first step of this work, we reproduce in numerical simulations the narrowing of the emission spectrum across the lasing threshold found in experiments, see  Cao03r and Fig. 1. We, further, demonstrate that the narrowing of the spectrum takes place concomitantly with an ergodicity-breaking transition, the glass transition, known to occur in random lasers Antenucci15a ; Antenucci15b ; Antenucci15f ; Ghofraniha15 . We put in evidence two salient features of the glassy phase of light: 1) the breaking of equipartition between the fundamental degrees of freedom of the system [Fig.1]; 2) the rise of non-trivial 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 does not only provide a strong connection between the physics of non-linear 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 easy 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 underline 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\to\infty$. 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 variables—only 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 shortly explain why. In the statistical physics framework complete graphs with continuous variables, pertaining to the set of $p$-spin models Kirkpatrick87b ; Crisanti95 , 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 presence of $p$-body interactions all the independent $p$-upletes of degrees of freedom are considered, and hence its homogeneity, that makes an exact solution available Mezard87 ; Kirkpatrick87b ; Crisanti95 . 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, 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 Antenucci15c ; Antenucci15d . 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 non-complete 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 non-trivial 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 SargentIII78 .

I Model and Observables

The disordered model studied here is the mode-locked $4$-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 Gordon02 ; Antenucci15a ; Antenucci15b ; Antenucci15e ; Antenucci16 . The dynamic variables of the system are the complex *slow * amplitudes $a_{k}(t)=A_{k}(t)~{}e^{i\phi_{k}(t)}$ of the electromagnetic field expansion in normal modes

[TABLE]

Being $A_{k}=|a_{k}|\in\mathbb{R}^{+}$ and $\phi_{k}=\arg a_{k}\in[0,2\pi]$, the dynamics of the stationary regime can be shown to be a stochastic potential dynamics whose Hamiltonian reads  Angelani06a ; Antenucci15a :

[TABLE]

An important ingredient of the model is then the implementation of gain saturation, formally rephrased into a global constraint on the intensity Gordon02 ; Antenucci16 :

[TABLE]

The diagonal part of the $g_{ii}$ effective linear coupling is the net gain profile, that plays a role mainly below the lasing threshold. The off-diagonal terms can interpreted as effective interactions between radiant modes Hackenbroich03 ; Viviescas03 , though only modes of very similar frequency can display a coupling. They are all zero when the frequencies are well 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 multi-mode mode-locked lasers Haus00 , 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. Actually the nature and behavior of modes in random media is still the subject of active research Andreasen11 . In some cases the solution can be found using some self-consistent procedures, starting from the solution obtained without the non-linear coupling Tureci06 ; Tureci08 ; Tureci09 ; Rotter14 . In particular, when the non-linear 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 non-linear term becomes non-perturbatively 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 far as it is well defined and complete, and it focuses on the onset lo 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$, $\forall\,i$, in the whole wave-length 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

[TABLE]

The sum termed FMC is generated by choosing the non-linear interactions according to a selection rule which depends on mode frequencies, the so-called Frequency Matching Condition (FMC) Antenucci15d :

[TABLE]

Here $\gamma$ is the typical line-width and we explicitly wrote the permutation that can satisfy the constraint in the ordered sum over $i_{1}<i_{2}<i_{3}<i_{4}$. 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 Mezard87 cannot be applied. Let us consider for instance the simple case of a linear dispersion relation with equispaced angular frequencies, $\omega_{j}=\omega_{0}+j~{}\delta$, with $\delta\ll\gamma$ and $i=1,\ldots,N$. Eq. (6) very simply reads as the constraint

[TABLE]

on the summation indices in Eq. (5), diluting by an order $N$ the number of interactions from the fully connected case, $N(N-1)(N-2)(N-3)/24$ Marruzzo18 . We call the resulting interaction network the mode-locked graph.

The values of the coefficients $J_{\vec{\i}=\{i_{1},i_{2},i_{3},i_{4}\}}$ in Eq. (5) are linked to the electro-magnetic normal modes and the nonlinear optical susceptibility as

[TABLE]

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 $\bm{E}(\bm{r})$ and the susceptibility of $\chi^{(3)}(\bm{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 Gordon02 ; Angelani06a ; Antenucci15a ; Antenucci16 one can assume an effective equilibrium distribution for the amplitudes:

[TABLE]

where $\beta$ is some effective inverse temperature (related to the rate of spontaneous emission of photons) and $\epsilon$ measures the optical power per mode available to the system. Rescaling $A_{k}\to A_{k}/\sqrt{\epsilon}$ in Eq. (9), the new variables are constrained on the same hypersphere at the cost of a rescaling of the effective temperature as

[TABLE]

where $\mathcal{P}$ is the so-called pumping rate parameter. In these rescaled variables the Boltzmann weight reads

[TABLE]

making explicit the role of the pumping as effective heat bath for the stationary regimes of the lasing medium 111We have investigated how the system behaves varying the pumping rate $\mathcal{P}$. According to Eq. (9) and Eq. (10), in numerical simulations it is identical to fix the constraint $\epsilon$ and change the effective temperature $T=\beta^{-1}$ or work at fixed temperature varying the value of $\epsilon$. 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 will often discuss our results in terms of pumping rate $\mathcal{P}$. The reader has just to bear always in mind that $\mathcal{P}\sim 1/\sqrt{T}$..

II Numerical algorithm

We have studied systems with $N$ complex variables $a_{k}=|A_{k}|e^{i\phi_{k}}$ interacting with Hamiltonian in Eq. (5).

II.0.1 Exchange Monte Carlo

The sampling of the probability distribution in Eq. (9) was done by means of a Parallel Tempering Monte Carlo algorithm (PT) Hukushima96 . 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_{\vec{\i}}$. 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^{-\beta_{i}E}$.

We have run $K$ independent simulations at temperatures $T_{i}=\beta_{i}^{-1}$, from $T_{0}=T_{\textrm{min}}$ to $T_{N_{PT}}=T_{\textrm{max}}$. Each $64$ steps of the local Metropolis algorithm an exchange of configurations among simulations running at neighbouring temperatures is proposed. That is, if the Metropolis algorithm for ${\bf a}_{i}$ runs with $\beta_{i}$ and that of ${\bf a}_{j}$ with $\beta_{j}$ one tries the following attempt: $\{({\bf a}_{i},\beta_{i});({\bf a}_{j},\beta_{j})\}\Longrightarrow({\bf a}_{j},\beta_{i});({\bf a}_{i},\beta_{j})$, which is accepted with probability $p_{\textrm{swap}}$:

[TABLE]

The replica exchange update is proposed sequentially for all pairs of neighbouring temperatures $\beta_{i}$ and $\beta_{i+1}$. For all simulations the $N_{PT}$ temperatures where taken with a linear spacing in $\beta$, i.e. $\beta_{i+1}=\beta_{i}+\Delta\beta$.

II.0.2 The Monte Carlo update under spherical constraint

In the local Metropolis algorithm the configuration of complex “spins” $(a_{1},\ldots,a_{N})$ is updated with the requirement of keeping $\sum_{k}|a_{k}|^{2}=\textrm{const}$. In order to fulfill such constraint and the detailed balance each update is realized choosing at random two spins $a_{i}=A_{i}e^{i\phi_{i}}$ and $a_{j}=A_{j}e^{i\phi_{j}}$ and proposing an update to randomly chosen values of $a_{i}^{\prime}$ and $a_{j}^{\prime}$ such that:

[TABLE]

This is simply achieved by extracting three random numbers $\phi_{i}^{\prime}$, $\phi_{j}^{\prime}$ and $\theta$ with uniform probability in the interval $[0,2\pi]$ and proposing the four simultaneous updates $\phi_{i}\rightarrow\phi_{i}^{\prime}$, $\phi_{j}\rightarrow\phi_{j}^{\prime}$, $A_{i}\rightarrow A_{i}^{\prime}=C\cos(\theta)$ and $A_{j}\rightarrow A_{j}^{\prime}=C\sin(\theta)$. We have simulated systems at four different sizes: $N=32,48,64,102$. For the parallel tempering we used $N_{\rm PT}=32$ temperatures for all sizes.

II.0.3 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_{4}^{(i,j)}=\mathcal{O}(N^{2})$ quadruplettes, i. e.,

[TABLE]

The calculation of each $\Delta E_{\kappa}$ is realized in parallel on a distinct kernel on GPU.

Moreover, the execution of the $K$ evolutions at different temperature 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 K$20$, achieving an overall speedup, with respect to the sequential code on CPU, up to a factor $8$.

III Breaking of equipartition

Let us first compare the dependence on the pumping rate $\mathcal{P}$ of the emission spectrum measured in experiments (panel a. of Fig. 1, data taken from Cao98 ) with the emission spectrum in our numerical simulations averaged over the random coupling distribution (panel b.): the increase of $\mathcal{P}$ is accompanied by the same kind of narrowing. We stress that the symmetry of the spectral distribution around the central wave-length is due to the fact that the gain is taken as constant in this case, cf. Eq. (5), whereas in reality a material dependent, non-constant, gain curve $g(\lambda)$ is present.

The equipartition breaking obtained by raising $\mathcal{P}$ is, then, even clearer if we consider the behaviour 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 $\mathcal{P}$ a pattern of random sharp peaks emerges in the spectrum at high $\mathcal{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 non-linear systems Livi85 ; Cretegny98 we know that the degree of equipartition in the spectrum can be measured by the so-called spectral entropy:

[TABLE]

where $\hat{\mathcal{I}}_{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 Livi85 ; Cretegny98 :

[TABLE]

where $n_{\text{eff}}=1$ for perfect equipartition and $n_{\text{eff}}=1/N$ when the total energy is concentrated in a single mode Livi85 ; Cretegny98 . For our system the behaviour of its average, $\overline{n}_{\textrm{eff}}$, as a function of the pumping rate $\mathcal{P}$, is reported in Fig. 2, panel $a)$ [data are plotted as a function of $\mathcal{P}^{-2}=T$]. At low pumping rate we find $\overline{n}_{\textrm{eff}}\approx 1$, which signals a good degree of equipartition for all the sizes $N$ studied. On the contrary, by increasing $\mathcal{P}$ we find that $\overline{n}_{\textrm{eff}}$ rapidly decreases and the decrease depends on $N$. We identify as the pumping value at which $\overline{n}_{\textrm{eff}}$ starts decreasing below $1$ with a size-dependent lasing threshold $\mathcal{P}_{c}(N)=\sqrt{T_{c}(N)}$. The figure clearly shows that the larger the $N$, the steeper the $\overline{n}_{\textrm{eff}}$ decrease.

A stronger indication that we are dealing with a first-order transition with respect to $\overline{n}_{\textrm{eff}}$ comes from the study of the Binder, $\mathcal{B}$, and the bimodality, $b$, parameters of the distribution $P(n_{\text{eff}})$, displayed in panels $b)$ and $c)$ of Fig. 2. The Binder parameter, which is simply defined as the ratio between the fourth and second moment of the distribution $P(n_{\text{eff}})$,

[TABLE]

with $\Delta n_{\text{eff}}=n_{\text{eff}}-\overline{n}_{\text{eff}}$, is 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 valued 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 behaviour, pointed out first by Binder Binder84 : it is a non-monotonic reversed bell behaviour, 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

[TABLE]

where $n$ is the number of data composing the histogram of the probability distribution, $\kappa$ 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 $\gamma$ is the skewness, defined as

[TABLE]

Indeed, what we find is precisely what Binder presents as the phenomenology of a first-order transition: the value $\mathcal{P}_{c}(N)$ at which $\mathcal{B}$ is mostly far from a Gaussian corresponds to the peak of the bimodality indicator $b$, as well as to the breaking point of $\overline{n}_{\text{eff}}$. Eventually, in panel $d)$ of Fig. 2 we compare the spectral profiles obtained for pumping rates either well below or well above the lasing threshold. The evidence that $\overline{n}_{\text{eff}}$ behaves as the order parameter is very important for the purpose of the whole discussion. The main point of presenting here data on equipartition breaking is to compare them with the standard way to detect ergodicity breaking in disordered system, 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 the information that also the “order parameter” for equipartition breaking behaves as the parameter of a first-order transition: this puts on solids 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) Derrida81 ; Kirkpatrick89 ; Lubchenko12 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 Mezard87 . 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  Mezard87 . This is the typical behaviour 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 Sherrington75 , where the probability distribution of the overlap $P(q)$ is expected to take a nontrivial shape, different from a bimodal one, even for $N\to\infty$. 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 $4$-phasor random laser model displays these RFOT features at the critical pumping power. We stress that

Let us start from the specific heat anomaly. In panels $a)$ and $b)$ of Fig. 3 are shown the specific heat $C_{V}$ curves for different sizes of the system, where $C_{V}$ per phasor is

[TABLE]

The $\langle\phantom{\sharp}\rangle$ and $\overline{{\phantom{\sharp\sharp}}}$ 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 quadruplettes (see Methods for details). At all sizes $N$, $C_{V}/N$ has a characteristic non-monotonic behaviour, with the position of the peak depending on $N$. We identify this point as a finite size “critical temperature” $T_{c}(N)=1/\mathcal{P}_{c}^{2}$.

The good collapse of the curves at different $N$, shown in the inset of both panel $a)$ and $b)$ of Fig. 3, 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 $\tau\sim N^{-1/2}$, where $\tau=T/T_{c}-1$: the scaling exponent $1/\nu=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 Derrida81 . On the contrary, the scaling region for our 4-phasor mode-locked model follows the behaviour $\tau\sim N^{-3/2}$, with an exponent $1/\nu=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 pre-asymptotic 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 Mezard87 . Labeling with the Greek indices $\alpha$ and $\gamma$ two replicas, we have initially measured (cf. Methods) the phasor overlap $q^{\alpha\gamma}$

[TABLE]

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 non-trivial distribution in presence of a glass phase Mezard87 . In panels $c)$ and $d)$ of Fig. 3 the phasor overlap distribution $P(q)$ is shown, respectively for the randomly diluted 4-phasor model [panel $c)$] and the mode-locked model [panel $d)$] at different values of $\mathcal{P}$. In both cases, in the low-$\mathcal{P}$ ergodic phase, it turns out to be a symmetric Gaussian. Then, for $\mathcal{P}>\mathcal{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 quadruplettes overlap probability distribution $\mathcal{P}(\mathcal{Q})$ (see Methods). The quadruplette overlap $\mathcal{Q}^{\alpha\beta}$ between the energy stored in the same set of four coupled modes

[TABLE]

at equilibrium in state $\alpha$ and in state $\gamma$ reads

[TABLE]

where $\kappa$ runs over the ordered list $\{\vec{\i}\}$ of all $N_{4}$ non-zero four-body couplings. While the phasor overlap is computed as the average over $N$ correlated random variables, i.e. the local phasors overlaps $\bar{a}_{k}^{\alpha}a_{k}^{\gamma}$, the quadruplette overlap $\mathcal{Q}_{\mu}^{\alpha\gamma}$ is the average of $N_{4}=\mathcal{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 $\mathcal{P}$ the quadruplettes overlap has a very peaked distribution at $\mathcal{Q}\simeq 0$, for $\mathcal{P}\approx\mathcal{P}_{c}$ we find a clear signature of a secondary peak at $\mathcal{Q}>0$. In the inset of panel $a)$ of Fig. 4 the corresponding multimodality parameter $b$.

The definition of $b$ is the same given in Eq. (18) above, where skewness $\gamma$ and curtosis $\kappa$ are now defined as follows:

[TABLE]

with $\Delta q=q-\overline{\langle q\rangle}$. We find that the region where the parameter $b$ signals a bimodal distribution of the overlap is precisely the interval of pumping rates around $\mathcal{P}_{c}$. For the study of the plaquette overlap $\mathcal{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 Trebino94 ; Iaconis98 applies only to high-power directional impulses: unfortunately this is not the operating regime of standard random lasers Wiersma08 ; Leonetti13a . We cannot therefore rely on observables which require the measurement of phases. The phasor and quadruplette overlaps defined in Eqns. (20) and (22) depend on phases, so they are not quantities that can be measured in current experiments. A further overlap $\mathcal{C}_{\alpha\beta}$ between intensity fluctuations $\Delta_{k}\propto\mathcal{I}_{k}-\langle\mathcal{I}_{k}\rangle$ can, then, be introduced:

[TABLE]

where $\mathcal{I}_{k}=A_{k}^{2}$. At the mean-field narrow-band level the Intensity Fluctuation Overlap (IFO) $\mathcal{C}_{\alpha\beta}$ is proved to be in a one-to-one correspondence with the standard overlap, i.e., $\mathcal{C}_{\alpha\gamma}\propto q_{\alpha\gamma}^{2}$ $\forall~{}\alpha,\gamma$ Antenucci15f . The advantage of $\mathcal{C}_{\alpha\beta}$ with respect to $q_{\alpha\beta}$ is that it can be measured in real random lasers Ghofraniha15 ; Gomes16 ; Pincheira16 ; Supratim16 ; Tommasi16 ; Lopez18 ; Tommasi18 . In panel $b)$ of Fig. 4 the distribution $\mathcal{P}(\mathcal{C})$ is shown for $N=64$ at four different values of the pumping rate $\mathcal{P}$, two above and two below the critical $\mathcal{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 $\mathcal{B}$ dependence on $\mathcal{P}$ which reveals how $\mathcal{P}(\mathcal{C})$ brings the signature of a first-order transition: data are shown in the inset of panel $b)$ of Fig. 4. The behaviour of $\mathcal{B}$ as a function of $\mathcal{P}$ (plotted as $\mathcal{B}$ vs $\mathcal{P}^{-2}=T$ to have a clearer term of comparison with the literature) is the one characteristic of first-order transitions Binder84 . This behaviour 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 non-linear 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 a 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 Kirkpatrick89 .

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 equipartition - typical of ergodicity breaking in ordered systems with non-linear interactions Fermi55 ; Livi85 ; Chaos05 - and the breaking of ergodicity as described within the paradigm of replica-symmetry-breaking in disordered systems Mezard87 .

The mode-locked $4$-phasor model is the first example of a system where ergodicity breaking manifest 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 new phenomenon is due to the following two salient features of the model: the $\mathcal{O}(N)$ dilution of the non-linear couplings with respect to the corresponding mean-field model of the p-spin class, induced by the Frequency Matching Conditionand, and the local unboundedness of the light-modes 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 random laser presented here it 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 suggest a possible way to overcome the intrinsic difficulty usually encountered in the measure of $q$, which is not a single-experiment 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 jergon of disordered systems one needs more replicas of the same system. As we have shown, the occurrence of a non-trivial $P(q)$ is simultaneous with the loss of spectral equipartition, and one can, thus, simply take advantage of latter to detect the ergodicity breaking glass transition.

In conclusion, as weird as it can 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.

VI Ackowledgements

The authors thank D. Ancora, G. Benettin, L. Biferale, A. Crisanti, G. Parisi, A. Ponno and A. Vulpiani for useful discussions. The research leading to these results has received funding from the Italian Ministry of Education, University and Research under the PRIN2015 program, grant code 2015K7KK8L-005 and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, project LoTGlasSy, Grant Agreement No. 694925. G.G. acknowledges the financial support of the Simons Foundation (Grant No. 454949, Giorgio Parisi).

VII Appendix

VII.1 How to build the Mode-Locked graph

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 Angelani06a ; Conti11 ; Antenucci15a ; Antenucci15f , in particular when the non-linear interactions is the one of Eq. (5). Operatively, it is easier to describe the structure of the mode-locked interaction network as a bipartite graph where interaction nodes $\mathcal{J}_{\mu}$ 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 $4$ variable nodes and is defined by the ordered list of their indices $\mathcal{J}_{\mu}(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 non-equivalent permutations. The steps to generate the mode-locked graph are as follows:

a virtual complete graph with $\binom{N}{4}$ interaction nodes is generated; 2. 2.

for each interaction node the three non equivalent index permutations in Eq. (5) should be considered: $\mathcal{P}_{\mu}(i_{1},i_{2},i_{3},i_{4})$, $\mathcal{P}_{\mu}(i_{2},i_{1},i_{3},i_{4})$ and $\mathcal{P}_{\mu}(i_{1},i_{2},i_{4},i_{3})$, even though, for the index order $i_{1}<i_{2}<i_{3}<i_{4}$ only $\mathcal{P}_{\mu}(i_{2},i_{1},i_{3},i_{4})$ can be satisfied at most; each time that the mode frequencies (or indeces) satisfy the condition (6), the corresponding interaction of the virtual graph is added to the real graph. 3. 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)$ of all interacting quadruplets Marruzzo18 . Operatively, for a system with $N$ complex variables we have drawn a bipartite graph with a number of interactions scaling as $\mathcal{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 Tab. 1.

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 non-local in the basis of normal modes. The only quantity which depends on the real-space dimensionality is the spatial overlap between the normals modes, the information about which is stored in the disordered coefficients $J_{\vec{\i}}$ of the Hamiltonian in Eq. (5) as Conti11 ; Antenucci15a :

[TABLE]

where $\chi_{\vec{\nu}}^{(3)}$ is the non-linear susceptibility of the system. We do not consider contributions from the first non-linear term $\chi^{(2)}$ in the polarization expansion because of the relatively limited band width 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, that would read

[TABLE]

$g_{k}=g(\omega_{k})$ being the net gain profile (gain taken away losses due to the open cavity) that can be expressed as a function of $\chi^{(1)}(\omega,\bm{r})$. In our model though, where gain saturation is expressed by the global constraint on the overall intensity distributed in the system

[TABLE]

and considering for simplicity a constant net gain profile $g_{k}\simeq g$, this contribution amounts to a constant $\mathcal{H}_{2}=-g$, with no role in the dynamics.

In random lasers couplings as expressed in Eq. (25) are, in general, disordered because modes display different spatial shape and extension Conti08a ; Fallert09 . The constituents of the integrals in Eq. (25) are very difficult to calculate from first principles. The only specific form of the non-linear susceptibility has been computed by Lamb Lamb64 ; SargentIII78 for few-modes ordered lasers and no analogue study for RLs has been performed so far, to our knowledge. Integrals like Eq. (25) 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 Zaitsev10 . 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_{\vec{\i}}$:

[TABLE]

The scaling of variance with $N$, $\langle J^{2}\rangle\sim N^{-2}$, guarantees energy extensivity. There are very many couplings but each one is vanishingly small. Macroscopic phenomena like lasing occur because the system undergoes a transition to a collective behavior. Even in presence of randomness cooperativity is the leading mechanism. We further stress that, from the perspective of probing RFOT, considering correlated $J$’s leads to qualitatively analogue phase diagram as it is well known in spin-glass systems such us, e.g., the Random Orthogonal Model Marinari94a ; Parisi95b .

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\to\infty$. That is why we expect the mean-field glass transition scenario drawn in Antenucci15a ; Antenucci15f to be quite robust for the mode-locked $p$-phasor, even if the narrow-band hypothesis Gordon02 is removed.

Last but not least, the non-locality of interactions between light modes also guarantees that phenomena like energy localization, a pathology of sparse networks Antenucci15c ; Antenucci15d , are avoided.

VII.2 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 MC dynamics. Different istances of the PT dynamics characterized by the same set of quenched couplings $J_{\vec{\i}}$ and the same interaction network between modes are different replicas.

For $N=32$, $N=48$ and $N=64$ for each disorder instance we simulated $4$ replicas, which gave us the availability of $6$ independent values of $q^{\alpha\beta}$: $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^{\alpha\beta}$ comparing replicas at the same iteration of the PT dynamics, each $640$ iterations. Since the distribution $P(q)$ is not self averaging Mezard87 , for each size of the system we have sampled the equilibrium measure for $N_{\text{sample}}\approx 100$ istances of the disorder.

It is useful to clarify also how we measured in practice all thermal averages indicated with angular brackets, $\langle\mathcal{O}[{\bf A}]\rangle$. Once that the system reaches equilibrium at a given temperature/power, at a time in Monte Carlo sweep units of $N_{\rm term}<N_{\text{sweep}}/4$, thermal averages were measured as time averages along the dynamics, along the second half of each run:

[TABLE]

VIII References