The glassy random laser: replica symmetry breaking in the intensity fluctuations of emission spectra

The behavior of a newly introduced overlap parameter, measuring the correlation between intensity fluctuations of waves in random media, is analyzed in different physical regimes, with varying amount of disorder and non-linearity. This order parameter allows to identify the laser transition in random media and describes its possible glassy nature in terms of emission spectra data, the only data so far accessible in random laser measurements. The theoretical analysis is performed in terms of the complex spherical spin-glass model, a statistical mechanical model describing the onset and the behavior of random lasers in open cavities. Replica Symmetry Breaking theory allows to discern different kinds of randomness in the high pumping regime, including the most complex and intriguing glassy randomness. The outcome of the theoretical study is, eventually, compared to recent intensity fluctuation overlap measurements demonstrating the validity of the theory and providing a straightforward interpretation of qualitatively different spectral behaviors in different random lasers.

The behavior of a newly introduced overlap parameter, measuring the correlation between intensity fluctuations of waves in random media, is analyzed in different physical regimes, with varying amount of disorder and non-linearity. This order parameter allows to identify the laser transition in random media and describes its possible glassy nature in terms of emission spectra data, the only data so far accessible in random laser measurements. The theoretical analysis is performed in terms of the complex spherical spin-glass model, a statistical mechanical model describing the onset and the behavior of random lasers in open cavities. Replica Symmetry Breaking theory allows to discern different kinds of randomness in the high pumping regime, including the most complex and intriguing glassy randomness. The outcome of the theoretical study is, eventually, compared to recent intensity fluctuation overlap measurements demonstrating the validity of the theory and providing a straightforward interpretation of qualitatively different spectral behaviors in different random lasers.
Light amplification and propagation through random media have attracted much attention in recent years, with present-day applications to, e.g., speckle-free imaging and biomedical diagnostics 1 , chip-based spectrometers [2][3][4] , laser paints 5 and cryptography 6 . Whatever the amplifying medium, ordered or random, in a closed or in an open cavity, two are the basic ingredients to produce laser in any optically active system: amplification and feedback. In closed cavities the electromagnetic modes straightforwardly depend on the cavity geometry. In cavity-less random media, instead, some kind of modes are established by spontaneous emission and are localized in closed photonic trajectories by means of multiple scattering. Indeed, the phenomenon of amplified spontaneous emission (ASE) can occur even in systems without any optical cavity, whose fluorescence spectrum is simply determined by the gain curve of the active medium [7][8][9][10][11][12] . When, because of an external pumping, the multiple-scattering feedback process becomes strong, amplification by stimulated emission is established in the random medium and we have a Random Laser (RL) 13 . The feedback is, here, associated to the existence of well-defined long-lived modes, characterized by a definite frequency and a spatial pattern of the electromagnetic field inside the material. Modes are expressed as slow amplitude contributions to the electromagnetic field expansion in terms of spatial mode eigenvectors ( ) The complex amplitudes a k (t) of these slow modes turn out to be the fundamental degree of freedom in the statistical mechanical modeling of interacting modes 14,15 , while the irregularity of their spatial profiles results into quenched disordered mode interactions. By quenched we mean that the interaction strengths are time independent 16 , as it occurs, in practice, when they change on time-scales much longer then the inverse strength of the nonlinear interaction coupling with respect to the off-diagonal linear coupling α = / = / ∈ , . In a closed cavity the linear dumping is absent and it corresponds to α = 1. In a standard semiclassical approach, the field is expressed in the slow amplitude basis, equation (1), where each mode displays a determined frequency. The lifetimes of these modes are assumed to be much longer than the characteristic times of population inversion and amplification processes, so that the atomic variables can be adiabatically removed and result in an effective interaction between the electromagnetic modes. The nonlinear couplings are, indeed, nonzero only for the terms ⁎ ⁎ a a a a j k l m that meet the frequency matching condition [49][50][51] , γ being the finite linewidth of the modes. The mean-field approximation of the model equation (2) is exact when the probability distribution of the couplings is the same for all the mode couples (j, k) and tetrads (j, k, l, m). This is true, e.g., when mode extensions scale with the volume occupied by the active medium and their spectrum has a narrow-bandwidth around some given central frequency ω 0 , i.e., ω ω γ | − | < j 0 , ∀ j, so that the frequency matching condition, cf. equation (3), always holds.

