# Perturbed Nonlinear Evolution of Optical Soliton Gases: Growth and Decay in Integrable Turbulence

Loic Fache, François Copie, Pierre Suret, Stéphane Randoux

► To cite this version:

Loic Fache, François Copie, Pierre Suret, Stéphane Randoux. Perturbed Nonlinear Evolution of Optical Soliton Gases: Growth and Decay in Integrable Turbulence. Physical Review Letters, 2025, 135 (15), pp.157201. 10.1103/hwpr-628r. hal-05382173

HAL Id: hal-05382173

https://hal.science/hal-05382173v1

Submitted on 25 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Perturbed Nonlinear Evolution of Optical Soliton Gases: Growth and Decay in Integrable Turbulence

Loic Fache, François Copie, Pierre Suret, and Stéphane Randoux*

Univ. Lille, CNRS, UMR 8523 - PhLAM - Physique des Lasers Atomes et Molécules, F-59 000 Lille, France

(Dated: November 25, 2025)

We present optical fiber experiments investigating the perturbed, non-integrable evolution of soliton gases (SGs) under weak linear damping and gain. By measuring the amplitude and phase of the optical field in a recirculating loop, we determine the spectral distribution of SGs at various propagation distances. We demonstrate that a SG, initially prepared as a soliton condensate with a Weyl distribution, undergoes significant changes in its spectral distribution due to dissipation. Specifically, weak gain leads to the emergence of a bichromatic SG, while weak damping results in a SG with a semicircular spectral distribution at long propagation distances. These observed dissipation-driven changes in the SG's spectral characteristics are not explained by existing hydrodynamic or kinetic theories.

Integrable turbulence (IT) refers to both a physical phenomenon and a theoretical framework for studying the dynamics of nonlinear random waves in systems governed by integrable equations, such as the Korteweg–de Vries (KdV) equation or the one-dimensional nonlinear Schrödinger equation (1D-NLSE) [1–9]. In contrast to wave turbulence [10–13], which addresses the behavior of out-of-equilibrium, nonlinear random dispersive waves in non-integrable systems, the theoretical framework underlying IT is the inverse scattering transform (IST) method [14, 15]. This distinction provides a relevant perspective on the analysis of large-scale, emergent hydrodynamic phenomena in integrable systems, where soliton gas (SG) serves as a central concept [16–20].

SG can be defined as an infinite ensemble of interacting solitons with randomly distributed amplitudes, velocities, and positions [16, 17, 20]. The integrable dynamics, constrained by the existence of an infinite number of conserved quantities, prevent SGs from achieving relaxation to the classical Gibbs ensembles. Instead SGs exhibit local nonthermal stationary states, the so-called generalized Gibbs ensembles, which play a fundamental role in generalized hydrodynamics (GHD), the hydrodynamic theory of many-body quantum and classical integrable systems [21–25]. Currently, the out-of-equilibrium evolution of SGs represents an active area of research, explored using the spectral kinetic theory of SGs and GHD [18, 20, 22, 24].

Despite several experimental observations of SGs across various physical systems [3, 26–31], only a few experimental studies have questioned the validity of the spectral kinetic theory of SGs in a quantitative way [32–34]. Several experiments have demonstrated that while perturbative higher-order effects are inevitably present, they break integrability to such a negligible degree that the predictions of kinetic theory still align closely with experimental results [32, 33]. In contrast, other experiments have demonstrated that even small perturbative

higher-order effects can significantly alter integrable dynamics, leading to unexpected behaviors not anticipated at the theoretical level [34]. For example, a KdV soliton condensate has been shown to emerge spontaneously as a result of weak dissipation in a nonlinear electrical transmission line [34].

In this Letter, we present experiments investigating the perturbed, non-integrable evolution of optical SGs, which reveal unexpected and remarkable changes in their spectral, IST properties. Using a recirculating optical fiber loop that allows precise control of either optical gain or losses, we examine the growth and decay of IT in SG experiments, accurately described by the focusing 1D-NLSE with a weak linear gain or damping term [35–37]. Our experiments reveal that the SG, initially prepared as a soliton condensate characterized by the Weyl distribution, undergoes significant changes in its spectral distribution due to dissipation. We find that a small amount of optical gain leads to the emergence of a bichromatic SG. On the other hand, small damping is found to alter the geometric shape of the IST spectrum, resulting in a semicircular spectral configuration at long evolution time. The dissipation-driven changes observed experimentally in the spectral characteristics of the optical SGs are not explained by existing hydrodynamic or kinetic theories of SG.

Our experimental setup is depicted schematically in Fig. 1. It comprises a recirculating fiber loop constructed from ~8 km of single-mode fiber (SMF), which is closed upon itself using a 90/10 fiber coupler. This coupler is configured to recirculate 90% of the optical power. The optical signal circulates in the clockwise direction, and at each round trip, 10% of the circulating power is extracted and directed to a photodetector (PD) connected to a fast oscilloscope with a sampling rate of 160 GSa/s and a bandwidth of 65 GHz. The overall detection bandwidth of the oscilloscope and the PD is 32 GHz. The periodic extraction of light at each round trip within the recirculating fiber loop enables stroboscopic monitoring of the wavefield's evolution every 8 km. The signals recorded by the oscilloscope are subsequently processed numerically to generate space-time diagrams that illustrate the

<footnote>* stephanerandoux@univ-lille.fr</footnote>
dynamics of the wavefield over hundreds of round trips inside the fiber loop [33, 38, 39].

<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental setup. A SG, initially in the form of a 100 ns-long flat-top pulse at 1550 nm, perturbed by optical noise, propagates through a recirculating fiber loop where the optical gain or loss rate can be precisely controlled through Raman amplification using a 1450 nm pump laser. Measurement of the amplitude and the phase of the optical signal is made using a heretodyne technique where the optical signal at the output of the fiber loop is beated against a single-frequency laser at 1550.361 nm.</figcaption></figure>

A typical experiment involves circulating a long, flat-top pulse (~100 ns duration) for 400 to 800 round trips in the fiber loop, corresponding to total propagation distances of ~3200 km to ~6400 km. The flat-top pulse is generated using a fast electro-optic modulator (EOM) and is perturbed by a small amount of optical noise introduced by an Erbium-doped amplifier prior to the pulse's injection into the fiber loop, see Fig. 1. As a result of the process of modulation instability, the perturbed flat-top pulse evolves into a fully randomized bound state SG [29, 33, 40–43].

The mean gain or loss of the optical signal over hundreds of round trips in the fiber loop can be precisely controlled through Raman amplification, using a counter-propagating 1450 nm pump laser coupled into and out of the loop via wavelength division multiplexers (WDMs), as shown in Fig. 1. By carefully adjusting the optical power of the pump laser using a variable optical attenuator (VOA), we performed experiments where the mean power of the optical SGs is set to exponentially increase or decay at a given rate, allowing us to explore both growing and decaying stages of IT [35, 36].

Capturing the changes in the IST spectrum experienced by SGs during propagation is a significant challenge in optics, as it requires simultaneous measurement of both the amplitude and phase of the optical field. In this work, we address this challenge using a heterodyne technique. Specifically, the optical signal at the output of the recirculating fiber loop is mixed with a single-frequency reference field that is detuned by approximately 45 GHz (equivalently 0.361 nm) from the central frequency of the optical signal. As illustrated in Fig. 1, the optical signal and the detuned reference field are combined using a 50/50 fiber coupler, with the output directed to a fast photodetector (PD). The reconstruction

of the optical field's amplitude and phase employs conventional signal processing techniques that are described in the Supplemental Material [44].

