The focusing nonlinear Schrödinger (fNLS) equation is a ubiquitous model describing nonlinear wave propagation arising in a variety of physical settings, including deep water waves [1], optics [2], plasmas [3], and Bose-Einstein condensates (BECs) [4]. The fNLS equation is also a completely integrable infinite-dimensional Hamiltonian system [5][6][7][8][9][10], endowed with a deep mathematical structure, including the existence of infinite families of exact solutions with both zero and nonzero background describing the elastic interactions of N solitons. In addition, its initial value problem can in some cases be solved by the inverse scattering transform (IST) [11][12][13][14][15][16][17][18][19]. The formulation of the IST is based on the representation of the nonlinear evolution equation as the compatibility condition of two linear equations called a Lax pair. The first half of the Lax pair of the fNLS equation, namely the Zakharov-Shabat (ZS) scattering problem, is equivalent to an eigenvalue problem for a non-self-adjoint one-dimensional Dirac operator.
While the classical theory and applications of the fNLS equation are mostly concerned with the description of regular, deterministic wave structures, the inherent statistical nature of some physical wave phenomena in focusing media (e.g., rogue wave emergence) calls for the study of stochastic fNLS solutions, characterized in terms of the probability density function, correlation function, etc. Establishing a connection between the IST spectra of random fNLS solutions and their statistical properties in physical space represents a challenging problem, which has recently been formulated in the context of integrable turbulence [20,21]-the general theoretical framework for the description of a broad spectrum of stochastic wave phenomena in physical systems modeled by integrable equations. A particular type of integrable turbulence, termed soliton gas (SG), has recently attracted considerable attention [22,23] due to its appearance in many physical systems including water waves [24][25][26][27], nonlinear optics [28][29][30], and BECs [31,32]. Notably, SG dynamics has been shown to underpin some fundamental physical phenomena such as spontaneous modulational instability [33] and the rogue wave emergence [34] in focusing media.
The concept of a SG was introduced in Ref. [35] as an infinite collection of randomly distributed solitons with small spatial density and with a certain amplitude distribution. Soliton interactions, accompanied by well-defined phase shifts, result in a modification of the effective velocity of a "tracer" soliton in a gas over large propagation distances, enabling an approximate description of the emergent, large-scale hydrodynamics or kinetics of a weakly nonuniform/nonequilibrium SG. The kinetic description of a SG was generalized in Ref. [36] to a dense SG using finite-gap theory, and a general phenomenological construction of SG kinetic equation for a broad class of integrable systems was proposed in Ref. [37]. A systematic spectral theory of fNLS SGs was developed in Ref. [38], where it was also extended to the case of a SG on a nonzero background, i.e., a breather gas (BG). For recent advances on the spectral theory of soliton and breather gases and its relation to the generalized hydrodynamics of integrable FIG. 1. Density plots of the time evolution of the potential in Eq. ( 4) according to fNLS equation ( 1) with = 1/20, m = 0.7, and (a) without noise, (b) with zero-mean white Gaussian noise with standard deviation σ = 10 -2 . In both cases, the spatial window was taken to be the period of the IC. many-body classical and quantum systems, see Ref. [23]. For further rigorous analysis of SGs and BGs of the fNLS equation, including periodic SGs and BGs, see Ref. [39].
While BGs represent a natural generalization of SGs, possible mechanisms for their generation have remained largely unexplored. In this work we present an analytically tractable model of BG generation via fission, based on the semiclassical limit of the fNLS equation with initial data in the form of a periodic elliptic "dn" potential with elliptic parameter m ∈ [0, 1], augmented by a small random noise. It was shown in Ref. [40] that the ZS spectrum of such elliptic potentials can be characterized analytically. It was also shown in Ref. [39] that, in the semiclassical limits, one hump, nonnegative, periodic potentials [including the Jacobi elliptic "dn" potential] are compatible with the so-called thermodynamic spectral scaling [38] and give rise to bound-state SGs or BGs. The addition of a small noise to the initial condition (IC) facilitates the "phase mixing" of the finite-gap semiclasical dn potential, giving rise, in the long-time limit, to a spatially uniform and statistically stationary integrable turbulence, associated with a BG (cf. Fig. 1).
The fundamental property of long-time effective "thermalization," or relaxation to a statistically stationary state, in an integrable system was established numerically in several scenarios of the evolution of random waves in the fNLS equation [41][42][43]. Here we extend this result to the qualitatively new framework of the BG fission, which encompasses a broad range of scenarios of transition to stationary integrable turbulence depending on the value of m in the semiclassical elliptic potential: from the nonlinear development of the spontaneous modulational instability for m → 0 to the rarefied soliton gas fission for m → 1. When the noise is sufficiently small, the ZS spectrum remains essentially unchanged from that of the semiclassical elliptic potential. The isospectrality of the BG fission enables one to take advantage of the results of Ref. [21] and evaluate statistical measures of the BG, such as the mean intensity and kurtosis, in terms of the spectral density of states of the initial elliptic potential. In Ref. [21], it was predicted that the kurtosis doubles in the long-time fNLS fission of the so-called partially coherent waves into a SG. Here we show analytically that this result generalizes to the BG fission from semiclassical elliptic potentials, and we confirm our results by comparison with direct numerical simulations of the time evolution of the noise-augmented potential.
An important ingredient of our construction is the interpretation of a BG as a "composite SG" comprising of two spectrally distinct components: a regular SG and a soliton condensate, defined as a critically dense SG [38] and describing the modulationally unstable background in a BG. This composite SG provides a natural extension of the "pure" fNLS soliton condensate framework used in Refs. [33,44] to model the development of spontaneous (noise-induced) modulational instability of fNLS plane wave and genus one elliptic solutions.
The cubic fNLS equation is, in normalized and dimensionless form and in the semiclassical scaling,
where subscripts x and t denote partial differentiation, q(x, t ) describes a complex-valued envelope of oscillations, and the physical meaning of the variables x and t varies depending on the physical context. (For example, in nonlinear fiber optics, x is a retarded time and t is the propagation distance through the medium.) Equation ( 1) is the compatibility condition v xt = v tx of its Lax pair, namely the overdetermined linear system
is the third Pauli matrix, and
the asterisk denoting complex conjugate. The first half of Eq. ( 2) is the ZS scattering problem, and q and z are referred to respectively as the potential and the scattering parameter, or eigenvalue, since the first half of Eq. ( 2) can be rewritten as the eigenvalue problem Lv = zv, where
The spectrum of the dn potential (4) in the complex z plane with m = 0.7 for = 1/7. Right: Periodic (red) and antiperiodic (blue) eigenvalues (vertical axis) as a function of the elliptic parameter m in (0, 1) (horizontal axis) for = 1/7. For every fixed value of m, as → 0, the number of bands grows like 1/ . The yellow line separates the range o = (0, iq min ) (in white) from the range g = (iq min , iq max ) (in gray). The dashed line shows the value m = 0.7 corresponding to the spectrum in the left panel. For each fixed value of m, inside z ∈ g the band widths decrease exponentially in , giving rise to point spectrum in the limit → 0. Conversely, inside o the band widths decrease only algebraically, giving rise to a continuous band in the limit → 0.
In Ref. [40] we characterized the Lax spectrum (L) of the above ZS problem with the potential
where dn(x; m) is a Jacobi elliptic function [45,46] and m ∈ (0, 1) is the elliptic parameter. (Recall that the Lax spectrum (L) is defined as the set of values z ∈ C such that the focusing ZS scattering problem admits solutions v(x, z) that are bounded for all x ∈ R.) Note the spatial period is 2L = 2K m , where K m = K (m) is the complete elliptic integral of first kind.
Two distinguished limiting cases are m = 0, 1. The case m = 0 reduces to the constant potential q(x, 0) ≡ 1, which is solved trivially, and for which (L) = R ∪ [-i, i]. Conversely, as m → 1, q(x, 0) → sech x and K m → ∞. This problem was studied in Ref. [47], where it was shown that (L) comprises the real z axis plus a set of purely imaginary discrete eigenvalues uniformly distributed in the interval (-i, i). Moreover, it was shown that sech x is a reflectionless potential (giving rise to pure soliton solutions) if and only if = 1/N, with N ∈ N. This work connects these two limiting cases by characterizing the class of periodic potentials (4) for all 0 < m < 1. Figure 1 shows for = 1/20 and m = 0.7 the time evolution of the dn potential according to fNLS without noise (left), and with Gaussian noise (right).
The structure of the Lax spectrum can be studied via Bloch-Floquet theory using the monodromy matrix M(z) = -1 (x, z) (x + 2L, z), where (x, z) is a fundamental matrix solution of the scattering problem, and 2L is the minimal period of the potential. [Since (L) is invariant when the potential q(x, t ) evolves in time according to the fNLS Eq. (1), we will often omit the time dependence from all quantities below.] In turn, the Floquet discriminant is defined by (z) = tr M(z)/2. Bounded solutions of the scattering problem then exist for z ∈ C such that Im (z) = 0 and -1 Re (z) 1. Importantly, (z) is an entire function of z that is also Schwarz-symmetric: (z * ) = * (z). As a result, it is sufficient to consider the upper half plane, Im z 0. Since L is non-self-adjoint, (L) is not confined to the real z axis.
Nonetheless, the zero-level curves of Im (z) define a countable set of analytic arcs n . Along each arc, the requirement -1 Re (z) 1 then defines a spectral band. With these definitions, one can talk about bands as in a self-adjoint problem, the difference being that in the non-self-adjoint case the bands are not restricted to lie along the real z axis, but lie instead along the n .
The endpoints of the bands are periodic and antiperiodic eigenvalues of L, namely, the values z for which (z) = ±1, which give rise to periodic and antiperiodic eigenfunctions, respectively. Moreover, the real z axis is an infinitely long band, and any band that intersects the real z axis transversally is called a "spine." Any potential whose spectrum is comprised of a finite number of bands is called a "finite-band" potential [8,9,48]. These special potentials correspond to the finite-genus solutions of all equations in the Ablowitz-Kaup-Newell-Segur (AKNS) hierarchy.
In Ref. [40] it was proved that for the dn potential (4) and = 1/N, N ∈ N, the Lax spectrum of the ZS scattering problem comprises 2N Schwarz-symmetric bands along the interval (-i, i), and produces a genus 2N -1 solution of the fNLS equation (see Ref. [40] and the Appendix for details).
The spectrum of the dn potential (4) in the complex z plane for for m = 0.7 and = 1/7. is shown in Fig. 2(left). Figure 2(right) shows, as a function of the elliptic parameter m (horizontal axis), the location of the periodic eigenvalues (red curves) and antiperiodic eigenvalues (blue curves) along the imaginary z axis (vertical axis), together with the value of
1m (yellow curve) for m ∈ (0, 1). For the potential (4),
The periodic and antiperiodic eigenvalues of L along the imaginary z axis were computed using finite truncations of Eqs. (A10) (see Ref. [40] and the Appendix for details). Notice all gaps are closed when m = 0 and they open as m > 0 and remain open for all m ∈ (0, 1). In the singular limit m → 1, the band widths tend to zero, and the periodic and antiperiodic eigenvalues "collide" to form the point spectrum of L on the line.
Even though L is not self-adjoint, its spectrum in the limit → 0 can be effectively characterized using WKB methods [49,50], which allows one to obtain precise asymptotic estimates for the location of band edges, band widths, and gap widths (see Appendix). The results of the WKB analysis show that all bands lie along the finite segment = (-iq max , iq max ) in the complex z-plane. As → 0, however, the behavior of bands and gaps in the range o = (0, iq min ) differs from those in g = (iq min , iq max ). These two portions of the spectrum are those respectively located below (in white) and above (in gray) the yellow curve in Fig. 2(right) for each value of m. Inside g , the band widths decay exponentially in , while the gap widths decay algebraically. Hence, the ratio of band width to gap width tends to zero, and the bands accumulate according to some limiting density ϕ(z) found in Ref. [39] (see the Appendix for details). Conversely, inside o , it is the relative band gaps that tend to zero [cf. Eq. (A13) in the Appendix], giving rise to a continuous band in the limit. Thus, in the semiclassical limit, the spectrum of the elliptic potential ( 4) is compatible with the thermodynamic spectral scaling of a SG or a BG [38,39]. As m → 1, the range o tends to zero, and one gets a pure SG spectrum. Conversely, as m → 0, it is the range g that tends to zero, and one obtains the spectrum of the soliton condensate. For 0 < m < 1 we get the spectrum of a BG, which can be interpreted as a "composite," or partially condensed, SG. Next we characterize in detail the statistical properties of this gas.
We can derive the two key quantities in the spectral theory of soliton and breather gases, namely, the density of bands, and the normalized logarithmic band width. In particular, the density of bands ϕ(z) is (see the Appendix for details)
where x o (z) = K m for z ∈ (0, iq min ) and x o (z) > 0 is a simple turning point [at which dn(x; m) = |z|] for z ∈ (iq min , i). (Note that K m 0 dn(x; m) dx = π/2, which ensures that 1 0 ϕ(z) d|z| = 1. Note also that for z ∈ g , the above expressions follow directly from the analysis of Ref. [39].) Following Ref. [39], we can obtain the density of states (DOS) f (z) (see the Appendix for details) f (z) = ϕ(z)/(2K m ). (Note that, strictly speaking, f (z) is a "reduced" DOS. The "full" DOS for the bound-state potential is given by f (z) = f (η)δ(ξ ), where ξ = Re z, η = Im z, and δ(x) is the Dirac δ function; see Ref. [21].)
As m → 0, g collapses, and f (η) tends to the Weyl distribution of a soliton condensate-a critically dense SG (see the Appendix for details). Conversely, as m → 1, L → ∞, o collapses, and f (η) tends to zero, as appropriate for a rarified SG. Thus, as we already mentioned, the BG interpolating between these two extremes, when 0 < m < 1, can be viewed as a composite SG. The interpretation of a BG as a composite SG is also supported by the results of Ref. [51], where basic fNLS breathers were generated numerically with high accuracy by appropriately configuring N-soliton fNLS solutions with large N.
Kurtosis is an important characteristic of a random wave field, indicating its deviation from Gaussianity. With a slight abuse of conventional terminology here we will use the definition of kurtosis κ as the fourth normalized moment of the probability distribution for the wave amplitude |q|. With this definition the Gaussian wave field will have κ = 2; see, e.g., Ref. [52]. The spectral DOS allows one to evaluate the kurtosis of the SG wave field. In Ref. [21], the following general expression for the kurtosis of spatially uniform bound-state SG was obtained:
where the angle brackets denote ensemble average, η = Im z, and the overbar, denote average over the DOS,
i.e., the moments of the DOS f (η), where γ is the spectral support of the SG and where the last equality follows from Ref. [39], with c 1 = 1 4 and c 3 = 3 16 . Thus,
where κ 0 is the normalized fourth moment of the IC, namely,
This result is consistent with the kurtosis doubling obtained in Ref. [21] for the SG fission from partially coherent waves, with a crucial difference: here, κ 0 is computed from a purely deterministic IC. A calculation of the kurtosis for a periodic BG based on the results of Ref. [39] leads to the same expressions (8a) and (8b). Moreover, and importantly, a straightforward analysis of the above expression using Jensen's inequality (see the Appendix for details) allows one to conclude that any periodic SG or BG of the fNLS equation ( 1) generated by a deterministic real and even single-lobe potential q(x, 0) has kurtosis κ 2. Moroever, κ = 2 if and only if q(x, 0) is constant. Recall that κ > 2 implies heavy-tailed non-Gaussian statistics, which is an indication of the possible presence of rogue waves [53][54][55].
In particular, for the potential (4), a straightforward calculation allows us to obtain the following explicit expression for κ 0 :
where E m is the complete elliptic integral of the second kind [45].
Next we discuss the numerically computed time evolution of the fNLS solutions produced by the potential (4).
Figure 1 demonstrates how the addition of a small amount of zero-mean white Gaussian noise (with standard deviation σ = 10 -2 ) to the IC effectively randomizes the solution. as before, as well as their ensemble average (i.e., the kurtosis), clearly demonstrating the establishment of a statistically stationary integrable turbulence state. Finally, Fig. 3 (right) shows how the dependence of the kurtosis on the elliptic parameter m matches the theoretical prediction from Eqs. ( 6)-( 9), providing a strong validation of the analytical predictions.
In summary, we have presented an analytically tractable model describing a mechanism for the formation of integrable turbulence via BG fission of a semiclassical elliptic potential augmented by weak noise. Our analysis introduces a natural interpretation of a BG as a "composite SG," consisting of two distinct components: a regular SG plus a soliton condensate. The analytical model comprises a one-parameter family of such BGs, which interpolates between a pure soliton condensate (as m → 0) and a rarified SG (as m → 1) [56][57][58][59]. Intermediate values of m give rise to a mixed regime that interpolates between the above two extremes. We validated the theoretical results with direct numerical simulations of the BG fission and demonstrated the establishment, at large t of a statistically stationary integrable turbulence field characterised by the kurtosis value which is shown analytically and confirmed numerically to be twice as large as the fourth normalized moment of the initial elliptic potential, implying, in particular, the presence of rogue waves for all nonzero values of m.
We expect that similar results will hold for arbitrary real single-lobe periodic potentials. Our results also open up a number of interesting avenues for further research. Since the ZS problem is common to all equations of the AKNS hierarchy, the results of this work will generate SGs/BGs for all such equations, which although spectrally equivalent, exhibit qualitatively different dynamics. The study of their resulting SGs/BGs is threfore an interesting open question. For example, despite some recent work [60,61], the spectral theory of SGs for the focusing modified Korteweg-de Vries (mKdV) equation is still open, and we expect that the mKdV soliton gas phenomenology will be different than that of both the KdV and focusing NLS gases, potentially involving solitons of both polarities in the same gas. 4) is equivalent to studying the large-A limit of Eqs. (A1) and (A2).
Figure 4(left) depicts the potential ( 4) for various values of the elliptic parameter m ∈ (0, 1). Recall that the real period of the dn potential is 2K m , where K m = K (m) is the complete elliptic integral of the first kind [45]. Figure 4(right) shows the contours {z : Im (z) = 0} (thin black curves) together with the Lax spectrum (thick blue curves) for a generic (noninteger) value of the semiclassical parameter . Figure 5 shows the noise-augmented initial conditions (ICs) for two different values of m and Fig. 6 compares the Lax spectrum with and without noise, demonstrating that the noise does not affect it appreciably. In particular, the number and locations of the solitonic excitations remains unchanged.
In Ref. [40] we proved that the Lax spectrum of the focusing ZS scattering problem associated with Eq. (A1) and dnoidal potential (A2) and A ∈ Z is comprised of 2A Schwarz-symmetric bands along the imaginary axis of the spectral variable. As a result, the corresponding solutions of focusing NLS are special "finite-gap" solutions of genus 2A -1. Next, we provide some details on how the spectrum can be efficiently computed numerically.
PT-symmetric periodic Schrödinger potentials with real spectrum. When the potential is real-valued as with Eq. ( 4), (z) also possesses an additional symmetry: (-z * ) = (z), and thus the eigenvalues come in quartets, i.e., {z, z * , -z, -z * }. Moreover, the transformation w = v, where w(x, λ) = (w + , w -) T and
maps the ZS problem into the time-independent Schrödinger equation with a complex periodic potential, namely,
where
Recall that the Schrödinger equation with a periodic potential is referred to as Hill's equation [62]. If q is real and even, then the potentials are PT-symmetric: V ± (-x) = V * ± (x). In particular, for the elliptic potential Eq. ( 4), H ± is
xdn 2 (x) ± i m sn(x) cn(x). (A5)
Since (H + ) = (H -) it is enough to consider only H -. Below we discuss that, even though H ± are not Hermitian, all of their eigenvalues are real. A further simplification is obtained via the change of variable y = 2 am(x; m) [where am is the Jacobi amplitude], which maps Eq. (A5) into a complex perturbation of Ince's FIG. 6. Left: Numerically computed spectrum of the scattering problem for the dn potential with m = 0.9 via the Fourier-Hill method [63]. Right: Same but with noise added to the potential with μ = 0 and σ = 10 -2 . equation [62]: 4 2 [1m sin 2 (y/2)] w yy -2 (m sin y) w y + λ + (1m sin 2 (y/2)) + i m 2 sin y w = 0. (A6) Bloch-Floquet theory implies that any bounded solutions of Eq. (A6) can be written as w(y) = e iνy p(y), where ν ∈ R is the Floquet exponent and p(y + 2π ) = p(y), and can therefore be expanded in Fourier series, with Fourier coefficients given by a three-term recurrence relation. The eigenvalues of the ZS problem are then related to the (real or) complex values of λ for which the above ODE admits bounded solutions. Integer and half-integer values of ν yield respectively periodic and antiperiodic eigenfunctions. Three-term recurrence relation. Any solution of Eq. (A6) that is bounded for all y ∈ R can be represented by a Fourier series as w(y; λ) = e iνy n∈Z c n e iny , (A7)
with ν ∈ R. The coefficients {c n } n∈Z are given by the threeterm recurrence relation
where, for all n ∈ Z,
(A9c) Elliptic finite-band potentials. When = 1/N, some of the coefficents in Eqs. (A9) vanish, and as a result it is possible to decompose the doubly-infinite recurrence relation into two semi-infinite ones. Specifically, let
with ν = N/2 and j = 0 for the minus sign and ν = (1 -N )/2 and j = 1 for the plus sign. The periodic and antiperiodic spectrum of H -is the union of the spectra of B ± o and B ± ∞ . Specifically: when N is even, B - o and B - ∞ yield the periodic eigenvalues and B + o and B + ∞ the antiperiodic ones, and viceversa when N is odd. In Ref. [40] we proved that the eigenvalues of all four of these half-infinite matrices are real. This is the key to prove that, for any m ∈ (0, 1), the potential q in Eq. ( 4) for the focusing ZS scattering problem is finite-band if and only if = 1/N with N ∈ N. (The result is easily extended to N ∈ Z by phase invariance.) Moreover, if = 1/N, then q is a 2N-band potential, and
This implies that the spectral curve and the flow induced by each member of the Ablowitz-Kaup-Newell-Segur (AKNS) hierarchy has finite genus 2N -1. Importantly, the finite truncations of the matrices (A10) provide an efficient way to numerically compute the spectrum.
Semiclassical WKB analysis of the ZS problem. Next we consider the semiclassical limit of the spectrum, namely the limiting behavior of bands and gaps as → 0. It is useful to briefly recall the asymptotic analysis of the focusing ZS scattering problem via WKB methods from Ref. [49]. We emphasize that the results of this and the next few sections apply to a broad class of potentials, not just dn.
Suppose that q(x) is the 2L-periodic extension of a real, even, nonnegative single-lobe potential. Thus, q has one maximum and one minimum in (-L, L], which without loss of generality can be taken to be respectively at x = 0 and x = L. To avoid trivial cases, assume that q is not constant. Let q max = q(0), q min = q(L), (A12a) and
where z ∈ R ∪ iR, and x o (z) is a simple (real) turning point.
Without loss of generality we can limit ourselves to considering Im z 0 thanks to the Schwarz symmetry of the spectrum. Then, as → 0 (see Refs. [49,50]): (i) For z ∈ o ∪ R, where o = (0, iq min ), one has q 2 (x) + z 2 > 0 for x ∈ (-L, L). There are no turning points, and
where
q 2 (x) + z 2 dx. Therefore, in the WKB approximation, this range constitutes a single, continuous spectral band. This range corresponds to the region in white in Fig. 2.
(ii) For z ∈ g , where g = (iq min , iq max ), there are two real symmetric turning points ±x o , i.e., values at which q 2 (±x o ) + z 2 = 0. In this region (see Fig. 7), (z; ) ∼ cos(s 1 (z)/ ) cosh(s 2, (z)/ ), (A14)
where s 1 (z) = s(x o , z) and s 2, (z) = ln 2 -2s(-L, z). This range corresponds to the region in gray in Fig. 2.
(iii) For z ∈ (iq max , i∞) one has q 2 (x) + z 2 < 0 for x ∈ (-L, L). There are no turning points and
where
-q 2 (x)z 2 dx. Figure 7 shows a comparison between the predicted behavior of (z, ) and the numerically computed value.
Semiclassical limit of bands and gaps in g . When z ∈ (iq min , iq max ), Eq. (A14) implies that the bands are approximately centered at the roots {z n } of the equation
2 )π . Thus, the asymptotic number of bands for z ∈ (iq min , iq max ) is
where
and • is the floor function.
Next, consider the limit → 0 and n → ∞ satisfying 0 < n < J g /π and n is the band index. Denote the nth band width as
, where (z ± n ) = ±1, respectively. Similarly, denote the nth gap width as
together with Eq. (A14), one gets
where Im(z n+1 ) < Im(ξ n ) < Im(z n ). Finally, the band-to-gap ratio w n /g n → 0 as → 0 exponentially fast in this region.
Effective solitons in the fNLS equation. Recall that, for localized potentials (i.e., q(x, t ) → 0 rapidly as |x| → ∞), solitons are parameterized by the discrete eigenvalues of the ZS scattering problem. In the periodic problem, nonlinear excitations are considered to be "effective solitons" if the relative band width is less than some small threshold parameter [49,64,65], i.e., if W n < κ, where 0 < κ 1 and
Thus, rewriting the relative band width as
one can show that
Note that for z ∈ (iq min , iq max ) the relative band width W n is monotonic decreasing. Thus, the effective solitons are confined to the interval
where z s is the unique solution to W n = κ. For z ∈ (iq min , iq max ) let s 2 (z) = -s(-L, z) 0 (see Eq. (A12b)), then to leading order z s is given implicitly by
and the number of effective solitons is given by
Expanding s 2 (z) about z = iq min and evaluating at z s , for any κ > 0 it follows that z s → iq min as → 0. Thus, in the semiclassical limit the entire interval (iq min , iq max ) is comprised of an infinite ensemble of effective solitons. Semiclassical limit of bands and gaps in o . Recall that, when q = dn(x; m) with m ∈ (0, 1) and = 1/N with N ∈ N there are precisely N bands for z ∈ (0, iq max ). Hence, not all bands become effective solitons in the limit → 0 (see Fig. 2). For z ∈ o , following similar arguments as for z ∈ g , the WKB approach of Ref. [49] implies
Thus, at z * = iq min = i dn(K m ; m), which is the boundary between the intervals g = (iq min , iq max ) and o = (0, iq min ), there is a transition from exponentially decaying band widths (giving rise to a soliton gas), to algebraically decaying band widths (giving rise to a soliton condensate). Thus, we have the physically realistic scenario of a generalized breather gasnamely, a soliton gas on a soliton condensate background, which can equivalently be interpreted as a composite, or partially condensed, SG.
The number of bands with support in o can be computed by noting that Eq. (A13) holds with exponentially small corrections at (z) = ±1. Thus, bands in this region are approximately centered at the roots z n of the equation
2 )π . (The existence of a small gap for each z n ∈ o was proven in Ref. [40].) The asymptotic number of bands for z ∈ (0, iq min ) is
where
Periodic soliton and breather fNLS gases. The semiclassical WKB asymptotics of the periodic ZS problem from Ref. [49] was utilized in Ref. [39] to introduce periodic fNLS SGs and BGs. They could be viewed as the semiclassical limit of periodic fNLS solutions after a sufficiently long time evolution. Under the assumption of an even, nonnegative, continuous, one-hump 2L-periodic initial potential q(x) with the maximum value located at x = 0, the density of band centers ϕ(z) and the scaled logarithmic band width ν(z) for a periodic fNLS BG/SG were calculated as
where q min 0 is the minimum of q(x) attained at x = L, i.e., q(L) = q min . The case of q min > 0 corresponds to a BG, whereas the case q min = 0 corresponds to a SG. In both SG and BG cases the DOS f (z) was calculated to be f (z) = rϕ(z), where
The interpretation of a BG as a composite SG can be derived from the results of Ref. [39], where one can observe that the DOS f (z) of a periodic BG (q min > 0) in the region g coincides in g with the DOS of SG obtained by replacing q(x) with a modified potential q(x) by means of making an infinitesimally narrow cut at x = L that replaces the minimum value of q(L) = q min by q(L) = 0. In the region o , the DOS corresponding to q(x) is still given by the expression for rϕ(z) in Eq. (A28), but where now one has to take q min = 0.
The average conserved quantities I m , m = 1, 2, . . . , for such gases were rigorously obtained in Ref. [39]. In particular, I k = 0 for any even k ∈ N, while, for any odd k ∈ N,
where, for any odd value of m, d m is defined as
Spectral characteristics of the composite soliton gas. Recall that the density of bands ϕ(z) and scaled logarithmic band width ν(z) are two key quantities in the spectral theory of soliton gases. Following Eq. (A27), we have that, for z ∈ g , and the dn potential ( 4)
with J is a suitable normalization constant, determined below. Thus we can easily express the spectral scaling function as
TABLE I. Numerical simulation parameters and results, where N Fourier is the number of Fourier modes used, M is the number of simulations (each of which corresponds to an independent noise realization), and T ∞ is the numerically observed approximate thermalization time. m κ o κ ∞ T ∞ M N Fourier 0.05 1.007 1.92 4 500 1024 0.3 1.016 1.96 4.5 400 1024 0.5 1.059 2.07 4 400 2048 0.7 1.17 2.37 2.5 200 4096 0.8 1.29 2.66 2.5 200 4096 0.9 1.53 3.022 3.0 400 4096 Similarly, for z ∈ o , as it was discussed above, we have
(In o , the spectral scaling function σ (z) is zero as required for a condensate.)
Finally, ensuring that the integral of the density of bands ϕ(z) over the whole support (0, i) is one, one can then show that
Limits m → 0 and m → 1. Next we discuss the density of bands formulas in the limits m → 0 + , and m → 1 -.
Recall that dn(x; 0) ≡ 1 and K (0) = π/2. Thus, q min → 1 as m → 0 + and Eq. (A33a) has support z ∈ (0, i). Moreover, J o → π as m → 0 + . Thus, which corresponds to Weyl's distribution as expected (see Refs. [33,38]).
Next, recall that dn(x; 1) ≡ sech(x) and K (m) → ∞ as m → 1 -. Thus, q min → 0 as m → 1 -and Eq. (A31a) has support z ∈ (0, i). Moreover, J g → π as m → 1 -, and dn -1 (|z|) → arcsech(|z|) as m → 1 -. Thus, in this limit we get a uniform distribution:
which matches well known results of semiclassical distribution of discrete spectrum for the sech potential [57,59]. Nonlinear dispersion relations. Together with the set + = g ∪ o , the spectral scaling function σ (z) determines the integral equation (the first nonlinear dispersion relation, NDR) for the DOS f (z) for fNLS SGs:
(Note that in the portion of the spectrum corresponding to a soliton condensate, i.e., o , one has σ (z) = 0.) The second NDR involves the spectral flux density, but this quantity is zero in our case since all the nonlinear excitations have zero velocity. It was proven in Ref. [39] that, for periodic SGs and BGs, f (z) = rϕ(z) solves the first NDR for soliton/breather gases respectively. In the case of the dn potential (4),
As a result, expressions (A31a) and (A33a) also yield the DOS f (z) for the dn potential (4).
Kurtosis and Jensen's inequality. In the main text we have shown that the kurtosis κ of a fully developed SG generated by a deterministic real and even periodic single-lobe initial condition q(x, 0) is given by κ = 2κ 0 , where
By employing the rescaling x = x/L and q(x, 0) = q(L x, 0), we can rewrite κ 0 as
Then Jensen's inequality [45] implies that κ 0 1, with the equality holding only if q(x, 0) is constant. Numerical methods. We performed several numerical experiments related to the kurtosis of randomly generated solutions to the fNLS equation. In particular, we studied the time evolution of fNLS in the semiclassical limit with elliptic FIG. 10. Same as Fig. 9, but for m = 0.7.
IC plus a complex perturbation of independent and identically distributed normal random variables with mean μ = 0 and standard deviation σ = 10 -2 .
Due to modulational instability of the fNLS equation, the numerical computation of solutions is a delicate topic, especially in the semiclassical limit. To this end, all the numerical simulations of the fNLS equation (1) presented in this work were performed using an eighth-order split-step method [66]. The spatial and temporal discretization parameters were chosen to be small enough that none of the results presented are affected by numerical errors. Moreover, to ensure numerical accuracy the isospectral property of the scattering data was confirmed at several data points as the simulated solution evolved in time. The Lax spectrum was computed using the Floquet-Fourier-Hill method [63]. Finally, to reduce the number of simulations needed a double-averaging technique was used to compute the kurtosis. That is, at each time the spatial average of |q| 2n , n = 1, 2 was computed, and then the ensemble average was computed.
Further simulation results. For consistency, the numerical experiments presented in this work related to the kurtosis have = 1/20, and σ = 10 -2 fixed for all values of m ∈ (0, 1) considered. Table I provides key parameters used in the sim-ulations as well as some of the key quantities related to the kurtosis experiments.
Figure 8 illustrates the time evolution of the fourth normalized moment (i.e., kurtosis) of the solution for m = 0.05, 0.3, 0.5, 0.7, 0.8, and 0.9. The gray region is the envelope of 100 ensemble realizations. The red curves are ensemble averages of M realizations. Notice for each m the ensemble average of the kurtosis settles at large times (see Table I for approximate thermalization times).
Finally, to further illustrate the effective thermalization of the SG at large times, Figs. 9-11 show the numerically computed solution of fNLS in the semiclassical limit are provided. In each case, the left column depicts the time evolution of the dn potential (4) while the right column depicts the time evolution of the dn potential (4) with a small complex-valued random perturbation. Further, the top row depicts a density plot while the bottom row depicts the solution at a particular time. In Fig. 9 we have m = 0.05 which is a small sinusoidal perturbation of the constant background studied in Ref. [34]. In