Coexistence of turbulence-like and glassy behaviours in a photonic system

Coexistence of physical phenomena can occur in quite unexpected ways. Here we demonstrate the first evidence in any physical system of the coexistence in the same set of measurements of two of the most challenging phenomena in complex systems: turbulence and spin glasses. We employ a quasi-one-dimensional random fibre laser, which displays all essential ingredients underlying both behaviours, namely disorder, frustration and nonlinearity, as well as turbulent energy cascades and intermittent energy flux between fluctuation scales. Our extensive experimental results are theoretically supported by a newly defined photonic Pearson correlation coefficient that unveils the role of the intermittency and describes remarkably well both the spin-glass Parisi overlap parameter and the distribution of turbulent-like intensity increments. Our findings open the way to unravel subtle connections with other complex phenomena, such as disordered nonlinear wave propagation, Lévy statistics of intensity fluctuations, and rogue waves.

SCienTiFiC REPoRtS | (2018) 8:17046 | DOI: 10.1038/s41598-018-35434-z equation, bridges the gap between hydrodynamics and disordered systems by providing the system's partition function in terms of the velocity field, the relevant quantity in fluid turbulence. Addressing this problem for disordered systems was possible under an approach that applied Parisi's RSB mechanism 23 . More recently, evidence of RSB was reported 24 in disordered nonlinear wave propagation, which admits 4 a statistical mechanics description in terms of wave turbulence. These findings suggest the existence of an interplay between turbulence and spin glass that possibly hides a number of striking intertwines at a deeper structure level.
Here we employ a RFL system as the photonic platform to demonstrate the first experimental evidence in any physical system of the coexistence of turbulence-like and spin glass behaviours from the same set of measurements. Our theoretical analysis introduces a novel photonic Pearson correlation coefficient that describes remarkably well both the spin-glass Parisi overlap parameter and the distributions of turbulent-like increments arising from a hierarchical model for the multiscale intensity fluctuations. Moreover, the Pearson coefficient also unveils the role of the intermittence to the RFL properties.