Focusing first on experiments involving decaying IT, Fig. 2(b) illustrates the space-time evolution of the optical wavefield power, experimentally observed as the mean power $P(z)$ of the flat-top pulse decreases from $P_0 \sim 11.5$ mW to ~2 mW over a propagation distance of ~3200 km (see Fig. 2(a)). At a qualitative level, the experimentally reconstructed space-time diagram over a narrow 5 ns time window closely resembles the space-time patterns typically associated with the noise-driven modulation instability of plane waves in systems governed by the focusing 1D-NLSE [39, 45–50].

As shown in Fig. 2(e), the dynamical features observed in our experiment are quantitatively well described by numerical simulations of the following 1D-NLSE with a small linear damping term:

$$i \frac{\partial A}{\partial z} = \frac{\beta_2}{2} \frac{\partial^2 A}{\partial T^2} - \gamma |A|^2 A - i \frac{\alpha_{\text{eff}}}{2} A. \quad (1)$$

$A(z, T)$ represents the complex envelope of the electric field that slowly varies in physical space $z$ and time $T$. The Kerr coefficient of the fiber is $\gamma = 1.3 \text{ W}^{-1}\text{km}^{-1}$. The group velocity dispersion coefficient is $\beta_2 = -22 \text{ ps}^2 \text{ km}^{-1}$. $\alpha_{\text{eff}} > 0$ represents the effective power decay rate of the circulating field. In the experiment presented in Fig. 2(a)(b), its value is $\alpha_{\text{eff}} \sim 5.8 \times 10^{-4} \text{ km}^{-1}$ or equivalently ~0.0025 dB/km.

Although the space-time diagrams in Figs. 2(b) and 2(e) might suggest that dissipation has a negligible effect on the noise-driven destabilization of the initial flat-top pulse, we now show that the discrete IST spectra of the optical SG reveal distinct features introduced by dissipation to the dynamics. For the computation of the discrete IST spectra, we first introduce the following dimensionless form of the focusing 1D-NLSE [37, 48]

$$i\psi_t + \psi_{xx} + 2|\psi|^2\psi + i\epsilon\psi = 0 \quad (2)$$

which is obtained from Eq. (1) using the following transformations: $\psi = A/\sqrt{P_0}$, $x = T\sqrt{\gamma P_0/|\beta_2|}$, $t = \gamma P_0 z/2$, $\epsilon = \alpha_{\text{eff}}/(\gamma P_0)$ [33].

Within the IST formalism, the non-self-adjoint Zakharov-Shabat eigenvalue problem that is associated with Eq. (2) for $\epsilon = 0$ reads [51]

$$\hat{\mathcal{L}}\Phi = \lambda\Phi, \quad \hat{\mathcal{L}} = \begin{pmatrix} -i\partial_x & -i\psi \\ -i\psi^* & i\partial_x \end{pmatrix}. \quad (3)$$

$\Phi(x, \lambda)$ is a vector wave function. $\lambda \in \mathbb{C}$ represent the eigenvalues composing the discrete spectrum associated with the soliton content of the field $\psi = \psi(x, t)$ that is measured at some given evolution time $t$ (equivalently at some given propagation distance $z$ in the experiment). Each soliton in a soliton gas, consisting of an ensemble of $N$ solitons, is characterized by a discrete eigenvalue $\lambda_j = \xi_j + i\eta_j (j = 1, \dots, N)$ of the spectrum of the linear
<figure><img src="image_3.png" /><figcaption>FIG. 2. Experiments and numerical simulations showing the evolution of an optical SG in a recirculating fiber loop in the presence of small damping. (a) The black points represent the measured power of the SG as a function of the propagation distance $z$. The red line represents the evolution of the optical power computed in numerical simulations. The green points (resp. line) represent the number of discrete eigenvalues $\lambda_i$ measured in the experiment (resp. computed in numerical simulations) as a function of $z$. (b) Space-time evolution of the optical power of the SG measured in the experiment. (c) DOS of the optical SG measured at $z = 1920$ km in the upper complex plane. (d) DOS of the initial condition. (e) Same as (b) but from numerical simulations of Eq. (1) with $\gamma = 1.3 \text{ W}^{-1}\text{km}^{-1}$, $\beta_2 = -22 \text{ ps}^2 \text{ km}^{-1}$, $\alpha_{\text{eff}} \sim 5.8 \times 10^{-4} \text{ km}^{-1}$, $P(z=0) = P_0 = 11.5 \text{ mW}$. (f), (g) Same as (c), (d) but computed from simulations reported in (e). (h), (i) Same as (f), (g) but including noise and finite bandwidth effects in the heterodyne measurement, see Supplemental Material [44] for details.</figcaption></figure>

operator $\hat{\mathcal{L}}$, where $\xi_j$ and $\eta_j$ parameterize the velocity and amplitude of each soliton, respectively. Importantly, integrable dynamical evolution (implying $\epsilon = 0$ in Eq. (2)) entails an isospectral property, meaning that all the discrete eigenvalues $\lambda_j$ remain invariant with respect to the evolution variable $t$.

The discrete IST spectra of the optical SG are computed from the normalized experimental complex fields $\psi = A/\sqrt{P_0}$ by solving numerically Eq. (3) using the Fourier collocation method described in ref. [51]. The output of this numerical calculation is an ensemble of several hundreds of discrete, complex-valued eigenvalues $\lambda_j$ that are located in a bounded region of the upper complex plane. Given the large number of discrete eigenvalues that are found from nonlinear spectral analysis, we follow the approach used in the kinetic theory of SG and compute the so-called density of states (DOS) $f(\lambda; x, t)$ of the SG.

The DOS $f(\lambda; x, t)$, where $\lambda = \xi + i\eta$, represents the density of soliton states in the phase space, i.e. $fd\xi d\eta dx$ is the number of solitons contained in a portion of SG with the complex spectral parameter $\lambda \in [\xi, \xi + d\xi] \times [\eta, \eta + d\eta]$ over the space interval $[x, x + dx]$ at time $t$ (corresponding to the position $z$ in the fiber loop experiment). Given that the optical SG is spatially homogeneous, the DOS represents the probability density function of the complex-valued discrete eigenvalues normalised in such a way that $\int_{-\infty}^{+\infty} d\xi \int_{0}^{+\infty} d\eta f(\lambda) = N/\Delta x$, where $N$ represents the number of eigenvalues found in the upper complex plane and $\Delta x$ represents the spatial extent of the gas [19, 30].

In practice, the experimental DOS is computed not from a single flat-top pulse, but from three identical flat-top pulses, each evolving independently within the fiber loop during the same experimental run. Assuming minimal variation in the IST spectra over short propagation

distances, the DOS is calculated by averaging the IST spectra of the three pulses collected over a 80 km propagation span, as highlighted by the grey-shaded regions in Fig. 2(a). The statistical computation of the DOS involves an ensemble of $\sim 23200$ (resp. $\sim 18300$) discrete eigenvalues at $z = 0$ km (resp. $z = 1920$ km).

Fig. 2(d) shows the DOS of the flat-top pulse perturbed by small optical noise that is taken as initial condition in the experiment ($z = 0$ km). As it can be easily anticipated from previous works [40], we find that the measured DOS is very well fitted by the Weyl distribution $f_W(\eta) = \eta/(\pi\sqrt{1-\eta^2})$ located along the vertical imaginary axis $\xi = 0$ (see also Fig. 4(a)). Due to dissipation, the isospectrality condition characteristic of integrable dynamics is not expected to be satisfied in our experiment. However, since the damping effects are small ($\epsilon \approx 0.038$ in Eq. (2)), IST analysis remains a relevant tool for examining the evolution of the optical SG. Fig. 2(c), computed at $z \sim 1900$ km, shows that optical damping transforms the spectral support of the DOS in the upper complex plane, reshaping the vertical linear support of the initial Weyl distribution into a semicircular one. It is important to note that the concept of a circular spectral support for the DOS of a SG was introduced mathematically in Ref. [18], but without any expectation that this class of SGs could be observed experimentally nor that it could emerge from a non-integrable evolution. It is also worth noting that our experiments demonstrate a transition from a soliton condensate to a circular SG, rather than to a circular condensate.

