Simultaneous evaluation of intermittency effects, replica symmetry breaking and modes dynamics correlations in a Nd:YAG random laser

Random lasers (RLs) are remarkable experimental platforms to advance the understanding of complex systems phenomena, such as the replica-symmetry-breaking (RSB) spin glass phase, dynamics modes correlations, and turbulence. Here we study these three phenomena jointly in a Nd:YAG based RL synthesized for the first time using a spray pyrolysis method. We propose a couple of modified Pearson correlation coefficients that are simultaneously sensitive to the emergence and fading out of photonic intermittency turbulent-like effects, dynamics evolution of modes correlations, and onset of RSB behavior. Our results show how intertwined these phenomena are in RLs, and suggest that they might share some common underlying mechanisms, possibly approached in future theoretical models under a unified treatment.

Nd:YAG micron size powder RL. We employed a spray pyrolysis method to synthesize the particles of Nd:YAG (4.0% mol/mol to Y 3+ ), with average diameter of 0.6 μm, see "Methods" section for details.
A series of 5000 spectra for each excitation energy was acquired to study the intermittency effects, and for the analysis of RSB behavior and modes correlations 1000 spectra were obtained in the regimes below, close to, and well above threshold. The pump source was an optical parametric oscillator pumped by the second harmonic of a Q-switched Nd:YAG laser (see "Methods" section). Figure 1a shows the spectral evolution from an excitation pulse energy below the RL threshold (0.3 mJ) to well above threshold (1.3 mJ) (the energy threshold is 0.62 mJ), whereas Fig. 1b displays the emitted intensity (left y-axis) and FWHM (right y-axis) versus excitation energy. We notice that the emitted spectrum below threshold presents large bandwidth, with the highest intensity maximum occurring at = 1064 nm (main emission of Nd 3+ ). As the RL threshold is approached, abrupt narrowing and increasing of the slope efficiency are observed.