Results
The RFL features and experimental details are described in the Methods section. Figure 1 shows a pictorial illustration of the multimode RFL system pumped by a continuous-wave (CW) semiconductor laser. We remark that an exceptionally large number (~10 3 ) of random Bragg gratings were inscribed in the erbium-doped fibre core to strengthen the degree of quenched disorder 25 . The Bragg gratings play the role of quenched random scatterers. This feature makes the RFL to differ substantially from conventional (i.e., not disordered) fibre lasers, in which light performs regular round trips 3 . The nonlinearity arises from the gain provided by the erbium ions, whereas frustration occurs due to the nonlinear interactions of the modes in the strongly disordered fibre medium.
The CW character of the pump source allows for a continuous exploration of the free-energy landscape, particularly in the photonic spin-glass phase above threshold, with the rugged multi-valley structure yielding the correlated states to be sequentially generated through the jump dynamics over energy barriers. Such continuous path can be hampered in the case of pulsed pumping, in which correlations are among states resulting from distinct laser shots. In addition, the laser shots may also hamper the intermittency effect observed in the intensity increments between subsequent spectra that gives rise to a turbulent-like regime in a CW pumped RFL system 7 . Moreover, this erbium-based RFL system provided the first observation 7 of the statistical signatures of turbulent-like emission in a CW-pumped disordered random laser. Indeed, in ref. 4 the authors found turbulent emission in quasi-CW Raman fibre lasers in the absence of any form of built-in disorder.
An extensive number (N s = 150,000) of emission spectra were collected for a value of the excitation power above the random-lasing threshold (P/P th = 2.92). A very long series of output intensities {I α (k)}, α = ... N 1, , s , was thus generated for each one of the 512 wavelengths indexed by k in the interval [1489.0 nm, 1591.4 nm] around the peak emission at 1540.0 nm. We notice that α can be also related to the time variable t = αt 0 , where t 0 = 100 ms is the integration time between two consecutive spectra acquisitions.
For P/P th = 2.92 the RFL exhibits a photonic spin glass phase with RSB 11,12 and turbulent-like signature in the distribution of successive intensity increments 7 . The central novelty of the present work lies in the discovery and theoretical characterization, in remarkable agreement with the experimental data, of the coexistence of turbulence-like and spin glass behaviours in the same set of measurements in the RFL system, motivated by the long-standing question in the literature 23 .
Before performing the experimental analysis, we first describe our theoretical framework. The photonic RSB spin glass behaviour is characterized by the distribution P(q) of values of the Parisi overlap parameter 9 , where α β = ... N , 1, , s denote the replica labels, the average intensity at the wavelength indexed by k is = ∑ α α I k I k N ( ) ( )/ s , and the intensity fluctuation is Δ = − α α k I k I k ( ) ( ) ( ). Each emission spectrum defines a replica, i.e., a copy of the RFL system under fairly identical experimental conditions. In the photonic RSB spin Figure 1. The photonic random fibre laser system. Pictorial description of the 30 cm long erbium-doped RFL with a large number (~10 3 ) of specially-designed random fibre Bragg gratings inscribed in the core. The pump source was a semiconductor laser operating in the CW regime at 1480.0 nm. The multimode RFL emission peaks at 1540.0 nm. A rather extensive number of 150,000 spectra were collected at the excitation power above the RFL threshold, P/P th = 2.92. Below the threhsold, a replica-symmetric prelasing phase sets in with nonturbulent behaviour. In contrast, at P/P th = 2.92 the RFL system exhibits the coexistence of photonic RSB spin glass and turbulence-like behaviours in the distribution of intensity increments between successive spectra. glass state for P/P th > 1, P(q) displays a bimodal profile 8,9 with two maxima around q = ±1, whereas in the replica-symmetric prelasing phase for < P P / 1 th a single maximum occurs at q = 0. As replicas with any 'temporal' separation τ ≡ β − α are considered in (1), thus P(q) is not appropriate to describe photonic turbulence.
On the other hand, if intensity increments are defined for fixed τ and k as δI ατ (k) ≡ I α+τ (k) − I α (k), then the distributions P(δI τ ) can be used to infer the photonic turbulent-like state at time scale τ (here we focus on a wavelength around the peak emission), as explained in the following. In homogeneous and isotropic fluid turbulence, the distribution of velocities is Gaussian, but the distribution of velocity increments between different points in the fluid can be Gaussian or not, depending on whether their spatial separation is, respectively, large or short scale 1 . In the latter, the non-Gaussianity of the distribution of velocity increments, evidenced by the presence of heavy tails, provides the signature of the turbulent state. Similarly, in the present photonic RFL system the distribution of intensities is Gaussian 12 for P/P th = 2.92. Nevertheless, as we shall see below, the time scale τ separating the emission spectra is actually determinant to the behaviour of the distribution of intensity increments as being Gaussian (τ ≫ 1) or not (τ ≈ 1). In particular, we show in what follows that this non-Gaussian behaviour of the RFL increments distribution can be actually described remarkably well by a statistical turbulence model that accommodates relevant ingredients of the fluid turbulence phenomenon, such as intermittency and energy fluxes (see below).
In this context, we start by assuming 7 the existence of N relevant time scales of RFL intensity fluctuations, with an intermittent flux of electromagnetic energy between them. The increments distribution P(δI 1 ) between successive spectra at the shortest time scale τ = 1 can be calculated exactly from a stochastic model 7,26,27 that considers a hierarchy of these N coupled time scales (see Methods for details), similarly to Kolmogorov's turbulence mechanisms of energy cascade and intermittency 1 . Noticeably, at small intervals of the long time series {I α (k)} we observe that the intensity increments δI ατ are Gaussian distributed with variance ε τ in a local scale, thus defining a conditional Gaussian density P(δI τ |ε τ ). Since the intermittency effect leads to a slowly fluctuating local variance ε τ with distribution f(ε τ ), then P(δI 1 ) can be obtained from the superposition integral The solution of the N-scale hierarchical model is 7,26,27 where G denotes the Meijer-G function 28 , ε 0 is the variance at the largest scale, β essentially represents the ratio between the deterministic and stochastic model coefficients, We observe from equation (3) that the ratio (2ε 0 /ω) 1/2 sets a typical scale for the distribution of intensity fluctuations δI τ with τ = 1. As we shall see below, one feature of the turbulence-like behaviour is the large-δI 1 asymptotic behaviour of P(δI 1 ), which displays non-Gaussian decay for τ = 1. In contrast, for τ ≫ 1 a Gaussian P(δI τ ) emerges, consistently with the non-turbulent state observed at large time scales.
Although this approach is appropriate to describe photonic turbulence, it is, however, unsuitable to signalize the spin glass phenomenon. Therefore, in order to capture the interplay of both turbulence-like and spin glass behaviours, we first need to define the new variable y ατ (k) For τ = 0 this fluctuation coincides with that of the Parisi overlap parameter (1), i.e., δ ατ=0 (k) = Δ α (k). We thus introduce the photonic Pearson coefficient to measure the correlation between δ-fluctuations at a fixed τ: The introduction of this new Pearson correlation coefficient is actually a key contribution of our work, since it is simultaneously sensitive to both turbulence-like and spin glass behaviours. Indeed, for τ = 0 it recovers the Parisi overlap parameter of the photonic glassy phase, Q αβ,τ=0 = q αβ . On the other hand, for τ ≥ 1 the Pearson coefficient is able to infer the turbulence phenomenon. For example, the maximum turbulence-like signature (intermittency) is manifested at the shortest time scale τ = 1. In contrast, for τ ≫ 1 the short-time-scale effects and intermittency rapidly fade away and a crossover takes place to the non-turbulent (Gaussian) behaviour, as demonstrated below.
We now turn to the analysis of the experimental RFL data. In the following, we define the normalized intensity increments as Fig. 2(a) shows an excerpt of the large dataset, with 3000 (out of 150,000) intensity increments between successive spectra. The distribution P(x 1 ) at λ = 1539.8 nm around the peak emission, calculated from the whole dataset, is depicted in squares in Fig. 2(b). In this semi-log plot, the large deviation from the parabolic shape signalizes the non-Gaussian turbulence-like behaviour of P(x 1 ). The red solid line shows a remarkable agreement with the results of the hierarchical model. The best fit is obtained for N = 6 time scales of intensity fluctuations, with P(x 1 ) given by a statistical mixture of two functions, P a (x 1 ) and P b (x 1 ), given by equation (3) (2), is also given by a statistical mixture with the same parameters as P(x 1 ). Remarkably, as seen in Fig. 2(c) the effect of the small and large variances ε 1 are distinguishably incorporated by each mixture component. Statistical mixtures can be also found, e.g., in oceanic turbulence, due to the combination of active and quiescent modes in the distribution of kinetic energy dissipation rates 29 . In the present work we do not attempt to separate our experimental results according to some identified mechanism responsible for the statistical mixtures, but we infer from general grounds that such mechanism, in the presence of both nonlinearity and disorder, could arise from a subtle combination of stimulated and spontaneous turbulent emissions.
On the other hand, by considering a hierarchical model similar to that described in the Methods section for the intensity increments, the distribution P(Q 1 ) of the Pearson coefficient for τ = 1 can be also calculated via a superposition integral as equation (2), yielding , ω (Q) and β (Q) defined similarly to equation (3). In Fig. 2(d,e) the experimental data (circles) are plotted together with the best theoretical fits (red lines) for P(Q 1 ) and the associated variance distribution g(ε 1 ). A very nice agreement is observed for N (Q) = 2 time scales of the Pearson coefficient and parameters of a similar statistical mixture: p (Q) = 0.70, β = .
2 53 . In particular, we argue below that the unimodal behaviour of P(Q 1 ) is intrinsically related to the maximum intermittency effect underlying the turbulent-like state for τ = 1.
A remarkably distinct picture emerges for large separation time scales between the emission spectra, τ ≫ 1. Indeed, Fig. 3 displays the crossover to the non-turbulent behaviour for τ = 5000. The intensity increments at the wavelength λ = 1539.8 nm now assume a Gaussian distribution P(x 5000 ), Fig. 3(b), indicating that intermittency has been fully suppressed (i.e., N = 0 in the theoretical model). Concurrently, a bimodal profile of P(Q 5000 ) is observed in Fig. 3(c), showing a significant contrast with the unimodal profile of P(Q 1 ) that can be understood as follows. The intermittency effect tends to increase for short scales the probability of events that are rarer at much larger scales. As seen in Fig. 3(c), when intermittency is absent in the non-turbulent state for τ = 5000, values Q 5000 ≈ 0 are less probable than those around Q 5000 = ± 1. Consequently, when intermittency is maximum for τ = 1, values Q 1 ≈ 0 have their probability much increased, leading to the unimodal behaviour of P(Q 1 ) in the turbulent-like state. (The unimodal behaviour of P(Q 1 ) in the RSB spin glass phase of the RFL at P/P th = 2.92 should not be confused with the unimodality of P(q), which, instead, signals the prelasing phase with replica symmetry for P/P th < 1).
Finally, we also observe that the bimodal behaviour of P(Q 5000 ) resembles the spin-glass profile 11 of the distribution P(q) of the Parisi overlap parameter (1) for P/P th > 1. In fact, since the Parisi overlap parameter 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. As a consequence, P(Q τ ) for τ ≫ 1 should actually appear qualitatively similar to the distribution P(q) that characterizes the RSB spin glass state.
We lastly remark that Figs 2 and 3 are actually insensitive to the number of RFL modes collected, since the distributions were calculated for the single mode at λ = 1539.8 nm. Moreover, we have also confirmed the stability of the photonic spin-glass results regarding the number of modes around the peak emission in Eq. (1).

