Statistics of extreme events in integrable turbulence

We use the spectral kinetic theory of soliton gas to investigate the likelihood of extreme events in integrable turbulence described by the one-dimensional focusing nonlinear Schr\"odinger equation (fNLSE). This is done by invoking a stochastic interpretation of the inverse scattering transform for fNLSE and analytically evaluating the kurtosis of the emerging random nonlinear wave field in terms of the spectral density of states of the corresponding soliton gas. We then apply the general result to two fundamental scenarios of the generation of integrable turbulence: (i) the asymptotic development of the spontaneous (noise induced) modulational instability of a plane wave, and (ii) the long-time evolution of strongly nonlinear, partially coherent waves. In both cases, involving the bound state soliton gas dynamics, the analytically obtained values of the kurtosis are in perfect agreement with those inferred from direct numerical simulations of the the fNLSE, providing the long-awaited theoretical explanation of the respective rogue wave statistics. Additionally, the evolution of a particular non-bound state gas is considered providing important insights related to the validity of the so-called virial theorem.

Integrable evolution equations typically arise as leading order approximations of nonlinear dispersive wave systems and often provide a very good description of the core nonlinear dynamics in a variety of physical contexts ranging from water waves to quantum gases [17].The integrable nature of the equations enables analytical solutions via the inverse scattering transform (IST) method with both zero and non-zero boundary conditions at infinity [17][18][19].Often seen as a nonlinear generalization of the Fourier method, the IST method consists of three main steps: (i) the direct spectral transform which decomposes the scattering data of the wave field at t = 0 into the so-called solitonic (related to discrete spectrum) and radiative (related to continuous spectrum) components; (ii) the (simple) time evolution of the scattering data for both components; (iii) the inverse scattering procedure for reconstructing the nonlinear wave-field at t > 0 from the evolved combined spectral data.
While the above classical, deterministic IST framework has been remarkably successful in many problems of nonlinear physics, its "stochastic" counterpart addressing the evolution of random initial data in integrable systems (essentially the theory of IT) is still in its infancy, * thibault.congy@northumbria.ac.uk with the majority of physically significant theoretical results being reliant on heavy numerical simulations (see e.g.[2,3,7]) or short-time expansions [6].The challenge is, given the statistics (the probability density function, the correlations, etc.) of the initial random data for an integrable equation to determine the statistics of the solution at t > 0.
We propose a general theoretical approach to the analysis of the long-time integrable evolution of random wave fields whose typical realizations are dominated by the solitonic spectral component.We focus on the statistics of extreme events (a.k.a.rogue waves) in random wave fields developing from certain generic classes of stochastic initial data for the focusing nonlinear Schrödinger equation (fNLSE).This fundamental problem of nonlinear physics has been the subject of extensive experimental and numerical investigations for several decades in various physical contexts [2, 5, 9-12, 14, 20-28].Recently, statistical estimates for the probability of extreme events have been derived using large deviations theory [29,30].However, an exact analytical treatment of the extreme wave statistics remains an open problem.
In this Letter we consider two ubiquitous random waves settings that exhibit extreme events in the process of the fNLSE evolution: (i) the nonlinear stage of the development of the noise-induced (spontaneous) modulational instability (MI) of a plane wave [2,14], and (ii) the evolution of the so-called partially coherent waves whose amplitude is given by a slowly varying random function with a given statistics [7,[10][11][12].
In both settings the amplitude of the initial oscillations, negligible in the MI case and finite in partially coherent waves, dramatically increases with time as depicted in Figure 1.The numerical simulations in [2,7,9] showed that, remarkably, the developed IT is characterized by statistically stationary states in the long-time regime, but the properties of these states are drastically different for the two types of random input.This is a clear consequence of integrability of the wave dynamics exhibiting infinite number of conservation laws, and in sharp contrast with the properties of classical dissipative hydrodynamic or weak (wave) turbulence characterized, in the absence of damping or forcing, by the equipartition of energy and the universal Rayleigh-Jeans Fourier spectra as t → ∞, independently of the initial data [31].In particular, it has been observed that in case (i) the fourth normalized moment κ = ⟨|ψ| 4 ⟩/|⟨|ψ| 2 ⟩ 2 of the probability density function (PDF) of the random wave amplitude ψ-the kurtosis-evolves from the initial value of 1 to 2 [2,14,26], while in case (ii) it grows from 2 to 4 in the highly nonlinear regime [7].Although in both cases the kurtosis doubles as a result of the nonlinear random wave field evolution, the former case (κ = 2) corresponds to the Gaussian statistics of the asymptotic IT while the latter (κ = 4) implies enhanced probability of high amplitude events-often described as a "heavy tail" of the PDF and associated with the rogue wave formation [10,12,14,[20][21][22]26]. The value of the kurtosis has been derived for partially coherent waves in the weakly nonlinear regime (i.e. when solitons can be ignored) in the framework of wave turbulence theory [20].To the best of our knowledge, there is no theoretical description of the kurtosis doubling in the above scenarios (in the strongly nonlinear regime).Below we present an analytical description of this phenomenon using recent developments of the spectral theory of soliton gas.FIG. 1. Spatio-temporal plot of |ψ| 2 for the asymptotic development of MI (a) and partially coherent waves (b).Snapshots at t = 0 (dashed red line) and t = 100 (black solid line) for MI (c) and partially coherent waves (d); ℓ represents the typical width of the initial pulses composing the partially coherent wave.The details of the numerical implementation are given in [7,26].
Soliton gas (SG) can be seen as an infinite stochastic ensemble of interacting solitons randomly distributed on the whole line [32].It represents a prominent example of IT that has been attracting a great deal of interest recently due to the recognition of its ubiquity in various physical systems [8,16,[33][34][35].
The theory of SG was initiated in Zakharov's 1971 paper [36] by considering an infinite collection of wellseparated (weakly interacting) KdV solitons randomly distributed in space and having some given distribution over the IST spectral parameter {λ k }-the discrete spectrum.This theory of rarefied SG has been significantly expanded by considering dense (strongly interacting) KdV and fNLSE soliton gases within the mathematical framework of the thermodynamic limit for spectral finite-gap solutions and their modulations [37][38][39].The key observation is that, in both MI and partially coherent wave settings the dynamics are dominated by the solitonic component of the spectrum so that the long-time behavior of the IT can be approximated by appropriate SGs.In what follows we shall take advantage of the fNLSE SG theory [39] to infer important statistical characteristics of the developed, homogeneous IT.This will be done within the "stochastic" version of the IST schematically shown in Figure 2. Specifically, one extracts the spectral statistical distribution-the socalled density of states (DOS)-of the approximating SG from the direct scattering analysis of random initial data ψ(x, t = 0) (solid line in Figure 2), and then reconstructs the statistics of the long-time asymptotic IT wave field ψ(x, t → ∞) via the inverse transform of the SG DOS (dashed line) using the invariance in time of the global spectrum statistics for problems with macroscopically homogeneous random initial data.The time evolution step for scattering data of the traditional IST is replaced by the assumption of an effective stochastization of the soliton phases ensuring spatial uniformity and statistical stationarity of the SG at t → ∞.Although generally, the direct scattering transform of a random potential represents a complex mathematical problem [40], we show that the determination of the spectral DOS for the two cases considered here reduces to evaluating the Abel transform of the PDF of the random initial data.We then use the relations between the spectral DOS and the fNLSE conserved densities and currents [41] to determine the kurtosis κ ∞ of the developed IT.
We consider the fNLSE in the form The discrete, soliton spectrum {λ k } in the linear scattering problem associated with (1) lies in complex plane [42].
If the spectrum contains only one point λ = ξ + iη in the upper half-plane C + , the wave field ψ(x, t) is given by the (bright) soliton solution where x 0 , σ 0 ∈ R are constant parameters describing the soliton spatial position and the phase at t = 0 respectively; 2η > 0 represents the soliton amplitude, −4ξ-the soliton velocity.The fNLSE supports N -soliton solutions characterized by the spectral set {λ k = ξ k +iη k } N k=1 complemented by the set of N complex norming constants related to the initial "positions" x k 0 and phases σ k 0 of the individual solitons within the N -soliton solution [42] (see Supplemental Material [43] for details).The N -soliton solution is called a bound state if all ξ k ≡ 0.
The fNLSE SG can be formally defined via the limit as N → ∞ of N -soliton solution with discrete spectrum points λ k distributed with some density over a domain Γ + ⊂ C + and appropriately chosen random distributions for the norming constants ensuring certain non-zero spatial density of SG on R. The key aggregated characteristics of SG is the DOS f (λ; x, t), defined as the local density in the spectral phase space Γ + × R so that f (λ; x, t)dξdηdx is the number of soliton states contained in the element [ξ, ξ + dξ] × [η, η + dη] × [x, x + dx] of the phase space at time t.For a spatially homogeneous, equilibrium SG f ≡ f (λ) globally i.e. f t = f x ≡ 0. fNLSE soliton collisions are pairwise, elastic and are accompanied by certain position and phase shifts [42].As a result, the effective velocity s(λ) of a tracer soliton in a SG is different from its velocity −4ξ in a "vacuum" (free soliton velocity) and is defined by the SG equation of state [38,39] describes the asymptotic spatial shift in a two-soliton collision [42].In this Letter, integrations are written for a 1D curve Γ + .If Γ + is a 2D domain in C + , the arc integration Γ + . . .|dλ| should be replaced by Γ + . . .dξdη [39].
The fNLSE has an infinite number of conservation laws (p j [ψ]) t + (q j [ψ]) x = 0 (j ≥ 1), see for instance [17].We focus in the following on p 3 is commonly called the energy density, where H L = |ψ x | 2 dx corresponds to the linear kinetic energy of the physical system, and H NL = − |ψ| 4 dx to the nonlinear interaction energy.It was shown in [41] that ensemble averages of the densities ⟨p j ⟩ and currents ⟨q j ⟩ in SGs are given by the moments of the DOS such that where the spectral average h(λ) = Γ + h(λ)f (λ)|dλ| (see Supplemental Material [43] for an alternative derivation of relations (S23)).Relations (S23) are based on ergodicity of a homogeneous SG inherent in the finite-gap construction of [39], so ⟨. ..⟩ in (S23) can be seen as spatial or temporal averages over an infinite period.A linear combination of these averages yields the value of the kurtosis for a uniform SG in terms of the spectral DOS.
The special case of the bound state SG is described by the DOS f (λ) = f (η)δ(ξ) where δ(ξ) is the Dirac delta function.Since the corresponding free soliton velocity vanishes, the equation of state (3) has the solution s(λ) = 0, and the kurtosis expression (6) simplifies to where η k are moments of the reduced DOS f (η).
We now consider the application of the general result (6) to the two fundamental scenarios of the IT development.We first consider the problem of spontaneous MI of the plane wave solution ψ(x, t) = e 2it of (1).At t = 0 the plane wave with small random perturbation, a "noise" ϕ(x): with |ϕ| 2 ≪ 1 and zero average, ⟨ϕ⟩ = 0, where ⟨. ..⟩ stands for the ensemble average.The plane wave solution is unstable with respect to long-wave perturbations which grow exponentially with time for t ≪ 1 [46,47].The solution ψ(x, t) develops into an incoherent, strongly oscillating structure, and the wave field statistics becomes stationary in the long time regime [2].It was shown in [26] that the asymptotic dynamics of the spontaneous MI can be accurately modeled by the uniform, bound state SG with the DOS The distribution (9), sometimes called the Weyl DOS, corresponds to the "soliton condensate" associated with the spectral support Γ + = [0, i] [39].Substituting the Weyl DOS in (7), we obtain κ ≡ κ ∞ = 2 which is twice the initial kurtosis, in perfect agreement with the longtime limit computed numerically in [26].This result provides the long-awaited analytical proof of the striking statistical property of the spontaneous MI originally discovered in the numerical simulations of [2]: the long-time dynamics exhibits large amplitude fluctuations characterized by Gaussian single-point statistics-a counterintuitive result for strongly nonlinear IT.
We now consider a more general random initial condition-a "stochastically modulated" plane wavewith the following slow variation assumption for the amplitude and phase: where , see e.g.[6].Additionally, the random process ψ(x, 0) is assumed to be ergodic; an example of such initial condition is displayed in Figure 1d.The field ψ(x, 0) satisfying ( 10) is called a partially coherent wave and can be regarded as an infinite sequence of broad "humps" with random distributions for the width O(ℓ) (ℓ ≫ 1), the amplitude O(1), and the position.The randomness of such a wave is realized on the macroscopic scale L ≫ ℓ, while each slowly varying hump corresponds to a coherent structure (details of the numerical implementation can be found in [6,48,49]).At an early evolution time, each single hump exhibits a smooth evolution dominated by nonlinearity [6] which culminates in the emergence of a gradient catastrophe followed by the dispersive resolution via an ordered sequence of coherent structures locally well approximated by the Peregrine breather solutions of fNLSE [50,51].Eventually, at t → ∞, the solution decomposes into a statistically uniform SG as described below.
In an idealized partially coherent wave with u 0 (x) ≡ 0 each hump can be approximated at leading order by a non-propagating, bound state N -soliton solution in the semi-classical limit (N → ∞) [52,53].Within a physically realistic partially coherent wave satisfying (10) each soliton has a small but non-vanishing velocity component ξ k ̸ = 0 (see [54,55] for precise analytical estimates), as depicted by the trajectories of solitons in Figure 1b.However the real part of the soliton spectrum only contributes to a small correction to the averages (S23) and can be neglected in the computation of the wave field statistics.
Since the stochastic process ψ(x, 0) is ergodic, the statistics of the IST spectrum of partially coherent waves can be determined from one representative realization of ψ(x, 0).In the semi-classical setting ℓ ≫ 1, the spectral distribution of the initial condition for x ∈ [0, L], is approximated by the Bohr-Sommerfeld density φ [42].Since the fNLSE evolution is isospectral, the global spectrum statistics is invariant in time and the DOS of the homogeneous SG at t → ∞ is given simply by f (η) ∼ φ L (η)/L as L → ∞.Using the change of variable x → ρ = ρ 0 (x) and the standard geometrical definition of the PDF (see Supplemental Material [43] and [56]), we obtain that the DOS is given by the Abel transform of the PDF of the field ρ 0 (x), denoted P 0 (ρ): A similar result was derived for the KdV SG in [57].
Clearly the Weyl DOS (9) for the initial data in the MI scenario is obtained by taking P 0 (ρ) = δ(ρ − 1), corresponding to the plane wave solution ψ = 1 at t = 0. We can now replace the average over randomly distributed partially coherent waves by the average over different realizations of the SG described by (11).As an illustration, we generate numerically partially coherent waves with Gaussian single-point statistics implying the exponential PDF P 0 (ρ) = exp(−ρ) [31].Using formula (11) we obtain the Rayleigh distribution f (η) ∼ η exp(−η 2 )/ √ π which yields by ( 7) κ ∞ = 4 in the long-time regime, which is twice the kurtosis at t = 0. Our theory thus explains the largest value of the asymptotic value κ ∞ = 4 observed in the numerical simulations performed in the large nonlinearity regime in [7].Similar to the spontaneous MI scenario, we can infer that the probability of high amplitude waves drastically increases with time for partially coherent waves.
The doubling of the initial kurtosis is a general feature of partially coherent waves in the semi-classical limit, i.e. when the solitons velocity can be neglected to leading order.Indeed (11) yields a relation between the moments of the SG DOS f (η) and the moments of the initial PDF P 0 (ρ), in particular, |ψ(x, 0)| 2 = 4η and |ψ(x, 0)| 4 = 16 3 η 3 (see Supplementary Material [43]).Thus, the formula (7) derived for the bound state SG implies that, regardless of the expression for the initial PDF, the kurtosis of the IT developing as t → ∞ satisfies An alternative derivation of ( 12) can be found in the Supplemental Material [43].
The general kurtosis formula ( 6) is also valid for nonbound state SGs (ξ ̸ = 0).We consider the so-called circular soliton condensate defined in [39] by the DOS supported on a semi-circle in the complex plane: The substitution of the DOS (13) in the equation of state (3) yields the effective velocity s(λ) = −8ξ, which is twice the free soliton velocity, far from the bound state regime.Now Eqs. ( 6), ( 13) yields κ = 2, similar to the modulational instability induced IT, which compares very well with the value computed numerically for circular condensates (see Supplemental Material [43]).Although κ = 2 for both the bound state SG generated by MI and the circular soliton condensate, the energy averages ⟨H L ⟩ and ⟨H NL ⟩ are drastically different for the two SGs.The average current ⟨q 2 ⟩ vanishes for bound state SGs (see formula (S23)) implying the relation i.e. the average interaction energy is twice the average kinetic energy.A dynamical analogue of ( 14) with ⟨. ..⟩ corresponding to a spatial integration and known as the virial theorem, has been previously derived for 2D and 3D fNLSEs using spatial zero boundary conditions [58,59].In the bound state SG context, one can assume zero boundary conditions for any sufficiently large spatial interval due to the cancellation of the solitons velocity.
Summarizing, we have formulated a general theoretical framework for the IST analysis of IT and have shown that statistical moments of the long-time development of IT can be effectively computed for certain classes of random initial conditions using the SG approximation.
In particular, we have analytically explained the asymptotic doubling of the kurtosis for two ubiquitous nonlinear wave phenomena: the long-time evolution of spontaneous MI [2] and partially coherent waves [7].Concluding, our work paves the way to the determination of the full statistics (i.e. the PDF, the correlations, etc.) in IT and, ultimately, to the realization of the stochastic IST schematically shown in Figure 2.
Appendix for: Statistics of extreme events in integrable turbulence

I. N -SOLITON SOLUTION OF FNLSE
The numerical modeling of the fNLSE SG is achieved via the construction of N -soliton solutions with N large and appropriately chosen distributions for the discrete spectrum and norming constants.In our numerics we used the efficient approach developed in [44].
The N -soliton solution of the fNLSE is given by the ratio of two determinants (see e.g.[17,18]) where . ., λ N ) are the solitons spectral parameters, x k 0 ∈ R their initial "positions" and σ k 0 ∈ [0, 2π) their initial phases.Note that x k 0 coincides to the initial position of the soliton with spectral parameter λ k only when the solitons are well separated (e.g.realization of a rarefied SG).It is not the case for realizations of the circular soliton condensate presented thereafter.The norming constants at t = 0 are If the solution (S2) is not a bound state (Re(λ j ) ̸ = Re(λ k ) for j ̸ = k), it asymptotically reduces to a superposition of well-separated solitons in the limit |t| → ∞ [42] where ψ s is the one-soliton solution: with λ = ξ + iη, and where x k ± , σ k ± are the positions and phases of the k-th soliton at t → ±∞.Because of the interaction between the N solitons, occurring at finite time t, the position x k − and phase σ k − are different from the position x k + and phase σ k + .In numerical applications (e.g.realizations of the circular condensate), ψ N (x, t; λ) can be evaluated efficiently using the dressing method presented in [44].The numerical scheme is subject to roundoff errors during summation of exponentially small and large values for a large number of solitons N , and the implementation of high precision arithmetic routine is necessary to generate solutions with the number of solitons N > 10.We note that the fNLSE considered in [44] has a different normalization than (S1), and one should substitute t by 2t in the dressing method to obtain the N -soliton solution of (S1).

II. NUMERICAL IMPLEMENTATION OF RANDOM INITIAL CONDITIONS
The evolution of random initial conditions, noise (MI) and partially coherent wave, displayed in Figure 1 in the main text is obtained via direct numerical resolution of (S1) with periodic boundary conditions: ψ(−L/2, t) = ψ(L/2, t).The corresponding initial value problem is solved via a pseudo-spectral, adaptive fourth order Runge-Kutta method with the spatial discretization ∆x = L/N x , N x = 2 12 .
Both the noise and the partially coherent wave is implemented with a sum of incoherent, discrete Fourier components: with k j = 2π L j and σ j 0 uniformly distributed between 0 and 2π.The initial conditions implemented to generate Figure 1 in the main text are: • Noise: ψ(x, 0) = 1 + Ψ RP (x) with L = 200, N modes = 4096, ∆k = 11.3, and Ψ 0 = 2.4 × 10 −4 .
Both initial conditions develop into realizations of spatially uniform SGs in the long time regime.The numerical value of the kurtosis is obtained by averaging solutions of the initial value problem (with different randomly distributed initial conditions) in the long time regime.The averaging procedure corresponds to an ensemble average, over the different realizations, as well as a spatial average between x = −L/2 and x = L/2.
The evolution of the kurtosis in time has been precisely determined numerically in [2] for MI and in [7] for partially coherent waves.In the first case, after an initial exponential growth of the noise due to MI (similar to evolution depicted in Figure 1a in the main text), the value of the kurtosis oscillates around the expected value κ ∞ = 2. Authors in [2] have observed in numerical simulations that the amplitude of the oscillations decays with time as t −3/2 and a stationary regime establishes in the long time regime.In the second case, the stationary regime established after a finite time t = O(ε −1 ) where ε is the dimensionless, small parameter characterizing the semi-classical limit.Since the typical amplitude of the partially coherent waves is O(1) in this work, the semiclassical parameter is simply given by ε ∼ ℓ −1 where ℓ is the typical width of the humps.
Since the problem implemented numerically has periodic boundary conditions, one could expect Fermi-Pasta-Ulam-Tsingou (FPUT) recurrence phenomenon for which one realization ψ(x, t) evolves back to the initial condition ψ(x, 0).In that sense, the initial value problem solved numerically does not reach an asymptotic state described by stationary statistics.The stationarity is achieved for L → ∞ when the recurrence time also goes to infinity.In practice, the numerical simulations are performed with L = O(10 2 ) such that a "quasi-stationary" regime establishes well before the FPUT recurrence; the procedure is detailed in the reference [2].Note that we do not observe FPUT recurrence numerically in the range of parameters considered in this work.

III. CIRCULAR SOLITON CONDENSATE
The DOS of the circular soliton condensate is given by formula (13) in the main text [39]: (S7) The effective velocity is s(λ) = −8Re(λ).The spectral parameter λ ∈ Γ + can be parametrized by an angle: λ = ξ + iη with ξ = cos θ and η = sin θ such that the arc integration Γ + . . .f (λ)|dλ| becomes π 0 . . .sin(θ)dθ/π.Equations ( 4) and ( 5) from the main text yield: yielding the value of the kurtosis κ = 2.The special SG termed in [39] the "circular soliton condensate" is implemented numerically using the N -soliton solution with N ≫ 1 and the spectral parameters λ k distributed in a certain way along a circle in the complex plane.Since the gas of interest has stationary statistical properties, we distribute the spectral parameters and norming constants such that ψ N (x, t; λ) approximates a typical realization of the SG in a finite region of space x ∈ [−L, L]; we set here t = 0, such that the norming constants are given by (S3).Note that the homogeneous SG is considered here as established at the outset, and does not result from the long time evolution of random initial conditions as depicted in Section II.
The discrete spectral parameters λ k of the approximating N -soliton solution (S2) are sampled from the continuous distribution on Γ + with the normalized density f (λ)/κ where corresponds to the spatial density of solitons, i.e. the number of solitons per unit interval of x ∈ R. Norming constants' norm |C k 0 | and argument −σ k 0 are uniformly distributed in the intervals [1 − α, 1 + α] and [0, 2π) respectively, with α < 1.A typical N -soliton solution with N = 100 and α = 0.15 is shown in Figure S1a.If α is too small, solitons accumulate at the position x = 0 and the SG statistics is no longer uniform; in particular the N -soliton solution is symmetric if α = 0 [26].In practice we choose a value of α such that the SG statistics is uniform in the neighborhood of x = 0.
Moments of the wave field ⟨h[ψ]⟩ are obtained numerically with a realization averaging and a local, spatial averaging of the N -soliton solutions where M is the number of "realizations".The increase of the kurtosis above 2 in the neighborhood of x = 0 for a smaller value of α is a numerical artifact of the implementation as explained above, the kurtosis decreasing then to theoretical value 2 for |x| > 5.The kurtosis also increases outside the region [−L, L] as the soliton density decreases, and the N -soliton solutions no longer represent realizations of the dense, circular soliton condensate.

IV. AVERAGES OF DENSITIES AND FLUXES
In this section we compute ensemble averages of the densities p j [ψ] and currents q j [ψ] for a homogeneous SG with DOS f (λ).A similar physical derivation for defocusing NLSE and Kaup-Boussinesq SGs is presented in [45].The nonlinear wave field in a homogeneous SG represents an ergodic random process, both in x and t, and the ergodicity property implies that ensemble averages in the where we assume the property sgn(s(λ)) = sgn(−4ξ), with −4ξ being the free soliton velocity, see [45].

V. PARTIALLY COHERENT WAVE A. Derivation of the DOS
As indicated in the main text, the distribution of the discrete IST spectrum of a partially coherent wave for x ∈ [0, L] is approximated by the Bohr-Sommerfeld density The change of variable ρ = ρ 0 (x) yields where dρ] makes the contribution η/π ρ − η 2 , weighted by i |dx i |, to the Riemann sum in (S37).Thus φ L (η) can be seen as a "linear superposition" of Weyl's distributions (Equation (10) in the main text).
The total number of solitons "contained" in the portion x ∈ [0, L], L ≫ ℓ, of the partially coherent initial data is N L ∼ φ L (η)dη.The "global" (i.e.defined for large spatial scales ∆x ∼ O(L)) DOS corresponding to the partially coherent wave is then given by as L → ∞.The sum i |dx i |/L corresponds to the probability P 0 (ρ)dρ to have ρ 0 (x) ∈ [ρ, ρ + dρ].Isospectrality of the fNLS evolution implies conservation of the global DOS so that after the generation of the homogeneous SG via the stochastization of soliton phases at t → ∞ the conventional, locally defined, DOS f (η) = f 0 (η) since for the homogeneous SG the local and the global DOS coincide.This yields formula (11) in the main text.We emphasize that the global DOS f 0 (η) (30) is defined for the partially coherent initial data, while the expression of the kurtosis (Equation ( 7) in the main text) is valid only for a homogeneous SG, i.e. in the long time regime.

