Universality class of the mode-locked glassy random laser

We thank both referees for a thourough reading of our manuscript and theor comments and observations. They certainly helped improving the presentation of our work. In the following we reply in detail to each comment of the referees.

Report 2 REFEREE: As previously said, the main problem I encountered while refereeing is about the scaling laws provided by the authors. They claim that the critical exponents for a standard second order critical point (and also for the REM) are ν = 2 and α = 2/3. However, I naively thought that they should be ν = 1/2 and α = 0 (as for standard scalar φ4 -theory). To justify their claim, the authors use a scaling argument that I find hard to understand. In particular, they approximate the fuctuations below the critical point by the distance of the minimum of the Landau free energy from the origin, while the curvature of the Landau free energy at the minimum should be considered instead. REPLY: The referee is absolutely right and we apologize for being sloppy in the presentation. We unproperly used the symbol $\nu$ for $\nu\times d$, $d$ being the dimension of the underlying short-range interacting lattice. Indeed, if the theory if a $\phi^4$ theory, in the limit of validity of the mean-field theory it is $\nu = 1/2$, and $\nu d$ = 2 because the upper ctitical dimension is 4. Furhermore, we were also unclear in distinguishing the theoretical mean-field value of the specific heat exponent $\alpha=0$ 0 from the numerical estimates obtained from finite size scaling, as, e.g., $\alpha\simeq 2/3$ in a previous study of the ML random laser [Ref 20], or $\alpha\simeq 0.52$ in the numerical simulations of the REM (who is known to be in the universality class of the \phi^4 mean-field theory). The scaling argument, actually, directly concerns the susceptibility and the order parameter behavior around criticality, that is, a relationship for the β and γ exponents: 2 β +γ = 2 (in general 2 β +γ = n/(n-2) for a φ^n theory). Because of the hyperscaling relation 2 β +γ = ν d, in a φ^4 theory at the upper critical dimension this corresponds to what we errenously termed ν (and also referred to as “the correlation length critical exponent” in at least one occasion) in the old version of the paper. We have been rewriting the section 4, and part of the conclusions, making explicit the role of β and γ exponents in the scaling argument. We introduce the term ν_{eff} to sometimes shorten 2 β +γ. In the text of the revised version we explicitly report that the α exponent (estimated from the finite size peak of the specific heat) suffers of strong pre-asymptotic effects, as far as the simulated sizes are small. This can be observed already in the simpler mean-field cases, as the REM, where the estimated 2β +γ ~ 2 is compatible with 2 but the estimate for $\alpha$ is non zero. This finite size pre-asymptotic effect for α is detected also in the mode-locked $4$-phasor model.

— REFEREE: I do not understand the last paragraph of the left column page 3, where it is written that the Fourier transform of ak should be roughly a Dirac delta function at ωk . Indeed, ak should be a slowly-varying amplitude, and thus its Fourier transform should only display low frequencies centered around 0, and not around ωk . I guess the authors discuss the Fourier transform of the electric field rather than the one of the slow amplitude ak(t), but maybe I am missing a point here. The referee is right also here. There is a mistake in the formula. When we perform the Fourier transform of the em field (3) on fast time scales (inverse of the frequencies of the normal modes waves oscillations), i.e., on a time window $2\tau$ large enough so that many oscillations take place ($1/\omega_k \ll \tau)$ we have
 \sum_k \int_{t-\tau}^{t+\tau} ds/2\pi e^{-\imath ω s} a_k(s) e^{\imath ω_k s} 
that in the slow amplitude approximation (a_k(t) does not change on time scale $\tau$) is approximately equal to 
 \sum_k a_k(t) \int_{t-\tau}^{t+\tau} ds /2\pi e^{-\imath (ω-ω_k) s} \simeq
 \sum_k a_k(t) \delta(ω-ω_k).
