Wave kinetics of random fibre lasers

Traditional wave kinetics describes the slow evolution of systems with many degrees of freedom to equilibrium via numerous weak non-linear interactions and fails for very important class of dissipative (active) optical systems with cyclic gain and losses, such as lasers with non-linear intracavity dynamics. Here we introduce a conceptually new class of cyclic wave systems, characterized by non-uniform double-scale dynamics with strong periodic changes of the energy spectrum and slow evolution from cycle to cycle to a statistically steady state. Taking a practically important example—random fibre laser—we show that a model describing such a system is close to integrable non-linear Schrödinger equation and needs a new formalism of wave kinetics, developed here. We derive a non-linear kinetic theory of the laser spectrum, generalizing the seminal linear model of Schawlow and Townes. Experimental results agree with our theory. The work has implications for describing kinetics of cyclical systems beyond photonics.


Supplementary Note 1. Basic dynamic equations
Here, we provide details of the formalism for describing wave kinetics in active cycled systems. We illustrate our method with an example of a random distributed feedback fibre laser. However, the approach is general. The random laser here is a span of optical fibre of length L that is optically pumped from fibre ends. The randomly backscattered light in the waveguide is amplified by the Raman effect, and the system starts to lase.
See more about random fibre lasers in the recent review [1].
A schematic distribution of the generated wave power along the fibre span is presented in Supplementary Figure 1. Two electromagnetic waves, -propagating to the right and the left, are generated in the fibre. Due to amplification, their amplitudes are increased during the propagation and achieve maxima near the ends of the fibre, before passing outside the fibre. Let us stress that nonlinear interaction between counter-propagating waves is weak, since their amplitude maxima are achieved at opposite ends of the fibre. Therefore, they can be considered independently.
We begin with the dynamic model describing the evolution of the envelope of the generation electromagnetic field, ψ, over the coordinate z within the fibre, at 0 < z < L, in which L is the span length. The equation for the generation wave propagating in the fibre to the right is in which t stands for time, γ is the Kerr nonlinear coefficient, and β is the quadratic dispersion coefficient. We consider the generation processes high above the threshold, so we neglect noise terms in Eq. (1). An equation similar to Eq. (1) describes the signal propagating to the left. The only difference is in the sign of the derivative ∂ z .
The effective gain operatorĝ is determined by the energy pumping and attenuation of the signal. In the frequency domain, it is a frequency-dependent function g = g R P (z) − α, in which g R is the Raman gain coefficient, P (z) is the power of the pumping wave, and α is the linear attenuation coefficient in the fibre. The distribution of energy pumping over the evolution coordinate z is defined by the factor P (z) that is assumed to be known (it can be consistently found from the power balance equations). Lasing occurs for frequencies near the maximum of the effective gain g. We can approximate spectral dependence of the effective gain by the Taylor expansion of the gain coefficient near its maximum: Further we use a notation g 0 = g R P (z) − α. Note that above the generation threshold, the condition g 0 > 0 must be fulfilled.
In reality, the power of the pumping wave P (z) depends on the generation wave ψ.
They are related via the power balance equations [2,1]. Therefore, the problem should be solved in two steps. First, we must find ψ at a given P and then self-consistently find P (z), using the balance equation. Here, we concentrate on the most challenging first step.
In a random fibre almost all generated radiation goes out from the fibre end. Only a small part of the energy is reflected back via weak Rayleigh back scattering processes.
Due to spatial distribution of the gain, amplitudes of generated waves are increased during evolution, making the scattering process feedback most effective only at the very end of the fibre. This implies an effective initial condition for the generation wave ψ + , propagating to the right, through the amplitude of the generation wave ψ − , propagating to the left.
Formally, the effective initial conditions for the waves can be introduced as in which R l and R r -reflection coefficients on the respective left end and the right end of the fibre are defined in the frequency domain. These coefficients may have different ω-dependence in different schemes. In the case of the random fibre laser, |R| 1, the reflection smallness leads to the conclusion that the signal is only weakly disturbed by the reflection, thus justifying the conditions (3).
The spectrum of the generated wave in the random fibre laser is relatively broad (in comparison to few-mode lasers) and consists of a high number of spectral components near the carrying frequency ω 0 (see [1]). The main challenge here is to describe an influence of non-linearity on the generation spectrum. To do that, we use the standard kinetic approach dealing with averaged quantities. We assume that the dispersion length (β∆ 2 ) −1 (in which ∆ is the spectral width) is small in comparison to the fibre length L. Then phases of the harmonics with different frequencies possess essentially different phases. Therefore, at averaging over a length larger than the dispersion length the harmonics can be treated as approximately independent.
The main object in kinetic theory is the pair correlation function in which angular brackets mean averaging over distances larger than the dispersion length and ψ designates complex conjugation. Due to assumed time homogeneity, the average (4) solely depends on the time difference t and is independent of t 1 . However, when examining real fibres, it is useful to average over time (integrate over t 1 ) to eliminate effects related to different fluctuations (noises) neglected in our formalism. Let us stress that due to z-dependence of the generation wave, the system is not homogeneous in space, in contrast to the time behavior. The function F is nothing else but the spectrum of the generated signal. Note that the signal intensity I can be expressed via the spectrum, it is the integral The boundary conditions (3) lead to the following relations for the averages in which F + and F − correspond to the generation waves respectively propagating to the right and the left. Further, we consider the symmetric stationary situation in which R l = R r and F + (z) = F − (L − z). Then we come to the condition for the signal propagating to the right. The condition relates values of the correlation function F taken at different ends of the fibre. Therefore, it can be considered as a closed loop.