Results and discussion
We start by defining our notation. In order to avoid confusion, we use Greek letters to represent the spectrum labels and small Latin letters to indicate the wavelength labels. We thus denote, for example, I γ i as the intensity emitted by the Nd:YAG RL at the wavelength i in the spectrum γ , with γ integer. Moreover, I i is the intensity at the wavelength i averaged over the spectra. We also define � γ i as the relative difference (fluctuation) with www.nature.com/scientificreports/ respect to this average, � γ i = � γ i / K � γ i 2 , where � γ i = I γ i − I i and capital Latin letters represent either the spectra (e.g., K = γ ) or wavelengths (e.g., K = i ). We then propose a modified Pearson correlation coefficient as follows, where above we use the Einstein summation convention over repeated indexes. The advantage to use this modified Pearson coefficient is that it comprises both the Parisi overlap parameter q γβ , which characterizes the photonic glassy RL phase with RSB, and the Pearson correlation C ij between RL modes. Indeed, on the one hand, by setting in Eq. (1) the spectrum indexes M = γ and N = β and the wavelength index K = i , we obtain the Parisi overlap parameter 33 , given by In the photonic context, each spectrum is considered a replica emitted by the system, i.e., a copy of the RL signature generated under identical experimental conditions. The probability density function (PDF) P q , analogue to the Parisi order parameter in spin glass theory 37,38 , describes the distribution of replica overlap values q = q γβ . It signals a photonic replica-symmetric paramagnetic-like phase or a RSB glassy phase if its maximum occurs at q max = 0 (no RSB) or at values |q max |≠ 0 (RSB), respectively 33 .
On the other hand, by considering in Eq. (1) the wavelength indexes M = i and N = j and the spectrum index K = γ , we write the Pearson correlation coefficient 44,45 between intensity fluctuations at wavelengths i and j , We note that, differently from Eq. (2) in which the summations are over the wavelengths, in the Pearson coefficient C ij the sums are over the spectra γ emitted at distinct times. In other words, Eq. (3) takes into account the dynamics evolution of the correlation between fluctuations of intensity at wavelengths i and j . Consider, for instance, that i and j represent the wavelengths of two modes of the laser system. Then a null value of C ij implies that these modes behave statistically in an uncorrelated way over the time interval of the spectra series. Conversely, a positive (negative) C ij signals that the statistical fluctuations in the intensity of one mode are positively (negatively) correlated to those in the other mode, so that the spatially overlapped coupled modes share (compete for) gain along the measurement.
(1) P MN = MK NK ,  . Each horizontal line corresponds to an excitation power, from below threshold (two upper lines: P/P th = 0.18 and P/P th = 0.56 ), to around (central line: P/P th = 1.32 ), and well above threshold (two lower lines: P/P th = 2.07 and P/P th = 2.45 ). Below threshold, a replica-symmetric paramagnetic-like phase identified by a unimodal central peaked P q coincides with weakly fluctuating spectra and uncorrelated intensity fluctuations with mostly blue heatmaps of the Pearson coefficient. Above threshold, the RSB glassy behavior sets in with two side peaks in P q along with stronger intensity fluctuations and increasingly correlated modes (orange and red regions in the heatmaps). www.nature.com/scientificreports/ between intensity fluctuations at distinct wavelengths, as well as self-correlations (thin diagonal red lines), are shown in color scale in the form of heatmap plots. Each horizontal line in Fig. 2 corresponds to a value of the excitation power, from the prelasing regime well below threshold, with P/P th = 0.18 and P/P th = 0.56 (two upper lines), to the RL phase around the threshold with P/P th = 1.32 (center line), and the RL regime well above threshold, P/P th = 2.07 and P/P th = 2.45 (two lower lines).
The simultaneous measure of the Parisi and Pearson parameters in Fig. 2, together with the analysis of the spectra profiles, helps to shed light on the role of the modes correlations to the emergence of the RSB glassy phase in the Nd:YAG RL. We first note in Fig. 2b,e that the spectra below threshold present rather weak intensity fluctuations in association with mostly blue heatmap patterns in the Pearson coefficient (Fig. 2c,f), evidencing that the fluctuations are fairly uncorrelated in the prelasing state. Consistently, the unimodal PDFs P q centered around q = 0 , seen in Fig. 2a,d, indicate the presence of a photonic replica-symmetric paramagnetic-like regime.
As the excitation power is increased and the threshold is approached, we observe in the central line of Fig. 2 the concomitant onset of correlations (orange and red in the heatmap of the Pearson coefficient C ij , Fig. 2i), along with a strong bandwidth narrowing in the spectra with intensity peak near = 1064 nm and much larger intensity fluctuations, Fig. 2h. A simultaneous emergence of the RSB phenomenon takes place, as indicated in Fig. 2g by the presence of the side peaks in the bimodal P q .
Finally, as the excitation power is enhanced further (last two lines of Fig. 2), intensity fluctuations nail down (Fig. 2k,n) in the RSB regime above threshold displaying a bimodal P q (Fig. 2j,m), with even stronger correlations between modes (larger red regions in Fig. 2l,o) that share the gain in an increasingly homogenous way. This phase is termed 28 a self-averaged gain regime, with Gaussian distribution of emitted intensities. It contrasts with the Lévy-like intensity statistics observed near the threshold and the first prelasing Gaussian regime below threshold 28,29 .
Since modes overlap spatially and stochastically compete for gain, we note in Fig. 2i,l,o that this competition can favor some subsets of modes, displaying stronger Pearson cross-correlation values, while undermining others with lower correlation (less intense heatmap regions).
We now turn to the joint analysis of intermittency effects along with modes correlations and RSB in the Nd:YAG RL system.
Turbulence and RSB spin glass phenomena are among the most elusive problems in physics 51 . A task even more challenging is to study their intertwined properties in the quite rare event when they coexist simultaneously in any physical system 51 . Photonic turbulence-like behavior has been recently reported 49 in a random fiber laser.
Here we also demonstrate this phenomenon in a RL, and further analyze the onset and fading out of intermittency effects with respect to the emergence and persistency of modes correlations and RSB behavior.
In fluid turbulence 51 , the relevant statistical quantities are the velocity increments between two points in the fluid, rather than the velocities themselves. In analogy, here we analyze the intensity increments, δI τ ,γ i = I γ +τ ,i − I γ i , between spectra γ and γ + τ separated in time by �t = τ t 0 , with integers γ and τ , and t 0 = 100 ms as the time interval between two consecutive spectra emitted by the Nd:YAG RL. For a proper statistical analysis including long separation times τ ≫ 1 , we had to take into account a 5 × larger number of spectra (5000). Moreover, we conveniently define a variable with null average and unit variance, x τ ,γ i = δI τ ,γ i /var δI τ ,γ i , where var δI τ ,γ i is the variance of the long time series of intensity increments. We also choose in our analysis the wavelength i of maximum intensity (once i is set, we drop the wavelength index i hereafter).
On the one hand, when nonlinearities are not relevant, as in the replica-symmetric regime with uncorrelated intensity fluctuations below threshold, the intensity increments are statistically independent and the PDF P(x τ ) of x τ α values is a Gaussian for any separation time τ (in units of t 0 ). On the other hand, as the excitation power is increased above threshold, dynamical optical nonlinearities give rise to turbulent-like emission in the intensity fluctuations of spectra close in time ( τ ≈ 1 ), causing P(x τ ) to develop a heavy tail in the asymptotic large-x τ regime 48 . However, for large separation time scales τ ≫ 1 a Gaussian P(x τ ) still holds even above threshold, in similarity to fluid turbulence 51 , in which a heavy-tailed PDF of velocity increments sets in between nearby points, but a Gaussian PDF always arises when uncorrelated distant points are considered.
From the theoretical perspective, a hierarchical stochastic model was proposed in 49 to explain the photonic turbulent-like behavior observed in an erbium-based random fiber laser. It was noticed that the PDF of intensity increments between consecutive spectra ( τ = 1 ) in the turbulent-like regime above threshold remains Gaussian but only for short time intervals in the long series of intensities, with a variance that fluctuates slowly in time. In this regime with τ = 1 , a heavy-tailed PDF P(x τ ) arises by integrating out these Gaussians along with the associated distribution of variances. The hierarchical model generally considered multiple time scales of variance fluctuations, in a way similar to Kolmogorov's hypothesis of energy cascades and intermittency in fluid turbulence 51 . The intermittency of the stochastic flux of energy between the relevant time scales is a key ingredient to set the statistical properties of P(x τ ) . Indeed, when the intermittent behavior is mitigated, for example, by decreasing the input excitation power and optical nonlinearity degree, the heavy tail of P(x τ ) fades away and a Gaussian PDF with a single scale variance emerges for any τ in the non-turbulent regime below threshold. Figure 3 shows in the right column and insets the PDF P(x τ ) of intensity increments in full agreement with the above discussion. Indeed, we notice a Gaussian P(x τ ) for τ ≫ 1 both below and above the threshold (insets of Fig. 3c,g), as well as for τ = 1 below threshold (Fig. 3d), consistently with a non-turbulent behavior. On the other hand, for τ = 1 above threshold (Fig. 3h) the PDF P(x τ ) is well fitted by a Lévy α-stable distribution 24 , with asymptotic heavy-tail power-law decay P(x τ ) ∼ 1/x α+1 τ and stability index α = 1.83 (red solid line). For comparison, we show in black dashed line in Fig. 3h the unsuccessful fit attempt by a rapidly decaying α = 2 Gaussian PDF, P(x τ ) ∝ exp −x 2 τ /2 . These results sign for a photonic turbulent-like behavior of intensity emissions in the Nd:YAG RL above threshold. www.nature.com/scientificreports/ The intermittency effects and RSB phenomenon can be analyzed jointly by defining another modified Pearson coefficient 50 which is simultaneously sensitive to both photonic glassy and turbulent-like behaviors, In Eq. (4) the new fluctuation variables are defined as d τ ,γ i = y τ ,γ i − y τ i and y τ ,γ i = I γ +τ ,i − cI γ i , with c = 0 if τ = 0 and c = 1 if τ > 0. We note that for τ = 0 the fluctuations d τ ,γ i coincide with those of the Parisi overlap parameter, � γ i in Eq. (2), used to infer the RSB phenomenon. On the other hand, for τ > 0 d τ ,γ i is equivalent to the intensity increments δI τ ,γ i in the analysis of photonic turbulence.
The first column in Fig. 3 displays the PDF P(Q τ ) for τ = 0 , indicating, just like in the first column of Fig. 2, a replica-symmetric (RSB) phase below (above) threshold with uncorrelated (correlated) intensity fluctuations, see Fig. 3a,e. Interestingly, the effect of intermittency can be also inferred from the plots in Fig. 3 by comparing the behavior of P(Q τ ) in the RSB regime for separation times τ = 1 and τ ≫ 1 as follows.
We first note in Fig. 3g above threshold that the PDF P(Q τ ) for τ = 1100 shows a bimodal profile similar to that of the Parisi overlap distribution P q in the RSB regime. In fact, since the Parisi overlap parameter in Eq. (2) considers all separation times τ between spectra, we thus expect that the statistical weight of the replica overlaps with τ ≫ 1 dominates over the long time series of spectra. Thus, P(Q τ ) should actually appear qualitatively similar to the bimodal PDF P q when the turbulence-like intermittency effects are suppressed for large separation time scales τ ≫ 1.
A remarkably distinct picture emerges for short separation times, as in Fig. 3f for τ = 1 . In this case, intermittency turbulent-like effects occur above threshold and a quite distinct statistical behavior of P(Q τ ) arises for τ = 1 , in contrast with the results for τ = 1100 in Fig. 3g. For short time scales, the intermittency tends to increase the probability of events that are rarer at much larger scales. In other words, the most likely values of Q τ when strong intermittency effects are present essentially become the least probable ones when intermittency fades away, and vice-versa. As a consequence, differently from the results for τ = 1100 in Fig. 3g, the PDF P(Q τ ) displays in Fig. 3f a unimodal profile when τ = 1 , in agreement with the above discussion.
In summary, the above discussion can be separated into the analysis of the related photonic behaviors of the Nd:YAG RL below and above threshold. In the latter both RSB and turbulent-like phenomena coexist. The first is indicated by the double-peaked Parisi overlap parameter (Fig. 3e), while the second is signaled, e.g., by the www.nature.com/scientificreports/ long-tail behavior of P(x τ ) for τ = 1 (Fig. 3h). As for the generalized Pearson coefficient, this implies a broad maximum in P(Q τ ) near Q τ = 0 for τ = 1 but also with non-negligible values for Q τ = 0 (Fig. 3f). On the other hand, below threshold RSB cannot sustain, as seen from the single central maximum of the Parisi parameter in Fig. 3a. Correspondingly, turbulent-like behavior is absent, as shown in the Gaussian profile of P(x τ ) for τ = 1 (Fig. 3d). In this case, P(Q τ ) for τ = 1 shows a quite pronounced maximum at Q τ = 0 (Fig. 3b).
Although the theoretical approach 38 to the RSB glassy phase in RLs assumes the thermodynamical limit in which the number N of modes tends to infinity, actual photonic systems in fact present values of N much away from this limit. For example, N ≈ 200 for our system and an erbium-based random fiber laser 34 and N ≈ 300 for a random Raman laser 52 and a RL based on a Rhodamine 6G dye solution with TiO 2 particles 34 , although it has been also estimated that N ~ 10 10 modes/mm 3 in a liquid-state dye RL 53 . The use of finite-size-scaling techniques standardly applied to magnetic spin glasses 37 thus emerges as an interesting perspective also for photonic systems with RSB glassy behavior. To our knowledge, no critical exponents associated with the glassy transition in these systems have been calculated so far and an ideal photonic glassy system for this purpose should be one in which the number of modes could be experimentally controlled.
In conclusion, among several striking features, RLs have been recently established as remarkable experimental platforms to study quite intricate phenomena in the multidisciplinary field of complex systems. Indeed, some of the most challenging current problems in physics find their counterparts and have been investigated in the context of RLs. Examples include photonic turbulence-like properties, spin-glass phase with RSB, and dynamics correlations between lasing modes. In the same pace, a number of statistical tools have been employed in these investigations, such as Lévy distributions [27][28][29] , Parisi overlap parameter 24,33 , and Pearson correlation coefficient [44][45][46][47][48] .
Our joint analysis of these complex photonic phenomena through a couple of modified Pearson correlation coefficients in the Nd:YAG RL showed how intertwined these behaviors are in a RL system, and indicate that they might share common physical underlying mechanisms. Although it is out of the scope of the present work to propose a theorical model to account for such phenomena, we suggest that a unique combination of key ingredients is subjacent to the mechanisms responsible for such behaviors, as the structural disorder in the distribution of gain and scattering elements, spatial overlap of multiple coupled modes, and suitable degree of optical nonlinearity, all of them present in the Nd:YAG RL analyzed here from a single set of measurements of intensity fluctuations. We hope our findings can stimulate the quest for the joint understanding of these complex behaviors under a unified theoretical approach, which can benefit from tools and concepts employed in the separate treatment of such challenging phenomena.