As shown in Fig. 2(f) and 2(g), the dissipation-driven changes of the geometrical support of the DOS are also found in numerical simulations of Eq. (1) that are made with the experimental parameters. It is noteworthy that the DOS obtained from numerical simulations of Eq. (1), shown in Fig. 2(f), exhibits a perfectly symmetric shape
<figure><img src="image_4.png" /><figcaption>FIG. 3. Experiments and numerical simulations showing the evolution of an optical SG in a recirculating fiber loop in the presence of small gain. (a) The black points represent the measured power of the SG as a function of the propagation distance z. The red line represents the evolution of the optical power computed in numerical simulations. The green points (resp. line) represent the number of discrete eigenvalues $\lambda_i$ measured in the experiment (resp. computed in numerical simulations) as a function of z. (b) Space-time evolution of the optical power of the SG measured in the experiment. (c) DOS of the optical SG measured at z = 6320 km in the upper complex plane. (d) DOS of the initial condition. (e) Same as (b) but from numerical simulations of Eq. (1) with $\alpha_{eff} \sim -1.1 \times 10^{-4} \text{ km}^{-1}$, $P(z = 0) = P_0 = 2 \text{ mW}$. (f), (g) Same as (c), (d) but computed from simulations reported in (e). (h), (i) Same as (f), (g) but including noise and finite bandwidth effects in the heterodyne measurement, see Supplemental Material [44] for details.</figcaption></figure>

that is not found in the experiments, as seen in Fig. 2(c). This slight discrepancy between the experimental and numerical DOS, along with the blurrer appearance of the experimental DOS, is attributed to noise and finite bandwidth effects in the heterodyne measurement of the optical field's phase and amplitude, as illustrated in Fig. 2(h) and 2(i) and discussed more in details in Supplementary Material [44]. Let us also emphasize that, with the constant loss rate, the optical SG never attains a steady state and the circular support of the DOS continuously undergoes homothetic evolution; see the Supplementary Material [44] for more details on the DOS evolution.

Next, we examine another dissipative experiment in which the optical SG undergoes a linearly amplified evolution instead of a damped one. Increasing the power of the 1450 nm pump laser, we now consider a situation in which the mean optical power of the SG increases from ~2 mW to ~4 mW over ~6300 km, as shown in Fig. 3(a). Fig. 3(b) reveals that the space-time evolution of the optical SG looks qualitatively different from the situation with damping reported in Fig. 2(b). Following the initial destabilization of the flat-top pulse, individual nonlinear structures with small relative velocities are observed to form and persist within the fiber loop. At a qualitative level, the observed space-time evolution shows remarkable similarities to that reported in previous numerical studies on growing IT, even though these studies considered a completely different initial condition [36].

Figures 3(c) and 3(d) illustrate that the DOS of the SG, initially in the form of the Weyl distribution, essentially evolves into two distinct clustered distributions. One of these clusters is centered around $\lambda_0 \simeq 1.8i$, while the other is centered around $\lambda_1 \simeq 0.3i$. Here, the statistical computation of the DOS involves an ensemble of ~10700 (resp. ~12100) discrete eigenvalues at z = 0 km (resp. z = 6320 km). The fact that $Re(\lambda_0) \sim Re(\lambda_1) \sim 0$

is the spectral signature of all solitons in the amplified SG exhibiting small velocities, in stark contrast to the scenario depicted in Fig. 2(c), where damping was shown to broaden the soliton velocity distribution. The difference between $Im(\lambda_0) \simeq 1.8$ and $Im(\lambda_1) \simeq 0.3$ indicates that the SG formed at z ~ 6300 km is predominantly composed of two distinct species, characteristic of a bichromatic SG [17, 32].

Fig. 3(e) shows that the space time evolution observed in the experiment is well reproduced by numerical simulations of Eq. (1) for $\alpha_{eff} \sim -1.1 \times 10^{-4} \text{ km}^{-1}$, corresponding to $\epsilon = -0.042$. Fig. 3(f)-(i) show that the features observed experimentally in terms of evolution of the DOS are also found in numerical simulations that include noise and finite bandwidth effects due to the heterodyne measurement of the optical field's phase and amplitude, see Supplementary Material [44] for details.

In our experiments, dissipation significantly alters the mass $M(t) = \int_{-\infty}^{+\infty} |\psi(x, t)|^2 dx$, or equivalently, the mean power, $P(z) = \int_{-\infty}^{+\infty} |A(T, z)|^2 dT$, of the optical SG, see Fig. 2(a) and 3(a). However, the total number of discrete eigenvalues, $N(z)$, does not vary in proportion to the changes in mass. In Fig. 2(a) the mass decreases by a factor ~6, while $N$ changes by only ~20% during the evolution. Similarly, in Fig. 3(a), the mass approximately doubles but $N$ changes by only 10% during the evolution. This indicates that the changes in the mass of the optical SGs are not primarily caused by the creation or annihilation of soliton states, but rather by the rearrangement of the DOS in the upper complex plane.

To substantiate this point in a more quantitative way, we consider the expression of the mass of SG given by the spectral kinetic theory:

$$M = 4 \int_{0}^{+\infty} \eta \tilde{f}(\eta) d\eta = \int_{0}^{+\infty} \rho(\eta) d\eta. \quad (4)$$
<figure><img src="image_5.png" /></figure>

FIG. 4. Experimental histograms showing the spectral density of mass, $\rho(\eta)$, along the vertical imaginary axis $\eta$. (a) The green (resp. grey) histogram represents the initial distribution in the experiment with linear damping (resp. amplification) of the SG. The red dashed line represents the spectral density of mass $\rho_W(\eta) = 4\eta f_W(\eta)$ associated with the Weyl distribution $f_W(\eta) = \eta/(\pi\sqrt{1-\eta^2})$. (b) Same as (a) but at $z \sim 1920$ km (green histogram) and $z \sim 6320$ km (grey histogram).

where $\tilde{f}(\eta) = \int_{-\infty}^{+\infty} d\xi f(\xi, \eta)$ represents the DOS $f(\xi, \eta)$ that has been integrated over the real axis $\xi$ [18, 19]. $\rho(\eta) = 4\eta \tilde{f}(\eta)$ represents the spectral density of mass with respect to $\eta$.

Fig. 4(a) illustrates that the spectral density of mass, $\rho(\eta)$, for the initial bound-state SG closely matches the spectral density $\rho_W(\eta) = 4\eta f_W(\eta)$, which corresponds to the Weyl distribution. In contrast, Fig. 4(b) shows how experimental damping significantly alters $\rho(\eta)$, transforming it into a bell-shaped distribution centered around $\eta \simeq 0.4$ (green histogram). On the other hand, linear amplification of the SGs results in the development of an asymmetric distribution of the spectral density (grey histogram in Fig. 4(b)), with its peak shifted to approximately $\eta \simeq 1.8$. The mass associated with the cluster of eigenvalues centered around $\lambda_1 \simeq 0.3i$ constitutes only a few percent of the total mass of the formed bichromatic SG, see Supplementary Material [44] for details.

In conclusion, we have presented optical fiber experiments exploring the perturbed evolution of SGs under the influence of small linear damping and gain. Our experiments reveal that dissipation-induced changes in

the mass (equivalently, in the mean power) of optical SGs are predominantly driven by a redistribution of the DOS in the upper complex plane, rather than by the creation or annihilation of soliton states. This feature contrasts with recent findings in the perturbed evolution of KdV SGs, see ref. [34] Previous studies investigating dissipative evolution in NLS systems have primarily focused on ensembles containing only a few solitons, employing approaches such as numerical simulations or perturbed IST theory [52–56]. We hope that our experiments, which explore large ensembles of interacting solitons with random characteristics, will inspire theoretical investigations using GHD, possibly building on previous studies of one-dimensional Bose gases with atom losses [57, 58]. Furthermore, our experiments demonstrate that higher-order perturbative effects can be exploited to manipulate the spectral (IST) characteristics of SGs, offering new experimental opportunities.