Discussion
The coexistence phenomenon addressed in this work relates two of the most challenging behaviours of complex systems: turbulence and spin glass. We highlight that the remarkable agreement between the RFL experimental data and the theoretical results was possible only due to the very large dataset analysed (>10 5 RFL emission spectra), combined with the suitable definition of a new photonic Pearson coefficient to simultaneously account for both phenomena.
We also remark that since both photonic regimes of spin glass and turbulence have been separately identified 7,11,12 in all excitation powers studied above the RFL threshold, therefore we believe that the coexistence phenomenon probed in the present work at P/P th = 2.92 is also robust for other values of P/P th above threshold.
We hope that our findings open the way to investigations on possible unforeseen coexistence of other related complex phenomena that share similar mechanisms, including Lévy statistics of intensities 16,30 , disordered nonlinear wave propagation 24,31 , and rogue waves 32,33 .

Methods
Experimental characterization. The RFL system is a 30 cm long polarization maintaining erbium-doped structure 25 fabricated by CorActive (peak absorption of 28 dB/m at 1530.0 nm, mode field diameter of 5.7 μm, 0.25 numerical aperture), see Fig. 1. By applying a specially-designed interferometric technique, a fibre Bragg grating with random phase was written over the whole fibre. Each fibre grating acts as a quenched scatterer and an exceptionally high number of scatterers (~10 3 ) was inscribed in the fibre to improve the disorder strength. The pump source was a home-assembled CW semiconductor laser operating at 1480.0 nm and presenting two longitudinal modes, as determined from speckle measurements. The linewidths were measured to be 0.70 nm using an optical spectrum analyzer (spectral resolution of 0.06 nm). Random laser emission occurs for input powers above the threshold value, P th = (16.30 ± 0.05) mW. We have also ensured 12 that the pump laser fluctuations do not affect the erbium-based RFL fluctuations. An instrument-limited laser linewidth of 0.2 nm was determined above threshold. Although this is a relatively narrow linewidth, the RFL was shown 11 to emit a large number of longitudinal modes (~200) from measurements using the speckle contrast technique. Although our RFL operates in a single transverse mode regime, the presence of a several longitudinal modes in the 30 cm length RFL confers its quasi-one-dimensional character 25 . A rather extensive number (N s = 150,000) of emission spectra were collected for the excitation power P/P th = 2.92. A spectrometer with a nominal 0.2 nm resolution at 1530.0 nm was employed, with a liquid-N 2 -cooled infrared CCD camera as the detector. A very long series of output intensities {I α (k)}, α = ... N 1, , s , was generated for each one of the 512 wavelengths indexed by k in the interval λ ∈ [1489.0 nm,1591.4 nm] around the peak emission at 1540.0 nm. The integration time between two consecutive spectra acquisitions is t 0 = 100 ms. As the lifetime of the active erbium ions in the RFL is of order of hundreds of μ s or less, then this integration time actually corresponds to averaging over many lasing emissions. We comment that this averaging process does not invalidate our method since it is the definition and characterization of the experimental replicas resulting from these many lasing emissions that actually matter to describe the emergence of distinct photonic regimes (see, also 9,14,15 , where similar averaging phenomena occur).
Hierarchical model. Here we review on the stochastic hierarchical model that has been successfully applied to describe turbulence-like properties in diverse systems, such as fluids 26,27 , financial markets 26,27 , and a photonic RFL 7 . Semi-log plot of the distribution P(x τ ), τ = 5000, of normalized increments, at the wavelength λ = 1539.8 nm around the peak emission, obtained from the whole set of 150,000 spectra (blue squares). The red solid line displays the nice fit to the results of the hierarchical model. At the large time scale τ = 5000, the parabolic shape indicates that the distribution P(x 5000 ) assumes a Gaussian form (N = 0 in the theoretical model), showing that intermittency has been fully suppressed. (c) Histogram with the distribution P(Q 5000 ) of values of the photonic Pearson coefficient calculated from 10,001 spectra displaying a spin-glass-like bimodal profile. The experimental histogram is also nicely fitted by the theoretical model (red solid line), with the mixture components of P(Q 5000 ) depicted by green dashed lines. Since the intermittency rapidly fades away for τ > 1, then P(Q 5000 ) appears qualitatively similar to the distribution P(q) of the Parisi overlap parameter (1) that characterizes the spin glass behaviour. In fluid turbulence, the significative statistical quantities are the velocity increments between two close points in the flow 1 . By analogy, in photonic turbulence we may consider intensity increments between successive optical spectra, δI ατ (k) ≡ I α+τ (k) − I α (k), with τ = 1 (shortest separation time scale between the spectra α = ... N 1, , s ) and k denoting the wavelength index. In the prelasing regime (P/P th < 1), in which nonlinearities are irrelevant, the intensity increments are statistically independent and the probability distribution for a given wavelength, P(δI τ ), τ = 1, is a Gaussian. In contrast, when the excitation power is increased beyond the threshold, nonlinearities give rise to a turbulent-like emission in which the Gaussian form of the intensity increments distribution remains valid only at a local level, with a slowly fluctuating variance ε τ for τ = 1. Therefore, the local conditional distribution at τ = 1 can be written for P/P th > 1 as δ ε δ ε πε | = − P I I ( ) exp[ ( ) /2 ]/ 2 1 1 1 2 1 1 . Consequently, in the statistical description of turbulence the non-Gaussian global form of P(δI 1 ) can be obtained by compounding the local Gaussian P(δI 1 |ε 1 ) with a background distribution of variance fluctuations f(ε 1 ), in the form of the superposition integral (2).
The background density f(ε 1 ) can be obtained exactly from a hierarchical stochastic model 7,26,27 that accommodates the basic concepts of Kolmogorov's turbulence theory 1 , namely energy cascade (the transfer of energy from large to small scales) and intermittency, which is caused by fluctuations in the energy transfer rates between relevant adjacent scales. Intermittency can be understood as the tendency of the distributions of the relevant variable, say the fluid velocity or intensity increments, to change substantially as one moves from a large scale, where inertial effects dominate, to a small scale, in which dissipation becomes dominant. Our hierarchical model is defined 7,26,27 by the following set of stochastic differential equations, The variable ε (i) represents the fluctuating variance parameter at the respective scale in the hierarchy. In our notation, the N-th hierarchy level is assigned to the shortest time scale τ = 1, that is, ε (i = N) ↔ε τ=1 , whereas ε 0 ≡ε (0) is related to the largest time scale. The coefficients γ i and κ i are positive and dW (i) denote independent Wiener processes. The first term in (6) describes the deterministic coupling between adjacent scales, whilst the stochastic term accounts for the non-linear couplings with the background variables at all scales, setting the ultimate source of intermittency. For convenience, we define the dimensionless parameter β γ κ = 2 / i i 2 , so that (2ε 0 /β N ) 1/2 determines a typical scale for the distribution P(δI τ ) of intensity fluctuations for τ = 1; see equation (3).