We are left with a slowly varying amplitude a_k(t) times a \delta (ω-ω_k). One may write, as well, a_k(t) = a(t,ω_k). We corrected the text accordingly. —
REFEREE: Second paragraph of page 4, I do not understand the end of the sentence ", and in the case of cavityless systems also compensate the leakages". I suspect one word is missing or replaced by another one. REPLY: We thank the referee to pint this out, we added a subject to the sentence:
“and in the case of cavityless systems also compensate the leakages.” becomes
“and in the case of cavityless systems energy pumping also compensate the leakages.” — REFEREE: Just above Section III, the authors say that "As already mentioned (and without loss of generality) we will consider the J's as real parameters". I could not find to which part of the paper the authors were referring. Does it mean that the authors assume to be in the purely dissipative regime? REPLY: Yes, we changed the sentence in “As already mentioned in section II when introducing the dissipative limit we will consider the J's as real parameters, without loss of generality.” —
REFEREE: At the beginning of the third paragraph page 5, the authors simplify the problem in the limit γ << \delta ω , and conclude that H2 is an additive constant to the Hamiltonian. I do not understand why, since the couplings are a priori different from one mode to another, and the spherical constraint cannot be factorised out. REPLY: The Referee is right also in this respect. If (16) holds we have N modes with N distinct frequencies and the H_2 is diagonal: \sum_J_k |a_k|^2. Our statement of additive constant only holds if J(k) is homogeneous in k. The diagonal part of the second order contribution to the Hamiltonian represents the gain profile of the optical random medium. As a working hypothesis we are assuming a uniform gain profile over the whole spectrum. 
We amended the revised text accordingly. —
REFEREE: Is there a reason why Eq. (19) is not just written
|ka - kb| = N - |ka - kb| if |ka - kb| \gtrsim N/2? Besides N/2 is not well defined for odd N , and so the integer part should be considered. REPLY: We agree that it is more immediate as the referee. Wa emended the text accordingly, including [N/2] (integer part of), in the condition (19). —
REFEREE: In Appendix B, I do not understand the argument for the choice of the scaling function [below Eq. (B3)]. Could you explain more? REPLY: The scaling function $\hat f(x)$ retraces the behavior of f_N(T) around T_c. That is, the linear term is not there because the maximum of the parabola is in $t_N=0$, by definition of $t_N=T-T_c(N)$. Both in $f_N$ and in $\hat f$ we might have added a third order term but we have tested that it becomes irrelevant in the critical region. —
REFEREE: At the end of Section IV, the authors conclude that the model belongs to a mean-field universality class because the value of 1/ν is consistent with this statement. Actually, the only rigorous statement that can be made is that the universality class of the model is not incompatible with a mean-field one. How do the authors conclude that the universality class is indeed a mean-field one? I think I am missing a point here. REPLY: We conclude that the universality class is compatible (i.e., non incompatible) with a mean-field universality class, though strong, possibly pre-asymptotic, finite size effects prevent us from saying what universality class. So we basically agree with the referee. Rather than “assess the mean-field nature” or “yield a mean-field universality class”, we have now limited ourselves to say that our best data, obtained with pbc, yield a picture that is compatible with a mean-field class of universality. 
We changed the text accordingly. For instance “we can assess the mean-field nature of the glass transition in the ML 4-phasor model.” becomes “we observe that the glass transition of the ML 4-phasor model is compatible with a mean-field transition”. —
REFEREE: Could you justify why the complexity for one step goes from N 2 to ln N via parallelisation on GPUs, as claimed in the Conclusion? REPLY: At each update proposal of a given couple of modes the new \Delta H is computed, that is all the contributions to the Hamiltonian involving the two modes are changed. These are O(N^2) and involve the product of the quadruplets mode amplitudes times the relative coupling constant.
If we split the computation of each one of these contribution on different kernels on a GPU, they are all computed in parallel, synchronously. If we increase N, the computation of the products yielding each quadruplet contribution to the energy does not scale with N because it is carried out in parallel. We are left with the operation of sum, only. This can be carried out summing - in parallel - couples of quadruplets contributions, then we sum the N^2/2 remaining numbers as couples, than the N^2/4 numbers left are summed in couples and so on until we reach a single number: the total sum. For instance, if N^2=16 it takes 4 steps, log_2 16, to sum the total energy. If N^2=64 it takes 6=\log_2 64 steps. In general it takes \sim \log N^2 = 2 \log N steps to sum O(N^2) terms. —
REFEREE: In the Conclusion, the authors claim that the RSB occurs at the lasing threshold. But from the data of Fig. 11, it seems that the the secondary peaks in the overlap distribution appear at much lower temperature (about 0.4) than the critical temperature Tc (about 0.6). Is it expected that the two temperatures coincide in the large N limit? REPLY: Numerically, at finite $N$, the critical region is around the temperature at which deviations from gaussianity occur in the P(q). Eventually, at low enough T, deep in the spin glass phase, side peaks arise but the transition already occurrs at higher T, when the P(q) at finite N is no longer a Gaussian. This is not easy to be detected by the eye and one has to rely on specific observables such as the Binder cumulant (basically the deviation of the kurtosys from 3). The analysis procedure was extensively presented in Ref [40] so that we did not report it here.
However, we feel that we have to clarify the point raised by the referee and we modified the text in the conclusions as follows: “We find clear evidence for the occurrence of a Replica Symmetry Breaking phase at low temperature. Studying the deviation of the overlap distributions from a Gaussian distribution by standard methods (e.g., the Binder cumulant), as performed in Ref. [40], the onset of such a spin-glass phase can be shown to occur at a temperature consistent with the laser threshold identified by FSS analysis of the specific heat peaks.”