## ACKNOWLEDGMENTS

We thank Dmitry Agafontsev, Thibault Bonnemain, Thibault Congy, Gennady El, Andrey Gelash, Giacomo Roberti, and Alexander Tovbis for their collaboration and the discussions we have had over the years on SG-related topics. This work has been partially supported by the Agence Nationale de la Recherche through the SOGOOD (ANR-21-CE30-0061) projects, the LABEX CEMPI project (ANR-11-LABX-0007), the Ministry of Higher Education and Research, Hauts de France council and European Regional Development Fund (ERDF) through the Nord-Pas de Calais Regional Research Council and the European Regional Development Fund (ERDF) through the Contrat de Projets Etat-Région (CPER Photonics for Society P4S). The authors would like to thank the Centre d'Etudes et de Recherche Lasers et Application (CERLA) for technical support and the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Emergent phenomena in nonlinear dispersive waves, where work on this paper was undertaken.

[1] V. E. Zakharov, Stud. Appl. Math. **122**, 219 (2009).

[2] J. M. Soto-Crespo, N. Devine, and N. Akhmediev, Phys. Rev. Lett. **116**, 103901 (2016).

[3] I. Redor, H. Michallet, N. Mordant, and E. Barthélemy, Phys. Rev. Fluids **6**, 124801 (2021).

[4] S. Randoux, P. Walczak, M. Onorato, and P. Suret, Phys. Rev. Lett. **113**, 113902 (2014).

[5] G. Michel, F. Bonnefoy, G. Ducrozet, G. Prabhudesai, A. Cazaubiel, F. Copie, A. Tikan, P. Suret, S. Randoux, and E. Falcon, Phys. Rev. Fluids **5**, 082801 (2020).

[6] T. Congy, G. A. El, G. Roberti, A. Tovbis, S. Randoux, and P. Suret, Phys. Rev. Lett. **132**, 207201 (2024).

[7] S. Randoux, F. Gustave, P. Suret, and G. El, Phys. Rev. Lett. **118**, 233901 (2017).

[8] D. S. Agafontsev and V. E. Zakharov, Nonlinearity **28**, 2791 (2015).

[9] Z.-Y. Sun, X. Yu, and Y.-J. Feng, Phys. Rev. E **108**, 054211 (2023).

[10] S. Nazarenko, *Wave Turbulence*, Lecture Notes in Physics (Springer Berlin Heidelberg, 2011).

[11] A. Picozzi, J. Garnier, T. Hansson, P. Suret, S. Randoux, G. Millot, and D. Christodoulides, Physics Reports **542**, 1 (2014).

[12] E. Falcon and N. Mordant, Annual Review of Fluid Mechanics **54**, 1 (2022).

[13] S. Galtier, *Physics of Wave Turbulence* (Cambridge University Press, 2022).

[14] P. G. Drazin and R. S. Johnson, *Solitons: an introduction*, Vol. 2 (Cambridge University Press, 1989).

[15] M. Ablowitz, *Nonlinear Dispersive Waves: Asymptotic*
*Analysis and Solitons*, Cambridge Texts in Applied Mathematics (Cambridge University Press, 2011).

[16] V. E. Zakharov, Sov. Phys.-JETP **33**, 538 (1971).

[17] G. A. El and A. M. Kamchatnov, Phys. Rev. Lett. **95**, 204101 (2005).

[18] G. El and A. Tovbis, Phys. Rev. E **101**, 052207 (2020).

[19] G. A. El, Journal of Statistical Mechanics: Theory and Experiment **2021**, 114001 (2021).

[20] P. Suret, S. Randoux, A. Gelash, D. Agafontsev, B. Doyon, and G. El, Phys. Rev. E **109**, 061001 (2024).

[21] R. Koch, J.-S. Caux, and A. Bastianello, Journal of Physics A: Mathematical and Theoretical **55**, 134001 (2022).

[22] T. Bonnemain, B. Doyon, and G. El, Journal of Physics A: Mathematical and Theoretical **55**, 374004 (2022).

[23] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Phys. Rev. Lett. **117**, 207201 (2016).

[24] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X **6**, 041065 (2016).

[25] B. Doyon, S. Gopalakrishnan, F. Møller, J. Schmiedmayer, and R. Vasseur, Phys. Rev. X **15**, 010501 (2025).

[26] I. Redor, E. Barthélemy, H. Michallet, M. Onorato, and N. Mordant, Phys. Rev. Lett. **122**, 214502 (2019).

[27] A. Costa, A. R. Osborne, D. T. Resio, S. Alessio, E. Chrivi, E. Saggese, K. Bellomo, and C. E. Long, Phys. Rev. Lett. **113**, 108501 (2014).

[28] S. M. Mossman, G. C. Katsimiga, S. I. Mistakidis, A. Romero-Ros, T. M. Bersano, P. Schmelcher, P. G. Kevrekidis, and P. Engels, Communications Physics **7**, 163 (2024).

[29] G. Marcucci, D. Pierangeli, A. J. Agranat, R.-K. Lee, E. DelRe, and C. Conti, Nature Communications **10**, 5090 (2019).

[30] P. Suret, A. Tikan, F. Bonnefoy, F. m. c. Copie, G. Ducrozet, A. Gelash, G. Prabhudesai, G. Michel, A. Cazaubiel, E. Falcon, G. El, and S. Randoux, Phys. Rev. Lett. **125**, 264101 (2020).

[31] L. Dieli, D. Pierangeli, E. DelRe, and C. Conti, Phys. Rev. Lett. **133**, 183801 (2024).

[32] L. Fache, F. Bonnefoy, G. Ducrozet, F. Copie, F. Novkoski, G. Ricard, G. Roberti, E. Falcon, P. Suret, G. El, and S. Randoux, Phys. Rev. E **109**, 034207 (2024).

[33] P. Suret, M. Dufour, G. Roberti, G. El, F. Copie, and S. Randoux, Phys. Rev. Res. **5**, L042002 (2023).

[34] L. Fache, H. Damart, F. Copie, T. Bonnemain, T. Congy, G. Roberti, P. Suret, G. El, and S. Randoux, Phys. Rev. Lett. **134**, 147201 (2025).

[35] D. S. Agafontsev and V. E. Zakharov, Low Temperature Physics **46**, 786 (2020).

[36] D. S. Agafontsev, A. A. Gelash, R. I. Mullyadzhanov, and V. E. Zakharov, Chaos, Solitons & Fractals **166**, 112951 (2023).

[37] F. Coppini, P. G. Grinevich, and P. M. Santini, Phys. Rev. E **101**, 032204 (2020).

[38] A. E. Kraych, P. Suret, G. El, and S. Randoux, Phys. Rev. Lett. **122**, 054101 (2019).

[39] A. E. Kraych, D. Agafontsev, S. Randoux, and P. Suret, Phys. Rev. Lett. **123**, 093902 (2019).

[40] A. Gelash, D. Agafontsev, V. Zakharov, G. El, S. Randoux, and P. Suret, Phys. Rev. Lett. **123**, 234102 (2019).

[41] F. Bonnefoy, A. Tikan, F. Copie, P. Suret, G. Ducrozet, G. Prabhudesai, G. Michel, A. Cazaubiel, E. Falcon, G. El, and S. Randoux, Phys. Rev. Fluids **5**, 034802 (2020).

[42] G. A. El, E. G. Khamis, and A. Tovbis, Nonlinearity **29**, 2798 (2016).

[43] R. Mullyadzhanov and A. Gelash, Phys. Rev. Lett. **126**, 234101 (2021).