Methods
Sample preparation. To synthesize the particles of Nd:YAG (4.0% mol/mol to Y 3+ ), we used 54,55 50 mL of 0.2 M boehmite sol co-doped with Nd 3+ and Y 3+ as a precursor and pyrolyzed in a Spray Pyrolysis (SP) system.
For the boehmite sol preparation, 12 ml of aqueous yttrium nitrate solution (0.5 M) and 4 ml of neodymium nitrate solution (0.1 M) were added to 20 ml of ultrapure water previously heated to 83 °C. Then, 2.43 g (0.01 mol) of the aluminum tri-sec-butoxide was added and hydrolyzed for 1 h under stirring, and 0.5 ml of nitric acid was added as a peptizing agent. The final volume was adjusted to 50 mL with ultrapure water after cooling the suspension. This 0.2 M sol was spray pyrolyzed into the SP system. In this setup, the aerosol was generated at an ultrasound chamber, where a 2.4 MHz frequency piezoelectric pellet was employed. Subsequently, the aerosol was transported by a gas flow (0.1 m 3 /h atmospheric air) through two heat treatment zones: the drying zone (150 °C) and the pyrolysis one (700 °C). In the first the solvent was evaporated and the initial solid particles were attained; and in the second the final material was obtained by a fast heat treatment. At the last step, the powder was collected in an electrostatic filter (150 °C) operating at 11 kV. Finally, the collected powder was further heated at 1100 °C for 12 h.
Micron size particles and optical characterization. The Nd:YAG particles (4.0%) were characterized structurally and morphologically by X-ray diffractogram (XRD), FTIR and SEM, and the luminescent properties were studied by photoluminescence spectroscopy. The XRD showed peaks related to the structural planes of YAG, thus confirming that the YAG phase was obtained 56 . Figure 4 shows the typical morphology and particles size distribution of the powder sample used in this work. The sizes of the spherical particles ranged from 0.1 μm to 2.0 μm, with average diameter of 0.6 μm.
The photoluminescence analyses were performed on a Fluorolog 3 Horiba Scientific (model FL3-22) equipment with dual excitation and emission monochromators and Hamamatsu R928 (visible) and H10330-75 (near infrared) photomultipliers, with results shown in Fig. 5. The photoluminescence spectrum obtained under continuous-wave excitation at 808 nm exhibits narrow and well-resolved bands characteristic of the Nd 3+ f-f transitions in the YAG crystalline structure, in which the Nd 3+ ions replace the Y 3+ ones in dodecahedral sites with D 2 symmetry 57 . Fig. 6  www.nature.com/scientificreports/ second harmonic of a Q-switched Nd:YAG laser and tuned to the 810 nm emission wavelength. A repetition rate of 10 Hz was employed and the source pulse duration was 5 ns. We used a high-resolution spectrometer (SpectraPro 500, Acton Research Corporation) coupled with a charge-coupled device (CCD), covering, in real time, the range from 1040.024 nm to 1070.836 nm. In the discretization of the spectra the bin width was 0.032 nm, a value that corresponds to the resolution of the acquisition system.

Random laser setup and characterization. The experimental setup shown in
To study the RSB behavior and modes correlations, we acquired a series of 1000 spectra for each input excitation energy, and for intermittency effects 5000 spectra were obtained in the regimes below, close to, and well above the threshold.   Figure 6. Experimental setup for measurements in the Nd:YAG RL. OPO is an optical parametric oscillator, M1 and M2 are 100% reflecting mirrors @800 nm, HWP is a half-wave plate @800 nm, and PBS is a broadband polarizer beam splitter. Together they control the pump energy. L1, L2, L3 are converging lens. SH is the sample holder, NDF are neutral density filters, and the data is acquired by a high-resolution spectrometer with CCD (see text for equipment details).