Eventually, on request of the referee, we have 1- Elaborated the scaling arguments. 

2- Indicated the assumptions behind the model used eventually for the simulations (only a quartic coupling, and real parameters).

3- Provided the simulation details for the REM in a new Appendix, B.

4- Discussed in the Appendices A and B where the error bars are coming from in Figs. 6, 7 and 8. 5- Clarified the FSS protocol in the Appendix C (former B). 6- Clarified the argument to conclude that the model is compatible with a mean-field theory.

7- Corrected the various typos.

Report 1

REFEREE
1- Just above Eq. (10), the authors introduce the ratio P between the pump rate and the spontaneous emission rate. According to what the authors say, the spontaneous emission rate is proportional to the temperature of the bath T , while the pump rate is proportional to the energy of the system E = \epsilon N . Consequently, I would expect P \propto \epsilon/T . However, Eq. (10) states that 1/P 2 = T /\epsilon2 , which seems to be inconsistent with the deRnition. Can the authors clarify this point? REPLY 
The sentence the referee is referring to is, indeed, confusing. We have corrected it in the revised version of the paper. In fact, with the notation P it is meant the pumping rate, which is defined as the inverse of the effective ("photonic") temperature squared. 
We modified the text as:
“However, a stationary regime can be described as if the system is at equilibrium with an effective thermal bath, whose effective temperature (a ``photonic'' temperature) accounts both for the amount of energy ${\mathcal{E}}=\epsilon N$ stored into the system because of the external pumping and for the spontaneous emission rate. The latter is proportional to the kinetic energy of the atoms, e.~g., to the heat bath temperature $T$. Eventually, the external parameter driving the lasing transition turns out to be [34,35,50]
 T_{\rm photonic} = \frac{T}{\epsilon^2}.
One can also introduce the pumping rate $\mathcal{P}$ [42] as the inverse of the square root of this ratio:
 \mathcal P^2 = \beta \epsilon^2 = \beta_{\rm photonic} “ REFEREE
2- In Sec. V the authors present numerical results for the overlap distribution across the glass transition. They say they used N equilibrium configurations of replicas to accumulate statistics, but only in Appendix A the authors relate the quantity N to the Monte Carlo steps of the algorithm and no reference to the equation in the Appendix is present in the main text. I would suggest adding a reference to the Appendix so that the reader can have a better understanding of the meaning of . Moreover, neither in the caption of Fig. (9) nor of Fig. (10) there is an indication of the typical value of N in the numerical results presented. Is it possible to add it? REPLY
Certainly we are adding this information and we thank the referee to point this out. The number of configurations $\mathcal N$ actually used from our data can be evinced from the fifth column of tables \ref{tab2}-\ref{tab3} in appendix \ref{app0}, in which the last half of the simulated Monte Carlo steps are surely thermalized and the correlation time was estimated to be $2^8$ Monte Carlo steps. Eventually, for each realization of the quenched random couplings we have $\mathcal N=2^{10}-2^{12}$, depending on the size. We added this information in the text with the reference to Appendix A. Furthermore in appendix A we added the valure of the correlation time (256 Monte Carlo steps).