[44] see Supplemental Material at [url] for details about the experimental setup, the experimental methodology, the evolution of the DOS in experiment and in numerical simulations.

[45] S. Toenger, T. Godin, C. Billet, F. Dias, M. Erkin-talo, G. Genty, and J. M. Dudley, Sci. Rep. **5** (2015), 10.1038/srep10380.

[46] M. Närhi, B. Wetzel, C. Billet, S. Toenger, T. Sylvestre, J.-M. Merolla, R. Morandotti, F. Dias, G. Genty, and J. M. Dudley, Nature Communications **7**, 13675 (2016).

[47] P. Suret, R. El Koussaifi, A. Tikan, C. Evain, S. Randoux, C. Szwaj, and S. Bielawski, Nature Communications **7**, 13136 (2016).

[48] O. Kimmoun, H. Hsu, H. Branger, M. Li, Y. Chen, C. Kharif, M. Onorato, E. J. Kelleher, B. Kibler, N. Akhmediev, et al., Scientific reports **6**, 28516 (2016).

[49] J. M. Dudley, G. Genty, A. Mussot, A. Chabchoub, and F. Dias, Nature Reviews Physics **1**, 675 (2019).

[50] G. Vanderhaegen, P. Szriftgiser, A. Kudlinski, M. Conforti, A. Armaroli, and A. Mussot, Phys. Rev. A **106**, 033519 (2022).

[51] J. Yang, *Nonlinear Waves in Integrable and Non-integrable Systems*, Mathematical Modeling and Computation (Society for Industrial and Applied Mathematics, 2010).

[52] K. J. Blow, N. J. Doran, and D. Wood, J. Opt. Soc. Am. B **5**, 381 (1988).

[53] T. Okamawari, A. Hasegawa, and Y. Kodama, Phys. Rev. A **51**, 3203 (1995).

[54] J. E. Prilepsky and S. A. Derevyanko, Phys. Rev. E **75**, 036616 (2007).

[55] M. Böhm and F. Mitschke, Phys. Rev. A **76**, 063822 (2007).

[56] A. Hause, C. Mahnke, and F. Mitschke, Phys. Rev. A **98**, 033814 (2018).

[57] I. Bouchoule, B. Doyon, and J. Dubail, SciPost Phys. **9**, 044 (2020).

[58] I. Bouchoule and J. Dubail, Phys. Rev. Lett. **126**, 160603 (2021).
Supplemental material for :
"Perturbed Nonlinear Evolution of Optical Soliton Gases: Growth and Decay in Integrable Turbulence"

Loic Fache,¹ François Copie,¹ Pierre Suret,¹ and Stéphane Randoux¹

¹ Univ. Lille, CNRS, UMR 8523 - PhLAM - Physique des Lasers Atomes et Molécules, F-59 000 Lille, France

The purpose of this Supplemental Material is to provide some details about the experimental setup, the experimental methodology, the evolution of the DOS in experiment and in numerical simulations. All equations, figures, and reference numbers within this document are prepended with “S” to distinguish them from corresponding numbers in the Letter.

CONTENTS

<table><tr><td>Acknowledgments</td><td>5</td></tr><tr><td>References</td><td>5</td></tr><tr><td>I. Detailed description of the experimental setup</td><td>7</td></tr><tr><td>II. Heterodyne measurement of the amplitude and phase of the optical field</td><td>8</td></tr><tr><td>III. Influence of finite bandwidth effects and of detection noise on the measurement of the DOS of the optical soliton gas</td><td>10</td></tr><tr><td>IV. Transition from the Weyl gas to the soliton gas with a circular DOS</td><td>12</td></tr><tr><td>V. Repartition of the spectral mass in the upper complex plane</td><td>14</td></tr></table>

I. DETAILED DESCRIPTION OF THE EXPERIMENTAL SETUP

<figure><img src="image_6.png" /><figcaption>FIG. S5. Detailed representation of the experimental setup. AFG: Arbitrary function generator. AWG: Arbitrary waveform generator. EDFA: Erbium-doped fiber amplifier. PBS: Polarizing beam splitter. HWP: Half-wave plate. AOM: Acousto-optic modulator. EOM: Electro-optic modulator. VOA: Variable optical attenuator.</figcaption></figure>

Fig. S5 shows a detailed representation of our experimental setup, including some technical details about the generation of the optical signals that are not described in the Letter. The light source used for generation of optical
soliton gases (SGs) is a single-frequency continuous-wave (CW) laser diode (APEX-AP3350A) centered at $\lambda = 1550$ nm which delivers an optical power of a few milliwatts. The optical SGs have the initial form of flat-top pulses with durations of 100 ns. These pulses are produced by using a 20-GHz electro-optic modulator (EOM, iXblue MX-LN-20) connected to an arbitrary waveform generator (AWG, Tektronix AWG70001A) having a bandwidth of 15 GHz.

The flat-top pulses are amplified at the Watt-level by using a Erbium-doped fiber amplifier (EDFA, Keopsys CEFA-C-BO-HP). Importantly the amplified spontaneous emission (ASE) from the erbium-doped fiber amplifier (EDFA) introduces a small amount of optical noise to the flat-top pulses. This noise drives their transformation into fully randomized soliton gases as they propagate nonlinearly within the recirculating fiber loop.

The mean power of the light signal at the output of the EDFA is regulated using bulk optics, which include a polarization beam splitter (PBS) and a half-wave plate (HWP). Following this power adjustment, the light signal is directed to an acousto-optic modulator (AOM, VSF MT110-IR25-Fio). The AOM functions as an optical gate, remaining open for approximately 500 ns to allow the injection of the light signal into the fiber loop. It then remains closed for significantly longer periods, typically around 5 ms, during which the signal circulates within the fiber loop. The gating signals for the AOM are generated by a 250-MHz arbitrary function generator (AFG, Tektronix 31252).

The recirculating fiber loop comprises approximately 8 km of single-mode fiber (SMF), configured as a closed loop using a 90/10 fiber coupler. The SMF, manufactured by Draka-Prysmian, has a measured second-order dispersion coefficient of $\beta_2 = -22$ ps$^2$ km$^{-1}$ and an estimated Kerr coefficient $\gamma = 1.3$ km$^{-1}$W$^{-1}$ at the working wavelength of 1550 nm. The 90/10 coupler is arranged so that 90% of the intra-loop optical power is recirculated. At each round trip, 10% of the circulating power is extracted and directed to a 50/50 fiber coupler, where it is mixed with a detuned CW reference laser for heterodyne measurement. The output of the fiber coupler is then routed to a photodetector (Finisar XPDV2120R), which is coupled to a sampling oscilloscope (Teledyne LeCroy LabMaster 10-65Zi, 160 GSa/s). This setup provides an overall detection bandwidth of 32 GHz.

The experimental data comprise a series of temporal sequences, each representing 500 ns per round trip, triggered via the second channel of the AFG. These sequences are processed numerically to construct single-shot space-time diagrams that reveal the wavefield dynamics. The evolution of the mean power within the loop is carefully controlled through Raman amplification via a counter-propagating pump laser at 1450 nm. This 1450 nm laser is a commercial Raman fiber laser (IPG FiberTech) capable of delivering several watts of optical power. It is coupled into and out of the loop through wavelength division multiplexers (WDMs). For the experiments, the power at 1450 nm is attenuated to approximately 200 mW using a 90/10 fiber coupler (not shown in Fig. S5). The fine control of the pump power injected within the recirculating fiber loop is made using an additional AOM together with a variable optical attenuator (VOA). Note that a slow photodetector (Thorlabs DET10C/M) connected to a 100-MHz oscilloscope (Agilent 54624A) is used to monitor the slow (ms scale) evolution of the power of the soliton gas within the fiber loop.