Supplementary Note 2. Kinetics
To examine non-linear effects in the random laser, a perturbation theory must be developed. The starting point for the theory is the basic equation (1) for the envelope ψ(z, t). We treat the non-linear term in Eq. (1) as a perturbation and expand the solution of the equation over the non-linearity. Then we use the expansion for calculating the average (4). Effectively, we apply a diagrammatic technique of the type developed first by Wyld [3] in the context of hydrodynamic turbulence. The diagrammatic technique enables us to calculate correlation functions of the field ψ.
Below, we aim to derive an evolution equation for the spectrum F . The equation enables us to analyze the form of the spectrum and its dependence on the system parameters. Our derivation is performed in the spirit of the derivation of the standard kinetic equation [4,5] for classic waves. However, our system is close to an integrable one since at g = 0, the basic equation (1) is the non-linear Schrödinger equation that is completely integrable and possesses an infinite number of integrals of motion. There are no kinetics in the system of waves described by the non-linear Schrödinger equation [6]. Therefore, non-trivial kinetics is related to the presence of the gain term g that makes the situation qualitatively different from the standard kinetic equation and requires a consistent derivation of a generalized kinetic equation.
A formal solution of the equation (1) can be written as in which z is an arbitrary point and is the Green function determining a linear response of the system to an external influence.
Analogously, it is possible to consider the "back" in z evolution. For example, in the linear As a next step, we derive an equation for the function F (4). Eq. (1) yields Here, as above, the angular brackets mean averaging over distance larger than the dispersion length. The equation (11) implies that both g and F are slowly varying at the length of averaging. We assume β∆ 2 g 0 . We can then choose the averaging length l much smaller than g −1 0 . In addition, the inequality aP ∆ 2 g 0 must be satisfied as a manifestation of the spectrum narrowness in comparison with the characteristic frequency range of the Raman scattering. Therefore, we come to the chain of the inequalities aP ∆ 2 g 0 β∆ 2 that must be satisfied for our theoretical scheme to be valid.
The phase randomization caused by dispersion leads to approximately Gaussian statistics of the field ψ, since it appears to be a sum of a large number of independent terms.
Therefore, by calculating averages like (11), we can use the Wick theorem (the presentation of the average of some product of ψ fields via its pair correlation functions). However, application of the Wick theorem to the combination on the right-hand side of Eq. (11) gives zero. Therefore, we must take into account a weak correlation between different harmonics caused by non-linearity. Technically, one should exploit the non-linear contri-bution to ψ from Eq. (8) in which we substitute the expression (10). The goal of the substitution is to express δψ in terms of ψ(z, t). Then averages in the right-hand side of Eq. (12) are expressed in terms of the function F (z, ω) (4).
Using the Wick theorem, substituting the expression (9) and taking integrals over time, we obtain from Eq. (12) Here, F = F (z, ω), F 1 = F (z, ω 1 ), so the following notations are introduced The equation (13) is the generalized kinetic equation derived for interacting waves in an unstable medium (due to pumping). In Eq. (13), the usual δ-functions, which ensure the wave vector conservation in the collision integral (right-hand side), are substituted by Lorentzians, in which the gain g is presented. It is a manifestation of the system nonhomogeneity in z that is caused by the gain. Other properties of the generalized kinetic equation are close to one of the usual wave kinetic equations. For example, an integral over ω of the collision integral is equal to zero. This is a consequence of the conservation law of wave action (number of waves), which is valid without gain.
Note that similar kinetic equation with Lorentzians was found in the paper [7] where finite wave damping was taken into account. However, there is a principal and fundamental difference. In our case one cannot substitute the Lorentzians by δ-functions (as it was made in the cited work) even for the case of narrow Lorentzians since then the collision integral becomes zero due to complete integrability of the 1D Schrodinger equation. This makes our result special and very different.
It is instructive to compare our generalized equation with usual kinetic equation for weak wave turbulence [5]. The latter has two types of solutions: Equilibrium solution and flux solution, both with power spectra. In our case, the collision integral is non-zero, since it should be balanced by some additional term appearing due to the space inhomogeneity of the system. That leads to a z-dependent characteristic spectrum width. Formally, it is a consequence of the "locality" property of our collision integral (it is determined by frequencies of the order of the external frequency). The "locality" property is also characteristic of the collision integral in usual weak wave turbulence. Our results are not limited to optics, but are relevant to many other areas of science in which the non-linear Schrödinger equation is a master model and four-wave mixing of a large number of waves occurs, such as in ocean waves science [14].
It is possible to substitute g → g 0 in the collision integral since we assumed aP ∆ 2 g 0 and explained relative shortness of the non-linear stage (where the collision integral is relevant). However, generally, we should generally keep the term aP ω 2 in the left-hand side of Eq. (13), since it is relevant at the linear stage. As a result, we find In the paper we examine the case of relatively strong dispersion (wide spectrum), in which β∆ 2 g 0 (here, ∆ is the spectrum width). The inequality β∆ 2 g 0 means that we can pass to the limit of small g in Eq. (13) or (14). However, we should be careful because of the above-mentioned cancellations. In the limit g 0 → 0, the Lorentzian in the collision integral (the right-hand side of Eq. (14)), turns to δ-function of Ω, acquiring the form of the usual collision integral [5]. However, in this limit, the collision integral turns to zero. It is a consequence of the complete integrability of the one-dimensional non-linear Schrödinger equation. An existence of an infinite number of integrals of motion in this case leads to an absence of kinetics in all orders in non-linearity [6].
Therefore, we should go beyond the zero order in g 0 (that gives δ-function) and keep the first order in g 0 . Thus we can neglect g 0 in comparison to Ω in the denominator in Eq. (14) and keep g 0 in the nominator to obtain Note the presence of the singular denominator in Eq. (15). This does not lead to any divergency because of integrability (any divergency would mean that the coefficient at the δ-function is not zero). The equation is a starting point of subsequent calculations.
As it follows from Eq. (15), in the linear approximation The expression describes the exponential growth of the signal amplitude. Besides, the relation (16) shows that in the linear regime, the laser spectrum becomes narrower following the gain spectral shape g(ω). If factor A > ∆ −2 0 , in which A = dz aP and ∆ 0 are the initial spectrum width at z = 0. The spectrum width ∆ at the end of the linear stage can be estimated as ∆ ∼ A −1/2 . Note that the spectral width in this case does not depend on the initial spectral width at z = 0.