Results
The Random Laser Transition. Given the quenched randomness of the J's, any observable depends on the particular realization of the disorder. Thus, the relevant quantity is the disorder averaged free energy β = − / F Z ln J , where the overline denotes the average over the distribution of quenched disordered couplings. This can be analytically evaluated using the replica method 27,37 , as reported in the Methods. In this procedure the evaluation of the relevant thermodynamic quantities is achieved considering n identical replicas (i.e., copies) of the system that act as probes exploring the multi-state phase space of the system. Further on, evaluating the distance between the replicas in terms of their similarity, termed overlap, one can retrieve the physical overlaps of the thermodynamic states 35 .
In the complex amplitude spherical model, equation (2), the order parameter of the replica theory turns out to be given by the overlap between amplitudes of replica a and replica b: where ε E = /N . This overlap Q identifies the onset of a RL regime at a given critical value of the pumping, or, otherwise, at a critical temperature at fixed pumping. We notice that the latter behavior is in qualitative agreement with the experimental results of refs 52-54 where random lasing appears to occur decreasing temperature, besides increasing pumping.
Any nontrivial structure of the values taken by Q ab implies that identical copies of the system, with the same interaction network and submitted to the same thermodynamic conditions, show different sets of values for microscopic observables at equilibrium and the ergodicity is broken in distinct equivalent states.
To our knowledge, from an experimental point of view, no phase correlation measurements, required for the evaluation of the complex amplitudes a k a and, consequently, of Q ab , is available so far in random media. Only magnitudes a k are measured and not their phases φ = ( ) a arg k k . The experimental reconstruction of the distribution of the values of equation (4) is, thus, unfeasible.

Real Replicas.
In recent experiments 44 , shot-to-shot fluctuations of intensity spectra in an amorphous solid RL, a functionalized thiophene-based oligomer named thienyl-S,S-dioxide quinquethiophene (T5COx), are measured and analyzed. Since the sample remains under identical experimental conditions shot after shot, n different shots of RL emission correspond to n real replicas and one can measure the overlap between intensity fluctuations of two real replicas. In these experiments, the set of the activated modes emitting after the shot = ,…, n a 1 , whose available coarse-grained degree of freedom is the intensity ( ) = I k a k a a 2 , is observed to change from shot to shot. When, during a single shot of the pumping source, the number of stimulated emission processes taking place is very large, the configurations of the mode dynamics can be considered as pertaining to a thermodynamic state. In terms of the photonic bomb language of Letohkov 10 , e.g., this is a situation in which the typical amplification time is much shorter than the photons lifetime inside the medium, i.e., of the lifetime of stochastic resonators supporting the localized optical modes. The possible observation of numerous different states from shot to shot is, consequently, an evidence of a thermodynamic phase described by a corrugated free energy landscape composed of many valleys separated by barriers.