The heterodyne measurement of the optical signal's amplitude and phase is performed by mixing it with a reference signal generated by a single-frequency CW laser diode (APEX-AP3350A) operating at $\lambda + \delta\lambda = 1550.361$ nm, as detailed further in Sec. II.

## II. HETERODYNE MEASUREMENT OF THE AMPLITUDE AND PHASE OF THE OPTICAL FIELD

As shown in Fig. S1, the fast PD placed at the output of the recirculating fiber loop detects the beating between the optical signal at 1550 nm and a frequency-detuned reference laser at a wavelength of 1550.361 nm (see Fig. S6(a)). The reference laser (APEX-AP3350A) is a single-frequency laser diode with an output power of a few mW and an intrinsic linewidth of 300 kHz (constructor datasheet). The wavelength difference (0.361 nm) between the two lasers corresponds to a frequency difference $\Delta f_0 \sim 45$ GHz at which the beating between the signal laser and the reference laser is detected by the fast PD (see Fig. S6(b)). The voltage detected by the PD is given by:

$$V(t) = k \left| A_{\text{sig}}(t) e^{i\phi(t)} + A_{\text{Ref}} e^{i(2\pi\Delta f_0 t + \phi_{\text{Ref}})} \right|^2 \quad (S1)$$

where k is a real constant whose value is associated with the efficiency of the PD. $A_{\text{sig}}(t)$ (resp. $\phi(t)$) represents the slowly-varying amplitude (resp. the phase) of the signal field. $A_{\text{Ref}}$ (resp. $\phi_{\text{Ref}}$) represents the amplitude (resp. the phase) of the reference field. Both $A_{\text{Ref}}$ and $\phi_{\text{Ref}}$ are assumed to remain constant during the measurement period,
<figure><img src="image_7.png" /><figcaption>FIG. S6. Demodulation scheme of the detected optical signal. (a) Typical experimental signal measured by the photodiode and the oscilloscope at $z = 400$ km over 2 ns. The red dashed line represents the mean value of the detected signal. (b) Fourier power spectrum of the signal plotted in (a). The red line represents the modulus of the transfer function of the passband Butterworth filter used to filter the frequency components centered around 45 GHz. (c) Black line: amplitude of the demodulated signal $A_{sig}(t)$. Red line: phase $\phi(t)$ of the optical signal. (d) Fourier power spectrum associated to the demodulated signal $A_{sig}(t) \exp(i\phi(t))$ shown in (c).</figcaption></figure>

which spans 100 ns—the duration of the flat-top pulses corresponding to the extent of the optical soliton gas in the experiment.

Expanding Eq. (S1), we get:

$$V(t) = k \left( A_{\text{sig}}^2(t) + A_{\text{Ref}}^2 + A_{\text{sig}}(t)A_{\text{Ref}}e^{i(2\pi\Delta f_0 t + \phi_{\text{Ref}} - \phi(t))} + A_{\text{sig}}(t)A_{\text{Ref}}e^{-i(2\pi\Delta f_0 t + \phi_{\text{Ref}} - \phi(t))} \right) \quad (S2)$$

This expression shows that the Fourier spectrum of the detected voltage is composed of three components: one DC component and two side-bands at $\pm\Delta f_0 \sim \pm 45$ GHz, as shown in Fig. S6(b). Considering that $\Delta f_0$ is sufficiently large that the overlap between the 3 spectral components can be neglected, we use a highly selective frequency filter to eliminate the DC component and the frequency components around -45 GHz. The filtering operation is purely numerical; it is applied at the post-processing stage on the signal recorded by the fast oscilloscope. We have used a causal 15th-order Butterworth filter with a bandwidth of $\Omega = 40$ GHz (full-width at half maximum) centered at $\Delta f_0 = +45$ GHz. The modulus of the transfer function of this filter is plotted in red line in Fig. S6(b).

The electrical signal $V_F(t)$ at the output of the filter reads

$$V_F(t) = k \left( A_{\text{sig}}(t) A_{\text{Ref}} e^{i(2\pi \Delta f_0 t + \phi_{\text{Ref}} - \phi(t))} \right) \quad (S3)$$

The next step in the signal processing is a demodulation operation in which the filtered signal $V_F(t)$ is multiplied by $e^{-2i\pi\Delta f_0 t}$, which results into a demodulated signal $V_{\text{dem}}(t)$ that reads:

$$V_{\text{dem}}(t) = k \left( A_{\text{sig}}(t) A_{\text{Ref}} e^{i(\phi_{\text{Ref}} - \phi(t))} \right) \quad (S4)$$

Removing the constant phase $\phi_{\text{Ref}}$ of the reference laser, the phase of the optical signal is extracted numerically by taking the argument of $V_{\text{dem}}(t)$ in Eq. (S4):

$$\phi(t) = -\operatorname{Arg}(V_{\text{dem}}(t)) \quad (S5)$$

The last step in the measurement procedure consists in converting the voltage $V_{\text{dem}}(t)$ into the power $A_{\text{sig}}^2(t)$ of the signal field. Following the experimental approach already used in ref. [33], we take advantage of the well-known process of noise-induced modulation instability (MI) of plane waves in the focusing propagation regime [39, 45–47, 49, 50]. After some propagation distance within the fiber loop, the flat-top square pulses perturbed by the optical noise of the EDFA are destabilized, which results in the emergence of a random train of coherent structures. In the Fourier spectrum, this phenomenon appears as the growth of two spectral sidebands symmetrically shifted by $\pm\Omega_{\text{max}}$ relative to the carrier frequency of the plane wave (see Fig. S6(d)). The frequency $\Omega_{\text{max}}$ is related to the power $P_0$ of the plane wave (or of the flat-top pulse) through the relation $\Omega_{\text{max}} = (2\gamma P_0/|\beta_2|)^{1/2}$.

In practice, measuring $\Omega_{\text{max}}$ is more precise when performed in the time domain by simply counting the peaks of the structures that emerge from modulation instability (MI) over a long observation window of several tens of nanoseconds. Using this approach, we determine the period of MI as $T_{\text{MI}} = 2\pi/\Omega_{\text{max}} = 190.1$ ps. This corresponds to an optical power of $P_0 = (2\pi/T_{\text{MI}})^2|\beta_2|/(2\gamma) = 9.24$ mW. However, this measurement, taken at $z = 400$ km (after
the full development of MI), needs to be corrected to account for the optical losses in the fiber loop. Since the power in the loop decays exponentially, the corrected value is $P_0 \approx 11.5$ mW. This calibrated value is crucial, as it enables the conversion of the voltage measured by the fast PD into an equivalent optical power.

### III. INFLUENCE OF FINITE BANDWIDTH EFFECTS AND OF DETECTION NOISE ON THE MEASUREMENT OF THE DOS OF THE OPTICAL SOLITON GAS

In this section, we employ numerical simulations to demonstrate how noise and finite bandwidth effects, intrinsic to our experimental detection system, introduce distortions in the measured amplitude and phase of the optical signal. These distortions, in turn, impact the measured IST spectrum and the density of states (DOS) of the optical soliton gas.

#### Influence of the heterodyne detection scheme

<figure><img src="image_8.png" /></figure>

FIG. S7. (a) The black line represents the power fluctuations of the optical SG obtained from numerical simulation of NLSE with damping at $z = 1920$ km (Eq. (1) of the Letter). The red dashed line represents power fluctuations that are reconstructed using the heterodyne scheme described in Sec. II. (b) Fourier power spectrum of the optical SG. (c) IST spectrum of the optical SG. (d) DOS computed from 10 realizations of the SG. (e) Same as in (a) but for the optical phase of the SG. (f) Fourier power spectrum of the heterodyne signal at 45 GHz. (g) IST spectrum of the optical SG reconstructed from heterodyne measurement. (h) DOS of the optical SG reconstructed from the heterodyne measurement.