Supplementary Note 3. Diagrammatic expansion
The generalized kinetic equation (13) is derived in the lowest order of the perturbation theory. To examine higher-order corrections, we need to use more powerful tool that is the Wyld diagrammatic technique. Here, we outline this technique, as applied to the considered problem. For illustration purposes, we will re-derive the generalized kinetic equation (13) using the diagrammatic technique.
Correlation functions of the field ψ governed by the equation (1) can be calculated using the diagrammatic technique of the type developed first by Wyld [3] in the context of hydrodynamic turbulence. In the technique, the correlation functions are presented as a series over the non-linearity in the dynamic equation for the field (over γ). We use the version of the formalism developed in works [8,9,10,11] in which the Wyld diagrammatic technique is reduced to a quantum field theory.
The correlation functions are then calculated as functional integrals over the fields ψ, ψ and auxiliary fields p and p with the weight exp(−A), in which A is the "action" (Without losing generality, this section uses β = 1.) Note that an integration over the auxiliary fields guarantees validity of the equation (1) and the complex conjugated equation.
The perturbation expansion is introduced for the pair correlation function and for the Green functions Due to space inhomogeneity, the pair correlation function (18) depends on both coordinates (the correlation function introduced in the previous section is taken at coinciding points, z = z ). The functions (19,20) determine a linear response of the system to an external influence. Therefore they should possess corresponding causality properties. We formulate dynamics in terms of z, therefore G(z, z , t) = 0 =Ḡ(z, z , t) if z < z . Note that the correlation functions of the type pp are zero.
In the framework of the perturbation theory, the functions F and G,Ḡ are represented as infinite series over γ, in which different terms are determined by bare functions (obtained at γ = 0). The bare values of the Green functions, G 0 andḠ 0 , are written as After substitution g = g 0 − aP ω 2 , the integrals (21,22) can be taken explicitly. Note that G,Ḡ grow exponentially as z − z increases. The bare value of the function F is determined by the boundary condition (3).
It is convenient to introduce a graphical representation. The Green functions are depicted by double lines directed to complex conjugated fields The bare Green functions G 0 ,Ḡ 0 are depicted by single lines The pair correlation function (18) is depicted by a wavy line To construct the perturbation expansion for D, G,Ḡ, we must expand the exponent exp(−A) in the respective functional integrals in the series over the interaction term and then express the terms of the expansion via the bare objects (which can be done explicitly). The terms of the expansion can be depicted as Feynman diagrams with vertices of the fourth order (corresponding to the factor γ).
A partial summation of the diagrams leads to the following formally exact relations gq=mwΠ ω mx w, which is characteristic of the Wyld technique. Here, the blocks Σ (self-energy function) and Π ω (polarization operator) are subjects of calculation. The quantities in the first non-trivial order are given by two-loop diagrams (that can be obtained in the second The equation (27) can be analytically rewritten as A similar relation is valid forḠ. Next, applying to the equation (28) the operator G −1 0 we find dy Σ(z, y, ω)F (y, z , ω) +Ḡ(y, z , ω)Π(z, y, ω) . (33) An analogous equation can be obtained for the second argument, z . It follows from Eq.
(33) that if the z-interval is small enough and if the non-linearity is weak then Summing the equation (33), and a similar one for z , we obtain in which c.c. means a complex conjugated quantity. The next step is substituting into Eq.
(35) expressions corresponding to two-loop diagrams (29) and (30) as Σ and Π. We then use the expressions (21,22) for the Green functions and reduce two-point pair correlation functions F to single-point ones by using Eq. (34). We can then calculate the integral over y in Eq. (35) to obtain the equation (13).

Supplementary Note 4. Solution
The right-hand side of Eq. (15) can be estimated as g 0 F (γI/β∆ 2 ) 2 . We first analyze the case in which γI β∆ 2 at the end of the fibre. That means that the inequality is satisfied everywhere because I monotonically grows as z increases. The inequality Then in accordance with Eq. (15), the non-linear correction to F can be written as in which all functions, F, F 1 , . . . , are taken at z = L.
To achieve a statistically steady state, we must satisfy the relation (7). We assume here that the signal scattering is produced by impurities. The reflection coefficient R then weakly depends on the frequency ω, because the impurity size is much smaller than the wavelength. In this case the only parameter coming into the game is κ = |R ω | 2 1. It then follows from Eq. (7) To satisfy Eq. (37), we must assume that the aω 2 -contribution to the law (16) is small. Therefore: in which η

1.
Using the relations (16,36,38), we find from the condition (37) in which all functions are taken at z = L and A = L 0 dz aP . As follows from Eq. (39), the spectrum width is determined by the balance of the terms in the left-hand side: Note the smallness of ∆ in η. Comparing different terms in Eq. (39), we find I ∼ βη 3/2 /(γA) and Therefore, ∆ ∝ I 1/3 in this regime.
The equation (39) admits a self-similar substitution in which The numerical solution of the equation gives the normalization factor dx φ(x) ≈ 23.8.
We see from Eq. (41) that the spectral width ∆ grows as the intensity I increases.
At some level of pumping, γI becomes of the order of β∆ 2 . At higher pumping levels the lasing regime completely changes. This regime needs a separate consideration. Our preliminary analysis shows that in this regime, the relation γI ∼ β∆ 2 is satisfied during the non-linear stage of the generation wave propagation (near the fibre span end). These results will be presented elsewhere.
Using equation (43) the generation spectrum shape can be numerically calculated and different ratios between gain and dispersion, λ = 2g/β∆ 2 , see Fig.2 of the main text.
Note that in the limiting case of small dispersion, the generation spectrum of random laser approaches the hyperbolic secant shape. Note that in the case of long Raman fibre laser with a conventional linear cavity, similar hyperbolic secant spectrum shape was observed and introduced in a phenomenological way in the limiting case of the dispersion much bigger than nonlinearity [2]. The fact that the hyperbolic secant shape of the wave spectrum being derived with the help of classical wave kinetics in [2] conventional fiber laser coincides with the spectrum shape derived within the new approach of active cycled wave kinetics for a very different system (random lasers with no cavity) is intriguing and should be further investigated.

Supplementary Note 5. Experimental verification
To experimentally verify the predictions of the developed nonlinear kinetic theory, we designed a random fibre laser (Supplementary Figure 2).
The random fibre laser is built using 850 meters of a phosphosilicate fibre [12]. We choose the phosphosilicate fibre because of specific Raman gain profile. There is a single Raman gain peak with a spectral shape close to be Gaussian (i.e. parabolic in logarithmic scale). The frequency shift between generation wave and pump wave is about 1,330 cm −1 .
The full width at half the maximum of the gain profile is around 8 nm. Under pumping at 1,115 nm, the laser generates at 1,308 nm. We use a random fibre laser configuration with a broadband mirror of a reflectivity close to 100% from one fibre end and only random feedback from other fibre end. The pump wave is coupled from the free fibre end. We use following fiber parameters (in accordance with notations in eq.(1) and (2) Thus, the main contribution to the backscattered wave also occurs only in the same region near fibre ends. As a result, the theoretical assumption that the feedback in the system is provided by an effective point-based reflector, I − (L) = R eff I + (L), with some random reflection coefficient R(ω), is well-justified under our experimental conditions. Here, I + and I − are intensities of forward-and backward-propagating waves (see Supplementary   Figure 1). The laser is designed in the way that the generation power satisfies the inequal-ity (γI out ) 2 / (g R P in ) 2 + (0.5β 2 ∆ 2 ) 2 < 1 in a sufficiently broad range of pump powers above the threshold, allowing use of the perturbation theory for kinetic equation.
Supplementary Note 6. Linear kinetic theory near the generation threshold: Schawlow-Townes limit The wave kinetic approach developed here for active cyclic systems provides an interesting possibility to expand the seminal result by Schawlow-Townes. It also describes how, after the initial linear spectral narrowing, the wave spectrum of random fibre laser becomes broader due to the nonlinear interactions, when power increases well above the generation threshold. It is common knowledge laser physics that the lasing spectrum becomes narrower while the pump power increases inversely proportional to the power [13]. This is true when there are no nonlinear interactions [13]. In our experiments, both spectrum narrowing at low powers and nonlinear spectrum broadening at high power are observed (see Fig. 3b in the main text).
At very low pump power g R P L 1, the emission spectrum shape repeats the spectral shape of the gain profile, which we define in this section as a generalized form of (2) g(ω) = g R P e −ω 2 /Γ 2 R (44) with a width Γ 2 R = g R /a, in which a is defined in (2). Note that we don't include linear losses α into the definition of g(ω) (44) for the sake of further simplicity.
with a zero initial condition. Here, ε is a Rayleigh scattering coefficient, ε = 0.001 · α, [1], 4g ω 0 is the term considering energy input from spontaneous emission with Planck constant and carrier angular frequency is ω 0 . We remind that frequency ω is frequency detuning from the center of the gain profile.
Note that the lasing threshold defined by random distributed feedback is clearly marked on this graph at generation power around 0.01 W by abrupt changes in the spectrum narrowing law. As a result, the spectrum becomes narrower inversely proportional to the generation power, similar to the Schawlow-Townes law [13]. Therefore, we extend Schawlow and Townes' approach to describe spectrum narrowing in random fibre lasers near the generation threshold.