Intensity Fluctuation Overlap (IFO
the intensity fluctuation of shot a around the average profile, one can define the overlap between the normalized intensity fluctuation of shots a and b as 44 : where the index k denotes now the frequency, i.e., the experimental accessible equivalent of a mode index, depending on the spectral resolution. The overlap is measured between the fluctuations of intensity, rather than the straight intensities, to exclude the effects due to the amplified spontaneous emission. From n measured spectra one can calculate the ( − )/ n n 1 2 values of the IFO  ab exp and its distribution exp e xp can be computed by repeated spectral measurements acquired on different samples. By different samples we, actually, mean different realizations of the microscopic disordered realization of scatterers positions as faced by the incoming pumping light beam. More precisely, one can realize a different realization by turning the material sample or, if the beam section is smaller than the random medium, by illuminating a different region with the pump laser spot.
If the variations of the normalization factors ∑ ∆ ( ) k k a 2 in Eq. (5) are neglected with respect to fluctuations ∆ ( ) k a , in the 2 + 4 complex amplitude spin-glass model given by equation (2), the matrix is the model equivalent of the IFO, up to an overall sign. Indeed, equation (6), defined in the dominion [0, 1], holds with the prescription that . To compare with experimental results we will use symmetrized  ( ) P ,  ∈ − [ 1:1], without any loss of generality. (6) can be carried out using the replicated action derived in the Methods, cf. equation (14). This leads to the following relationship between the IFO and the standard order parameters:

IFO vs. standard overlap relationship. The average in equation
where Q ab is defined in equation (4), and is the parameter of global coherence (cf. Methods). According to equation (7), if a RSB occurs in the standard overlap Q ab it propagates with the same structure to the IFO  ab . We, thus, have a theoretical well funded tool to detect RSB in experimental data. We stress that this analysis could have not been possible in XY models with quenched amplitudes considered in previous works [30][31][32] because there the intensities of the modes are kept fixed during the mode dynamics.
In equation (7) we have considered the most general case in which a high pumping regime can display both a global coherence ( ≠ m 0) and a multi-state non-trivial structure for the amplitude configurations ( ≠ Q 0 ab ). This mixing physically occurs for a degree of disorder R J next to the tolerance value beyond which standard mode locking (SML) breaks down, leaving place to glassy random lasing. This is displayed in the phase diagrams in the central panels of the triptych Figs 1 and 2, as the boundary lines between SML ( ≠ m 0) and glassy random laser (m = 0 but ≠ Q 0 ab ) at large . The low pumping regime is replica symmetric for any R J , with m = 0 and Q ab = 0 for ≠ a b 14 , implying a Dirac delta probability  ( ) P , peaked in zero, both in the incoherent wave (IW) and in the phase locking wave (PLW) in Figs 1 and 2 (cf. also Methods).
Replica Symmetric standard mode-locking laser. For weak disorder at high pumping , for every , with ≠ m 0, holds between the overlap and the (replica independent) global coherence parameter. In other words, the high pumping laser regime is replica symmetric, as well, and the  ( ) P is a Dirac delta function in zero, once again. The laser regime for negligible disorder (  R 0 J ) corresponds to a standard mode-locking laser in an ordered cavity 49,50 with spectral resonances equispaced in wavelength and smoothly distributed in intensity. For weak, though not very small, disorder, as, e.g., in the open cavity phase diagram of Fig. 2 , the laser regime could also correspond to a deterministic (i.e., non-glassy) random laser, with spectral resonances indeed at random wavelengths with random intensities, yet always with the same pattern for each experiment on the same sample under the same external conditions. From a classical statistical mechanics point of view these two cases are equivalent (in terms of the order parameters equations (4,8)) and they are referred simply as SML in the Scientific RepoRts | 5:16792 | DOI: 10.1038/srep16792 phase diagrams. We will come back to this kind of random lasers when discussing known experimental realizations.
Remarkably, though in terms of the parameter m the standard mode locked regime is clearly different from the fluorescence regime, and so is ( ) P Q , cf. left panels − a d of Figs 1 and 2 and refs 14,15, the IFO distribution does not change below and above the standard mode-locking transition. This has been observed in preliminary measurements on a Q-switched pulsed Nd-Yag standard laser in ref. 44. Indeed, the overlap of the model in equation (6) is between local fluctuations of intensity on different replicas, so a global ordering is invariably taken away (cf. Methods).
Onset of RSB across the random laser transition. For strong disorder the distribution of the coupling J's yields a non-negligible amount of both positive and negative values, inducing frustration in the modes interaction. Take, for instance and simplicity, three complex amplitudes a, pairwise interacting on a triangle. Each of the three bonds connecting the three modes can randomly acquire positive or negative values. Modes connected by positive bonds will tend to align (in the complex plane), whereas modes connected by negative bonds will tend to counter-align. If, e.g., in the triangle two bonds are positive and one negative, no single configuration of modes alignments will satisfy all interactions, minimizing the Hamiltonian. The system will then settle into one of different degenerate configurations with the lowest realizable energy. Frustration is, thus, the impossibility of finding a unique way of satisfying all bonds. Frustration is a necessary condition for the onset of glassiness. When the pumping increases above threshold,   > c , the replica symmetry is broken and the distribution of  ab becomes nontrivial, cf. panels f, g, h in  is invariant. In the right panels e to h the IFO and standard overlap distributions are shown for the RL transition: as  increases, we show that the low  solution is replica symmetric (e), while above threshold it becomes discontinuously 1RSB (f, g, h). In a closed cavity situation, α = 1, ( ) P Q is discontinuous at the standard mode-locking laser transition, while  ( ) P is unaffected. In this case the transition itself is discontinuous in the thermodynamic sense: the internal energy 14,50  , alternatively, the ( ) P Q , and similarly  ( ) P , changes in a nontrivial way: two different values, a zero and a nonzero one, are possible as the pumping is increased. In this situation the RL regime is one step RSB (1RSB) and the transition is a so-called random first order (RFOT) in glassy physics terming 16 . In the RFOT scenario the static (ideal) glass transition is preceded by a glassy dynamic arrest (drawn as a dashed line in the central panel of Fig. 1) 14,15 : a photonic system in this case should show the typical two-step dynamical relaxation for the time correlation function of light modes, in the same universality class of the mode-coupling theory for structural glasses. In this kind of transition there is no latent heat 15 , yet a new value for the overlap discontinuously appears at the transition, cf. right panel of Fig. 1.
In the cavity-less scenario α = 0.4, instead, ( ) P Q is continuous at the ordered ML transition. Indeed, a nonzero value for m increases continuously from zero as it can be observed looking at the peaks of ( ) P Q in panels a and b of Fig. 2, where = ± / Q m 2 2 . As in the closed cavity scenario, the  ( ) P of the SML does not change across the threshold. At the onset of the RL regime, illustrated in Fig. 2 is rather meaningful. At and just above the threshold,  ( ) P displays a continuous part between the central peak in  = 0 and the two side peaks, as displayed in panel h of Fig. 2. Here, the transition is thermodynamically continuous with a RL regime that is of the so-called full replica symmetry breaking (FRSB) kind, associated with a free energy landscape composed by a fractal hierarchy of valleys. As the pumping increases, the regime becomes 1 + FRSB, a combination of 1RSB and FRSB solutions, with both a continuous and a discontinuous contribution to the probability distribution, cf. panel g in Fig. 2. The continuous parts in the ( ) P Q and  ( ) P depend on the influence of the off-diagonal damping term in J ij in equation (2). For high enough pumping, well-above the does not change below and above threshold. In the right panels e to h the IFO and standard overlap distributions are shown for the RL transition. As  increases we show that the low optical power solution is replica symmetric (e), soon above threshold the solution is FRSB (f), further increasing  the solution becomes 1 + FRSB (g) and, eventually, for large pumping it is 1RSB (h). The transition is continuous in the order parameters  ( , ) P Q .
threshold, the non-linear term eventually becomes dominant 15 and the solution, cf. panel h in Fig. 2, eventually becomes 1RSB, as in the closed cavity case, cf. panels f, g, h of Fig. 1.
In the RL experiment of ref. 44 the distribution  ( ) P exp , with  exp defined in equation (5) is peaked in zero at low pumping, while it becomes nontrivial with a triple and, eventually, double peaked shape as the lasing threshold is overcome. Although in comparison with the theoretical predictions for N → ∞ the peaks of  ( ) P exp are smeared by noise effects and finite modes' number effects, in all regimes Fig. 3 we display a comparison between the analytic IFO distribution computed in our 2 + 4 complex amplitude spin-glass model, cf. equation (2), in an open cavity and the experimental measurements of  exp in refs 44,55.

Discussion
In this work we provide the theoretical analytical background for a recently introduced order parameter 44 that allows to probe the phenomenon known as replica symmetry breaking in random lasers by means of experimentally accessible observables. These are shot-to-shot intensity fluctuations and the order parameter is the distribution of the values of the overlap between intensity fluctuations in different shots, as analytically defined in equation (6). Replica symmetry breaking is a known property occurring in mean-field glasses, spin-glasses and hard optimization problems. The parameters of the theory have never been measured, though, in any real system in these fields. In, particular, no measurement of the overlap and its distribution has been provided. The only experimental measure, so far, of a quantity possibly related to the standard RSB overlap been recently carried out on a photonic system. The system is an amplifying and scattering random medium, the T5COx 44 displaying random lasing at high pumping. The parameter is the distribution of the shot-to-shot intensity fluctuations overlap (IFO). In the framework of a recently introduced general statistical mechanics theory of random photonic systems 4 , in equation (21) we give here an analytic proof of the relationship between the IFO and the standard overlap, equation (4), and we provide measurable predictions for its behavior in both ordered and random lasing systems below and above threshold and, furthermore, both in the cases of discontinuous and continuous transitions to the laser regime at the threshold. In particular, the transition in the probability IFO distribution of a random laser is shown to be discontinuous (cf. Fig. 1) for closed (or controllable, limited open) cavities while it becomes continuous (cf. Fig. 2) for highly open cavity nonlinear wave systems. In the cavity-less case, where experimental measurements are available in at least one case, in Fig. 3 we compare theoretical and experimental behavior of the distributions of the IFO  from low to high pumping.
According to our results, a RSB is to be expected only, though not always, in random lasers whose random configurations of scatterers are fixed, i.e. quenched, for all analyzed shots. That is, the dynamics of their positions evolves on time-scales much longer than the whole experiment and real replicas can be realized. This is the experimental case of the solid/powder samples of random lasers as GaAs R 1 1 J and for increasing pumping. Vertical lines represent Dirac's deltas, whose height is the probability of the argument value. Different regimes are represented from fluorescence to large pumping random lasing. They are chosen along the dotted line in Fig. 2 . Form left to right the first distribution is at point e in Fig. 2, the second between f and g, the third one between g and h and the following above h. In the bottom row the same regimes are reproduced in the IFO distribution experimentally measured and reported in refs 44,55 in an amorphous solid oligomeric random laser, T5COx.
Displaying fixed scatterers to realize real replicas is not sufficient to yield a glassy random laser, though. Indeed, as we previously discussed in the Results section, frustration is also necessary. A notable example of a frustration-less solid random laser might be porous gallium phosphide (GaP) filled with a solution of Rhodamine and methanol 61,62 , in which spectral fluctuations are reported to be minimal and the structure of the resonances, though random, appears to be reproducible from shot to shot. IFO measurements might yield, in this case, an ordered-like  ( ) P , peaked in zero both below and above threshold, as in the phase reported as SML in the diagram of Fig. 2 for On the other hand, experiments on optically active random media whose scatterer particles sensitively move between subsequent shots in a single experiment, as in liquid solutions of Rhodamine and methanol with particles of Titanium oxide 63 , Zinc oxide 64 , pure Titania 65 , or colloidal CdSe quantum dots 66 could establish no real replicas. Not having the same quenched disorder in all shots might prevent the observation of RSB. The overlap between copies of systems with different realizations of the disordered couplings, indeed, is known to be replica symmetric, as it has been shown in models with continuous spherical variables 67 , of which our model in equation (2) is a generalization. Similarly to what happens in the ordered ML case, cf. left panels of Figs 1 and 2, in that case the occurrence of a trivial single peaked  ( ) P in  = 0 is expected, both below and above  c . Such a behavior has been observed in a liquid system of TiO 2 scattering nano-particle suspensions in solution of Rhodamine and methanol 44 .
Eventually, we would like to stress that, besides a rigorous interpretation of recent experimental results for random lasers in terms of replica theory, our results provide an exciting and easily available test of spin-glass theory properties in continuous systems without local magnitude constraints, as disordered photonic systems.

Methods
Replica Theory and Order Parameters. The most complicated system that we are considering in our theory is a random system with disordered mode couplings that possibly display a high pumping/ low temperature phase with ergodicity breaking and the occurrence of very many states. By "very many" we mean that their number scales with the size of the system, i.e. the number N of optically active modes. These states are not related by any simple relationship among them. That is, e.g., no simple Z 2 spin reversal symmetry occurs between states, as in the Ising model, nor ( ) SU 2 symmetry as in the XY model. In the complex glassy case, to probe the multi-state disordered thermodynamic phase, one, thus, considers n copies of the system with exactly the same set of disordered couplings, the J's, and evaluates the disorder averaged partition function Z J n of the replicated system. A continuation to real n is, then, taken to evaluate As the system size becomes sufficiently large, the free energy sample-to-sample fluctuations die out and the free energy, equation (9), becomes independent of disorder, i.e., it is self-averaging. For N → ∞ the physical value of the matrices follows from the extremization of the free energy functional. Because of the fact that the number of independent elements of an overlap matrix is ( − )/ n n 1 2 (taken away the diagonal) in the limit n → 0 the usual minimization of the thermodynamic potential actually becomes a maximization in the space of the overlap matrices. To maximize F, a non-trivial Ansatz on the structure of Q and R is necessary. Indeed, it can be shown 14 that the most intuitive replica symmetric solution, with Q ab and R ab independent of a and b, does not lead to a thermodynamically stable solution in the whole phase space: beyond the critical point, in the glassy phase, one must, hence, resort to spontaneous RSB. Following the Parisi scheme 27 the overlap matrices are, then, taken -step RSB matrix, with  → ∞ for a continuous full RSB (FRSB). These are block matrices where the number of inner blocks  + 1 corresponds to the number of hierarchical levels in the multi-state phase space.
Depending on the value of J 2,4 the solution of the RL model Eq. (2)   For weak disorder (low R J ) the global coherence σ τ , m is non-zero above the lasing threshold and must be included into the description. If disorder is strong, though, in the frozen glassy phase, the global coherence is null: = σ τ , m 0. Because it turns out that σ τ = 0 a b 14,15 , the integrals in the σ, τ space factorize and the IFO  ab defined in equation (6)  The replicated action  given in equation (14) is quadratic. Thus, using the Wick's theorem, for the averages in equation (19) we easily obtain where Q ab is defined in equation (4),  14,15 . Equation (21) is one of our main results and is discussed in the main text, cf. equation (7).