As a preliminary step, we first demonstrate, using Fig. S7, that the heterodyne detection scheme described in Sec. II preserves the accuracy of the measurements of the IST spectrum and the density of states (DOS) of the soliton gas.

Fig. S7(a) represents numerical simulations of the damped NLSE (Eq. (1) of the Letter) showing typical power fluctuations of the optical SG at $z = 1920$ km. Fig. S7(b) and S7(c) represent the associated Fourier and IST spectra. Fig. S7(d) represents the DOS that has been calculated using 10 numerical realizations of the optical SG.

Fig. S7(g) and S7(h) in the bottom row of Fig. S7 demonstrate that the IST spectrum and the DOS remain unaffected by the heterodyne detection scheme described in Sec. II, which essentially involves shifting the Fourier spectrum around 45 GHz, as illustrated in Fig. S7(f). This invariance is also clearly evident in Fig. S7(a) and S7(e), where the red and black lines show that the amplitude and phase reconstructed from the heterodyne scheme closely match the numerical counterparts.
## Influence of the finite detection bandwidth

<figure><img src="image_9.png" /><figcaption>FIG. S8. (a) The black (resp. blue) line represents the power (resp. phase) fluctuations of the optical SG obtained from numerical simulation of NLSE with damping at $z = 1920$ km (Eq. (1) of the Letter). The red (resp. orange) dashed line represents power (resp. phase) fluctuations that are reconstructed using the heterodyne scheme including finite bandwidth effects. (b) Fourier power spectrum of the heterodyne signal. The red dashed line represents the modulus of the transfer function of the 32 GHz filter describing finite bandwidth effects. (c) IST spectrum of the optical SG reconstructed from heterodyne measurement with finite bandwidth effects. (d) DOS of the optical SG reconstructed from the heterodyne measurement with finite bandwidth effects.</figcaption></figure>

We now demonstrate that the finite detection bandwidth of our experimental system introduces a slight distortion in the IST spectrum, which in turn affects the measured DOS of the optical soliton gas.

As described in Sec. I, the power fluctuations of the SG are recorded using a fast photodiode (Finisar XPDV2120R) with a nominal bandwidth of 50 GHz, connected to a fast oscilloscope (LeCroy Labmaster 10-65ZI) with a nominal bandwidth of 65 GHz and a sampling rate of 160 GSa/s. Together, these devices behave as a first-order electrical filter having a transfer function given by

$$H(f) = \frac{1}{1 + i \frac{f}{f_c}}, \qquad (S6)$$

with $f_c = 32$ GHz.

The modulus of this transfer function is represented by the red dashed line in Fig. S8(b), illustrating the effect of the low-pass filter on the frequency-shifted signal at 45 GHz. As shown in Fig. S8(a), the power and phase of the signal reconstructed from the filtered heterodyne signal closely match the numerical profiles that would result without any finite bandwidth effects. While these differences due to the finite bandwidth appear minimal, the computed IST spectrum in Fig. S8(c) reveals an asymmetry, which becomes even clearer in the DOS of the optical SG, as depicted in Fig. S8(d).

## Influence of detection noise

In the last step, we consider the influence of the detection noise of the PD on the measured IST spectra and DOS.

Fig. S9(b) shows that we have modified the Fourier spectrum of the heterodyne signal by adding some noise having a level in the Fourier domain that is similar to the experimental detection noise, compare Fig. S9(b) with Fig. S9(f) and also with Fig. S8(b). As a result of the incorporation of this detection noise in our numerical simulations, we observe that the power and phase of the signal reconstructed from the noisy filtered heterodyne signal now deviate from their nominal numerical profiles, see Fig. S9(a) and Fig. S9(e).

By comparing Fig.S9(d) and Fig.S9(h) with Fig. S8(d), we conclude that the detection noise of the PD contributes to the blurred appearance of the measured DOS.

In summary, the symmetric semi-circular shape of the DOS predicted by numerical simulations of an optical SG in the presence of damping is not perfectly replicated in the experiment, primarily due to noise and finite bandwidth effects. The asymmetric shape of the DOS measured in the experiment is primarily due to the finite bandwidth of the detection system. Detection noise contributes to the blurred appearance of the measured DOS.
<figure><img src="image_10.png" /></figure>

FIG. S9. (a) The black line represents the power fluctuations of the optical SG obtained from numerical simulation of NLSE with damping at $z = 1920$ km (Eq. (1) of the Letter). The red dashed line represents power fluctuations that are reconstructed using the heterodyne scheme including finite bandwidth effects and detection noise. (b) Fourier power spectrum of the heterodyne signal perturbated by the added detection noise. The red dashed line represents the modulus of the transfer function of the 32 GHz filter describing finite bandwidth effects. (c) IST spectrum of the optical SG reconstructed from heterodyne measurement with noise and finite bandwith effects. (d) DOS of the optical SG reconstructed from the heterodyne measurement with noise and finite bandwith effects. (e) The black line represents the phase fluctuations of the optical SG obtained from numerical simulation of NLSE with damping at $z = 1920$ km. The red dashed line represents phase fluctuations that are reconstructed using the heterodyne scheme including finite bandwidth effects and detection noise. (f) Fourier power spectrum of the heterodyne signal recorded in the experiment at $z = 1920$ km. (g) Associated IST spectrum of the experimental field. (h) DOS experimentally measured at $z = 1920$ km.

## IV. TRANSITION FROM THE WEYL GAS TO THE SOLITON GAS WITH A CIRCULAR DOS

In this section, we provide some additional experimental results on the transition from the initial Weyl gas to the soliton gas with a circular DOS. Also, we use numerical simulations to clearly show that the DOS of the perturbed optical SG cannot achieve a steady state because of the continuous change in power (or mass), which fundamentally governs its dissipative evolution.

<figure><img src="image_11.png" /></figure>

FIG. S10. Experiments. DOS of the optical SG at various propagation distances from $z = 8$ km (a) to $z = 3200$ km (g).

Figures S10(a) and S10(b) illustrate that during the initial evolution from $z = 0$ km to $z = 400$ km, the imaginary parts of the discrete eigenvalues undergo a slight decay. At this early stage, no significant changes in the geometrical shape of the DOS are observed. However, a gradual broadening of the DOS emerges around $z = 800$ km, see Fig. S10(c). Between $z = 800$ km and $z = 3200$ km, the imaginary parts of the eigenvalues progressively decay, while the width of the distribution steadily expand, see Fig. S10(c)–(g).

The spectral changes observed in the experiment are more clearly illustrated through numerical simulations of the focusing 1D-NLSE (Eq. (1) of the manuscript) with a smaller dissipation parameter, $\alpha_{eff}$. Our numerical simulations have been performed by using a pseudo-spectral method working with a step-adaptative algorithm permitting to reach
a specified level of numerical accuracy. The numerical simulations are performed by using a box of size $\Delta T = 180$ ns that has been discretized by using 16384 points. The initial condition is under the form of a super-Gaussian pulse with full width at half maximum $T_0 \approx 100$ ns, mean power $P_0 = 11.5$ mW, and an exponent $p = 200$:

$$A(T, z = 0) = \sqrt{P_0} \exp \left[ - \left( \frac{T^2}{2\sigma^2} \right)^p \right], \quad \text{with} \quad \sigma = \frac{T_0}{2\sqrt{2\ln 2}}.$$

This envelope is perturbed by the addition of a small noise generated in the spectral domain from a Gaussian-shaped spectrum with Fourier components having random phases and amplitudes.

Figure S11(a)-(g) presents the simulated evolution of the DOS from $z = 0$ km to $z = 23200$ km for a value of the dissipation parameter $\sim 8$ times smaller than the experimental one: $\alpha_{eff} = 7.6 \times 10^{-5}$ km$^{-1}$ (equivalently $\epsilon = 0.005$ instead of $\epsilon = 0.038$ in Eq. (2) of the Letter). With a constant loss rate, the optical SG never reaches a steady state. Instead, the circular support of the DOS undergoes continuous homothetic evolution, characterized by a progressive decay in the imaginary parts of the discrete eigenvalues and a corresponding increase in the real parts.