REFEREE
3- In Eq. (2) the notation \sum_{\kappa | FMC(\kappa)} is not immediately clear as there is no reference to Eq. (4), where the meaning of FMC( ) is presented. I suggest adding a reference to Eq. (4) below Eq. (2) to ease the reading. REPLY
We aniticipated a reference to Eq. (4), and also to (1), soon after Eq. (2). REFEREE
4- In order to facilitate the understanding of the plots, I would add a "T " beside the color maps in Figs. (1), (3), (4), (5), (10), and (11), so that it is immediately clear what different colors stand for. REPLY
We added the label T in all figures and in their captions as well.

Page 3, left column, last paragraph after Eq. (2): “where the expression of the sum over the indices k satisfying a FMC, like (1), will be soon clarified in Eq. (4).” has been added.
Page 3, left column, last paragraph after Eq. (3): “taking the Fourier transform of Eq. (3), “ has been added and “implies that a (t) ≃ a (t, ω)δ(ω − ω ). “ has been corrected.
Page 4, first column, second paragraph: we changed “whose effective temperature (a “photonic” temperature) is related to the ratio P between the pump rate and the spontaneous emission rate. “
in
“whose effective temperature (a “photonic” temperature) accounts both for the amount of energy E = εN stored into the system because of the external pumping and for the spontaneous emission rate. “
Equation 10: we took away 1/P^2 here and after Eq. (10) we added
“One can also introduce the pumping rate P [42] as the inverse of the square root of this ratio: $$P^2=\epsilon^2/T=1/T_{\em photonic}$$.
Page 5, first column, fourth paragraph: we changed “and together with the spherical constraint Eq. (11), it becomes an irrelevant additive constant. “ with 
“If we assume that the diagonal part of the pairwise couplings does not depend on the modes, together with the spherical constraint Eq. (11), this term is an irrelevant additive constant. The diagonal part of the linear contribution to the Hamiltonian physically represents the gain profile of the optical random medium (possibly becoming a ran- dom laser at high pumping). As a working hypothesis we are assuming a uniform gain profile over the whole spectrum. ”
Page 6, Equation (19): we changed the formula as reported in the reply to referee 2
Page 7, Section IV, first paragraph substituted by the following three paragraphs
“In a φ4 mean-field theory (a Landau theory) the critical exponents characterizing the universality class are β = 1/2 for the order parameter ⟨φ⟩, γ = 1 for the sus- ceptibility χ and ν = 1/2 for the correlation length. They satisfy the hyperscaling relation 2β + γ = νd, holding for all dimensions d ≤ duc, the upper critical dimension, that is duc = 4 in a φ4 model. As an instance, this is the universality class of the Random Energy Model (REM), a reference simplified model for the glass transi- tion. This is also the universality class of the mean-field 4-phasor model representing a random laser in the so- called narrow-band approximation, both in a fully connected interaction network, where the solution can be analytically computed [34] and in a uniformly randomly diluted version of the model, analyzed by means of equilibrium Monte Carlo simulations in Ref. [40]. Moving to the more realistic random laser models, where the basic ingredient for mode-locking, the frequency matching condition (1) is implemented, it is more difficult to understand whether the universality class re- mains the same. In Ref. [40] an estimate of the value of the critical exponent νeff ≡ 2β + γ ≃ 2/3 was provided for the mode-locked random laser model. This result is quite different from the value 2β + γ = 2 which charac- terizes the REM model, even if we consider its numerical finite-size scaling analysis. 
As an instance, the REM specific heat behaviour for small sizes N = 16, 20, 24, 28 is reported in Fig. 6. Details about the numerical technique used are given in App. B. Even though the simulated N are not very large, from the interpolation of the cV (T) peaks it turns out that νeff = 2β + γ = 1.9 ± 0.2. Strong finite size effects are there, as one can observe from the estimate of the α exponent, displaying a value α = 0.52 ± 0.07, rather different from the mean-field exponent α = 0. Be- cause of preasymptotic effects, indeed, the scaling relation 2β + γ + α = 2 (independent from the system dimension) appears to be violated. “
Former Eq. (20) has been erased.
In IV.A, page 7, second column, in Eq. (21) (former 22) $\propro \chi…$ has been added.
Page 8, first colunm, after Eq. (29) replacing “This result reproduces the value 1/ν = 1/2 for the scaling exponent of the REM. “ the following paragraph and formula have been added
“Let us recognize that Eq. (27) is the susceptibility, cf. Eq. (21), whereas Eq. (28) is the scaling of the square of the order paramater ⟨φ⟩ = φ∗, which is scaling as φ^* ~|\tau|^β. Therefore 
 |τ| ∼ 1/N1/(2β+γ)