B. Relation between the PDF moments and the DOS moments
The formula for the DOS of partially coherent wave (Equation (11) in the main text) yields a relation between the moments of the reduced DOS, η j = ∞ 0 η j f (η)dη, and the moments of the PDF P 0 (ρ) of the initial wave field.Consider η 2k+1 :

C. Alternative derivation for the kurtosis doubling
The result κ ∞ = 2κ 0 for partially coherent waves (Equation ( 12) of the main text) can be inferred directly from the definitions of the conserved quantities and fluxes (Equation (4) of the main text) and the basic assumption of the long-time resolution of a partially coherent wave via the bound state SG.
First we observe that the relation ⟨q 2 ⟩ = 0, valid for the bound state SG, implies that for the long-time evolution of a partially coherent wave we have (S44) FIG. 2. Stochastic analogue of the nonlinear Fourier (IST) framework for the fNLSE with initial data in the form of a random potential dominated by solitonic content Figure S1b displays the variation of the kurtosis computed over M = 90000 N -soliton solutions (N = 100), and a local spatial average with ℓ = 30 for α ∈ {0.07, 0.15}.If α = 0.15, the kurtosis is approximately uniform in the region x ∈ [−50, 50] with L = 50, and we can assume that N -soliton solutions describe realizations of the circular SG in the region [−L, L].This assumption is confirmed by the agreement between the kurtosis determined analytically κ = 2, and the kurtosis evaluated numerically in the region [−L, L] which is approximately equal to 2.