<figure><img src="image_12.png" /><figcaption>FIG. S11. (a)-(g): Numerical simulations of the focusing 1D-NLSE (Eq. (1) of the Letter) with a value of the dissipation parameter $\sim 8$ times smaller than the experimental one: $\alpha_{eff} = 7.6 \times 10^{-5}$ km$^{-1}$ (equivalently $\epsilon = 0.005$ instead of $\epsilon = 0.038$ in Eq. (2) of the Letter). The numerical simulations show that the circular support of the DOS undergoes continuous homothetic evolution, characterized by a progressive decay in the imaginary parts of the discrete eigenvalues and a corresponding increase in the real parts. (h)-(n) Same spectral distributions as in (a)-(g), but after applying the exponential rescaling of the spectral parameters described in the text ($v = \eta \exp(\alpha_{eff}z)$, $u = \xi \exp(p\alpha_{eff}z)$, $F(u, v) = f(\xi, \eta) \exp(q\alpha_{eff}z)$). The bottom row illustrates that the semi-elliptic distributions shown in the top row can be approximately mapped onto a common semicircular distribution at large enough propagation distance.</figcaption></figure>

The semi-elliptic shape of the DOS shown in Fig. S11(g) at $z = 23200$ km can be rescaled to a semi-circular form by transforming the spectral parameter from $\eta$ to $v = \eta \exp(\alpha_{eff}z)$, as illustrated in Fig. S11(n). This scaling is physically intuitive, reflecting the simple idea that the soliton amplitudes -linked to the spectral parameter $\eta$ - decay exponentially with the propagation distance $z$. Applying the same scaling to the DOS at all other intermediate propagation distances eliminates the vertical decay of the spectral distribution along the $\eta$-axis.

Building on the previous observation, it is natural to introduce a corresponding exponential scaling for the spectral parameter $\xi$ and for the DOS $f(\xi, \eta)$:

$$u = \xi \exp(p\alpha_{eff}z), \quad F(u, v) = f(\xi, \eta) \exp(q\alpha_{eff}z)$$

where $p$ and $q$ are real parameters chosen such that the DOS at all intermediate propagation distances best fit the DOS shown in Fig. S11(g) at the longest propagation distance, $z = 23200$ km. The corresponding results, plotted in the second row of Fig. S11, show that all DOS beyond $z = 5800$ km can be mapped onto nearly the same spectral curve computed at the longest distance of $z = 23200$ km. Further studies are needed to determine to what extent this approximate self-similar evolution at long propagation distances may be connected to the existence of adiabatic invariants.
The spectral evolution presented in Fig. S11 is accompanied by the space-time evolution depicted in Fig. S12, where a gradual expansion of the optical SG is observed. The emergence of curved spectral contours in the DOS, driven by the power (or mass) drop as evidenced in Figure S12(d), results in the creation of soliton velocities visible at the edges of the square pulse, around $t = \pm 50$ ns, see Fig. S12(a) and S12(c).

<figure><img src="image_13.png" /><figcaption>FIG. S12. Numerical simulations of the focusing 1D-NLSE (Eq. (1) of the Letter) with a value of the dissipation parameter $\sim 8$ times smaller than the experimental one: $\alpha_{eff} = 7.6 \times 10^{-5} \text{ km}^{-1}$ (equivalently $\epsilon = 0.005$ instead of $\epsilon = 0.038$ in Eq. (2) of the Letter). The numerical simulations illustrate the spatiotemporal evolution of the optical SG across three distinct time windows presented in panels (a), (b), and (c). Panels (a) and (c) reveal the expansion of the SG due to dissipation, particularly at the edges of the pulse. Panel (d) depicts the longitudinal evolution of the power of the optical SG.</figcaption></figure>

The results illustrated in Figures S10, S11, and S12 show that the DOS of the perturbed optical SG cannot achieve a steady state. This is due to the continuous change in power (or mass), which fundamentally governs its dissipative evolution.

V. REPARTITION OF THE SPECTRAL MASS IN THE UPPER COMPLEX PLANE

Upon careful examination of Fig. S10 and Fig. S11, it is evident that the DOS measured in the scenario where a circular SG is formed due to loss is not uniformly circular. Instead, there exists a monochromatic component around $\xi = 0$. This component decays concurrently with the circular distribution over the course of the SG evolution, as illustrated in Fig. S10. We note that numerical simulations conducted with a dissipation rate lower than that observed experimentally reveal a reduction in the strength of this monochromatic component, as illustrated in Fig. S11.

In order to quantify this observation more precisely, we measured the ratio $R(z)$ between the mass of the monochromatic component around $\xi = 0$ and the total mass of the soliton gas at specific longitudinal positions $z$. Mathematically, the ratio $R(z)$ is defined as:

$$R(z) = \frac{\int_{-\alpha}^{+\alpha} d\xi \int_{\eta_0(z)-\delta}^{\eta_0(z)+\delta} d\eta \eta f(\xi, \eta, z)}{\int_{-\infty}^{+\infty} d\xi \int_{0}^{+\infty} d\eta \eta f(\xi, \eta, z)} \quad (S7)$$

where $\alpha$ and $\delta$ are small numbers ($\alpha = 0.05$, $\delta = 0.05$) that delineate a small area surrounding the monochromatic component. $\eta_0(z)$ represents the central position along the vertical imaginary axis of the cluster of eigenvalues associated with the monochromatic component of the SG. The areas taken for the measurement of the mass of the monochromatic component are depicted with green lines in Fig. S13, where the numerical values of $R(z)$ are also indicated in green. It can be seen that the mass of the monochromatic component decreases from 33% at $z = 0$ km to 7.5% at $z = 23200$ km.

Similar questions regarding the distribution of spectral mass in the upper complex plane can also be explored through numerical simulations for the scenario where a bichromatic SG is formed due to amplification. Figure S14 illustrates the evolution of the DOS between $z = 0$ km and $z = 12800$ km, a distance twice as large as that achieved
<figure><img src="image_14.png" /><figcaption>FIG. S13. Numerical simulations of the focusing 1D-NLSE (Eq. (1) of the Letter) with a value of the dissipation parameter ~8 times smaller than the experimental one : $\alpha_{eff} = 7.6 \times 10^{-5} \text{ km}^{-1}$ (equivalently $\epsilon = 0.005$ in Eq. (2) of the Letter). The small areas depicted with green lines indicate the regions where the calculation of the mass of the monochromatic component is made. The fraction of the mass $R(z)$ of this monochromatic component is indicated in green. It decreases from 33% at $z = 0$ km to 7.5% at $z = 23200$ km.</figcaption></figure>

in our experiment. The formation of two distinct clusters of eigenvalues at significant evolution distances indicates that the resulting SG exhibits characteristics of a bichromatic SG. It is important to note, however, that the mass associated with the second cluster of eigenvalues, which have small imaginary parts, is very minimal. Specifically, it does not exceed a few percent of the total mass of the soliton gas, as shown in Figure S14.

<figure><img src="image_15.png" /><figcaption>FIG. S14. Numerical simulations of the focusing 1D-NLSE (Eq. (1) of the Letter) with a value of the dissipation parameter identical the experimental one : $\alpha_{eff} = -1.1 \times 10^{-4} \text{ km}^{-1}$ (equivalently $\epsilon = -0.042$ in Eq. (2) of the Letter). The small areas depicted with green lines indicate the regions where the calculation of the mass of the monochromatic component is made. The fraction of the mass $R(z)$ of this monochromatic component is indicated in green. It decreases from 7.6% at $z = 0$ km to 2.4% at $z = 12800$ km.</figcaption></figure>