In the mean-field φ4 theory 2β + γ = 2. Since the upper critical dimension is duc = 4 this corresponds to νduc ≡ νeff = 2, i.e., ν = 1/2 for the mean-field critical correlation length exponent. “

Page 8, second column, in Eq. (30), former (31), 2n has been substituted by n, an consistently so in Eqs. (31,32,33). After Eq. (32) the expression 
β = 1/(n-2)
has been added.
Former Eqs. (34, 35) have been modified and substituted by current Eqs. (32,33)
Page 8. The paragraph “This result implies that, in order to be compatible with mean-field theory, the values of 1/ν must fall in an in- terval defined by taking n = 2 and n → ∞ in the previ- ous expression. Eventually, the critical exponent for the scaling of the specific heat width in a generic mean-field theory must take value in the interval 
1/2 ≤ ν < 1. 
The previous argument is exact only in the large-N limit where the saddle-point approximation holds. “
has been changed into
“This result implies that, in order to be compatible with mean-fieldtheory,thevaluesofνeff =2β+γ=n/(n−2) must fall in an interval defined by taking n = 4 and n → ∞ in the previous expression. Eventually, the critical exponent for the scaling of the specific heat width in a generic mean-field theory must take value in the interval 
1≤νeff ≤2.
Given the specific theory φ and its upper critical dimen- sion duc(n), the critical mean-field exponent ν is equal to ν=νeff/duc(n). In the model under consideration, though, we have a dense (though not fully connected) interaction network and we do not have a reference d-dimensional lattice un- derneath, such that a scaling relation of the number of modes to a characteristic length can be set, as, for instance N = Ld in a d-dimensional hypercubic lattice. Our analysis will, therefore, be limited to the estimate of the exponent α, β and γ. It is also worth noting that the previous argument is exact only in the large-N limit, where the saddle-point approximation holds. “

In Eqs. (40,41), in the figures and throughout the text when needed ν is replaced by ν_eff, or, equivalently, by 2 β+ γ.
Page 9, second column, last sentence of Section IV “we can assess the mean-field nature of the glass transition in the ML 4-phasor model. “ is replaced by “we observe that the glass transition of the ML 4-phasor model is compatible with a mean-field transition. ”
Page 10, second column, Section V. We added the sentence “The number of configurations N actually used from our data can be evinced from tables II-III in appendix A, in which the last half of the simulated Monte Carlo steps are surely thermalized and the correlation time was estimated to be 28 Monte Carlo steps. Eventually, for each realization of the quenched random couplings we have N = 210 − 212, depending on the size. “ in the third paragraph.

Page 12, second column, last sentence of the first paragraph has been modified in “In this case, the same analysis performed with periodic boundary conditions on the frequencies, with similar sizes and statistics of disordered samples leads to an estimate of νeff = 2β + γ = 1.2(2), that is compatible with a mean-field theory according to the condition (34). “
Page 12, second column, last paragraph: we modified “We find clear evidence for the onset of Replica Symmetry Breaking phase at the laser threshold studying the overlap distributions. “ in “We find clear evidence for the occurrence of a Replica Symmetry Breaking phase at low temperature. Studying the deviation of the overlap distributions from a Gaus- sian distribution by standard methods (e.g., the Binder cumulant), as performed in Ref. [40], the onset of such a spin-glass phase can be shown to occur at a tempera- ture consistent with the laser threshold identified by FSS analysis of the specific heat peaks. “
Page 14, appendix A, we added the last sentence “The statistical error on the average over disorder is much larger than the error on the thermal average and leads to the error bars on the observables displayed in the main text. “
New appendix B added.
Previous appendix B moved to appendix C.
Last equation of appendix C: absolute value of parameter C replaces C.