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

HAL
open science

# Manipulation of Strongly Interacting Solitons in Optical Fiber Experiments

Alexandre Mucci, Pierre Suret, François Copie, Stephane Randoux, Rustam

Mullyadzhanov, Andrey Gelash

► To cite this version:

Alexandre Mucci, Pierre Suret, François Copie, Stephane Randoux, Rustam Mullyadzhanov, et al. Manipulation of Strongly Interacting Solitons in Optical Fiber Experiments. Physical Review Letters, 2025, 134 (19), pp.193804. 10.1103/PhysRevLett.134.193804. hal-05099197

HAL Id: hal-05099197

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

Submitted on 5 Jun 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Manipulation of strongly interacting solitons in optical fiber experiments

Alexandre Mucci¹, Pierre Suret¹, François Copie¹, Stephane

Randoux¹, Rustam Mullyadzhanov²,³, and Andrey Gelash⁴*

¹University Lille, CNRS, UMR 8523 - PhLAM - Physique des Lasers Atomes et Molécules, F-59000 Lille, France

²Novosibirsk State University, Novosibirsk 630090, Russia

³Institute of Thermophysics SB RAS, Novosibirsk 630090, Russia and

⁴Laboratoire Interdisciplinaire Carnot de Bourgogne (ICB),
UMR 6303 CNRS-Université Bourgogne Franche-Comté, 21078 Dijon, France

The model underlying physics of guiding light in single-mode fibers – the one-dimensional nonlinear Schrödinger equation (NLSE), reveals a remarkable balance of the fiber dispersion and nonlinearity, leading to the existence of optical solitons. With the Inverse Scattering Transform (IST) method and its perturbation theory extension one can go beyond single-soliton physics and investigate nonlinear dynamics of complex optical pulses driven by soliton interactions. Here, advancing the IST perturbation theory approach, we introduce the eigenvalue response functions, which provide an intuitively clear way to manipulate individual characteristics of solitons even in the case of their entire overlapping, i.e. very strong interactions. The response functions reveal the spatial sensitivity of the multi-soliton pulse concerning its instantaneous perturbations, allowing one to manipulate the velocities and amplitudes of each soliton. Our experimental setup, based on a recirculating optical fiber loop system and homodyne measurement, enables observation of long-distance NLSE dynamics and complete characterization of optical pulses containing several solitons. Adding a localized phase perturbation on a box-shaped wavefield, we change individual soliton characteristics and can detach solitons selectively from the whole pulse. The detached solitons exhibit velocities very close to the theoretical predictions, thereby demonstrating the efficiency and robustness of the response functions approach.

Experimental observation of optical solitons in single-mode fibers reported in 1980 [1] and the following demonstration of their stable propagation over thousands of kilometers maintained by a loss compensation [2] opened avenues into fundamental studies and engineering applications in the field of nonlinear optics [3, 4]. The model underlying physics of guiding light in fibers – the one-dimensional nonlinear Schrödinger equation (NLSE), reveals a remarkable balance of the fiber dispersion and nonlinearity, leading to the effect of soliton shape preservation [3, 5–7]. Manipulation of nonlinear light pulses represents a long-standing challenge in fiber optics. Nowadays, thanks to modern ultra-fast lasers, wave-shapers, and electro-optic modulators, the reproduction and management of individual solitons in optical fibers become a routine procedure used e.g., in supercontinuum generation [8, 9] or as a tool for exploring various optical phenomena [10–12].

To go beyond single-soliton physics, one can rely on a powerful theoretical technique – the Inverse Scattering Transform (IST) method based on the integrable nature of the NLSE [13–15]. The IST realizes a mapping of nonlinear wavefield dynamics onto linear evolution of the IST spectrum, more rigorously referred to as scattering data. In this sense, the IST procedure and scattering data are nonlinear analogs of the conventional Fourier transform and the Fourier spectra respectively [14, 16]. Considering the focusing NLSE, the IST spectrum consists of discrete and continuous components corresponding to solitons and nonlinear dispersive waves. One can find the IST spectrum by solving an auxiliary

linear Zakharov-Shabat (ZS) eigenvalue system and vice versa – obtain the wave field from known IST spectrum. For example, soliton interactions can be described using exact multi-soliton formulas, parameterized by the IST soliton eigenvalues. In general case, the evolution of initial conditions containing non-zero continuous spectrum parts can be expressed via the integral Gelfand-Levitan-Marchenko (GLM) equations. In addition, the well-developed IST perturbation theory allows one to approximate the dynamics of perturbed solitons and generalize the IST approach to nearly integrable models, e.g., the NLSE with third-order dispersion, gain-loss or Raman scattering terms [17–22].

In this Letter, we consider theoretically and study experimentally nonlinear propagation of perturbed box-shaped light pulses in an optical single-mode fiber system used in the anomalous dispersion regime. The pulses consist of several solitons and a portion of nonlinear dispersive waves, which complicates the straightforward IST analysis because, in this case, the GLM equations can only be solved asymptotically [14, 15, 23]. Meanwhile, numerical simulations, always available for the NLSE model, do not give a general description of the problem. To cope with the analysis difficulties we advance the IST perturbation theory and introduce the eigenvalue response functions (RFs) which provide an intuitively clear way to manipulate strongly interacting solitons in nonlinear pulses. The RFs are computed individually for each of the solitons via the corresponding eigenfunction solutions of the ZS system and reveal the spatial sensitivity of the pulse concerning wavefield perturbations.
Being linked to the integral describing change of soliton eigenvalues in response to the perturbation, the RFs allow one to manipulate individual soliton velocities and amplitudes perturbing the pulse in a controllable way.

As a proof of concept, we design experiments for a recirculating optical fiber loop system with Raman amplifier and homodyne interferometric tools, enabling the recording of long-distance NLSE spatiotemporal dynamics together with the full characterization of complex optical pulses (phase and amplitude). Adding a localized phase perturbation of the light pulse we manage to change individual soliton characteristics and even detach a desired soliton or a soliton complex from the whole pulse. The detached solitons being measured in the fiber loop system exhibit velocities very close to the theoretical ones, thereby demonstrating the efficiency and robustness of the proposed approach.

The model of light propagation in our experimental setup is the NLSE with effective losses term [11, 24, 25],

$$i \frac{\partial \Psi}{\partial Z} = \frac{\beta_2}{2} \frac{\partial^2 \Psi}{\partial T^2} - \gamma |\Psi|^2 \Psi - i \frac{\alpha}{2} \Psi, \quad (1)$$

where $\Psi(Z, T)$ is the complex envelope of the electric field $E(Z, t) = \Psi(Z, t-z/v_g) \exp(i(k_0Z - \omega_0t) + \text{c.c.}$, that slowly varies in space $Z$ and the retarded time $T = t - Z/v_g$, measured in the frame propagating at the group velocity $v_g = v_g(\omega_0)$. The laboratory frame time is $t$, while $\omega_0$ and $k_0$ represent the carrier wave frequency and wavenumber. At our working wavelength 1550 nm, the group velocity dispersion coefficient is $\beta_2 = -22 \text{ ps}^2/\text{km}$, the nonlinear Kerr coefficient is $\gamma = 1.3 \text{ W}^{-1}\text{km}^{-1}$ and the recirculating fiber loop losses $\alpha = 4 \times 10^{-7} \text{ m}^{-1}$.

Our theoretical treatment is based on the dimensionless form of the NLSE,

$$i\psi_z + \frac{1}{2}\psi_{\tau\tau} + |\psi|^2\psi = 0, \quad (2)$$

which is obtained from (1) by neglecting the losses and using the rescaled variables $\psi = \Psi/\sqrt{P_0}$, $\tau = T\sqrt{\gamma P_0/|\beta_2|}$, $z = Z\gamma P_0$, where $P_0$ is a typical value of the optical power $|\Psi|^2$. The IST analysis of the NLSE wavefields relies on the linear auxiliary Zakharov-Shabat (ZS) system for a vector wave function $\Phi = (\phi_1, \phi_2)^T$:

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

where $\lambda = \xi + i\eta$ is the z-independent complex spectral parameter and T means transposition. The NLSE wavefield $\psi(\tau)$, playing the role of potential for the system (3), in general contains a finite number of solitons N and a portion of incoherent nonlinear dispersive waves [14, 15]. Each soliton in the wavefield corresponds to a discrete eigenvalue $\lambda_n = \xi_n + i\eta_n$ of the operator $\hat{\mathcal{L}}$, where $n = 1, \dots, N$ and $\eta_n > 0$. The eigenvalue defines asymptotic soliton amplitude A and velocity V as,

$$A_n = 2\eta_n, \quad V_n = 2\xi_n. \quad (4)$$

<figure><img src="image_3.png" /><figcaption>FIG. 1: Eigenvalue response functions (RFs) of the box-shaped pulse with $L = 8.85$ and manipulation of its solitons in numerical simulation. (a) The initial pulse (dashed gray) and $RF_n^{im}$. The pulse contains three solitons with the eigenvalues indicated on the inset and colored following the colors of the RFs. (b) Reference evolution of the unperturbed pulse in numerical simulation. (c) Evolution of the pulse having local imaginary perturbation; see red and green solid lines for power and phase modulation of the initial field in (a). Dashed lines in (c) show soliton trajectories predicted by theoretical expressions (8). The gray strip in (a) shows the area of the RF integration, which contributes to the soliton eigenvalue responses.</figcaption></figure>

When nonlinear dispersive waves are absent, the wavefield evolution obeys exact multi-soliton solutions, for $N = 1$ given by the well-known single soliton formula,

$$\psi_s = A_1 \frac{\exp\{iV_1(\tau - \tau_0) + i(A_1^2 + V_1^2)z/2 + i\theta_0\}}{\cosh A_1(\tau - V_1z - \tau_0)}, \quad (5)$$

where $\tau_0$ and $\theta_0$ are soliton position and phase.

In the presence of nonlinear dispersive waves exact analytical solutions of the NLSE (2) are not available; however analysis of soliton eigenvalues allows one to access the main features of the nonlinear wavefield behavior. Consider a box-shaped real pulse $\psi_{\perp}$ of amplitude A and width L, so that $\psi_{\perp}(\tau) = A$ for $|\tau| < L/2$ and $\psi_{\perp}(\tau) = 0$ otherwise, see Fig. 1(a). The pulse contains $N = \text{Integer}[1/2 + AL/\pi]$ zero velocity solitons with eigenvalues aligned on the imaginary axis: $\{\lambda_n\} = i\{\eta_n\}$ for $n = 1, \dots, N$, where $\eta_n$ obey the following transcendental equation [26]:

$$\tan(\chi_n L) = -\chi_n/\eta_n, \quad \chi_n = \sqrt{A^2 - \eta_n^2}. \quad (6)$$

For an arbitrary-shaped pulse, the eigenvalues can be found numerically using various available IST algorithms [22, 27–30].

The spatiotemporal dynamics $\psi_{\perp}(z, \tau)$ features the oscillating interactions of N solitons and the radiation,
<figure><img src="image_4.png" /><figcaption>FIG. 2: (a) Schematic experimental setup. (b) Measured power decay rate of the circulating light (red) and its exponential fit (black) with $\alpha = 4 \times 10^{-7} \text{ m}^{-1}$ in normalized units. (c, e) Spatiotemporal diagrams of the perturbed pulse propagation over 1200 km, with soliton trajectories predicted by theoretical expressions (8). (d, f) Initial pulse power (red), its phase profile (green) and the approximating box-shaped potential used to compute the RFs. (g, h) Discrete IST spectrum: direct numerical computation (circles) and theoretical eigenvalues computed from expressions (8) (dots). Panels (c, d, g) demonstrate the experimental data, while (e, f, h) show corresponding numerical simulations of the full model [Eq. (1)].</figcaption></figure>

see Fig. 1(b) obtained by numerical integration of the NLSE (2). An instant disturbance of the pulse $\psi(\tau) \to \psi(\tau) + \delta\psi(\tau)$ leads to a change of the initial soliton eigenvalues, that for relatively weak $\delta\psi(\tau)$ can be captured by the first-order IST perturbation theory [21, 31, 32]:

$$ \delta\lambda = \frac{\langle \Phi^\dagger, \delta\hat{\mathcal{L}}\Phi \rangle}{\langle \Phi^\dagger, \Phi \rangle}, \quad \delta\hat{\mathcal{L}} = -i \begin{pmatrix} 0 & \delta\psi \\ \delta\psi^* & 0 \end{pmatrix}. \quad (7) $$

Here $\Phi(\tau)$ represent solution of the ZS system with potential $\psi(\tau)$, $\Phi^\dagger = (\phi_2^*, \phi_1^*)^\text{T}$ and the scalar product $\langle f, g \rangle = \int_{-\infty}^{\infty} f^* g d\tau$.

For the rectangular potential $\psi(\tau)$ and its perturbation $\delta\psi(\tau)$ represented by a sum of Fourier harmonics, the integrals (7) were evaluated recently in [33]. Here, to simplify the analysis of the problem (7), we introduce the eigenvalue response functions (RFs) for arbitrary real-valued $\delta\psi^\text{re}(\tau) = f(\tau)$ and imaginary-valued $\delta\psi^\text{im}(\tau) = if(\tau)$ pulse perturbations as,

$$ \delta\lambda_n = \int_{-\infty}^{\infty} \text{RF}(\tau)_n^{\text{re}/\text{im}} f(\tau) d\tau, \quad (8) $$

$$ \text{RF}_n^{\text{re}} = -i \frac{\phi_1^2 + \phi_2^2}{\Delta} \Big|_{\lambda=\lambda_n}, \quad \text{RF}_n^{\text{im}} = \frac{-\phi_1^2 + \phi_2^2}{\Delta} \Big|_{\lambda=\lambda_n}, $$

$$ \Delta = \langle \Phi^\dagger, \Phi \rangle = 2 \int_{-\infty}^{\infty} \phi_1 \phi_2 d\tau. $$

Eq. (8) shows that the eigenvalue response to the pulse perturbation represents the integration of $f(\tau)$ with RF($\tau$). The latter plays the role of a weight function that provides an intuitively simple way to manipulate characteristics of individual solitons by managing $f(\tau)$. Remarkably, the functions $\text{RF}_n^{\text{re}/\text{im}}(\tau)$ are always either pure

imaginary or pure real-valued; see Supplemental Material [34], where we provide exact expressions for them. Thus, the real-valued perturbations of the box pulse can change only the imaginary part of $\lambda_n$, i.e., affect soliton amplitudes $A_n$, which is in line with general properties of the ZS system [35]. Meanwhile, imaginary-valued perturbations act vice versa and change only soliton velocities $V_n$.

Here, we focus on the manipulation of soliton velocities as they can be easily detected in experiments by tracking soliton trajectories; see also the previous theoretical developments on the application of IST perturbation theory on the decay of soliton complexes [36, 37]. Fig. 1(a) shows typical behaviour of $\text{RF}(\tau)_n^{\text{im}}$ for a pulse containing three solitons affected by imaginary perturbation. The RFs being formed by discrete spectrum eigenfunctions of the ZS system exhibit an oscillatory behavior inside the box-shaped potential with zero value in its center. The number of the RFs maxima/minima increases with decreasing $\eta_n$, i.e., soliton amplitude. Importantly, for an unperturbed box-shaped pulse, solitons are not individualized due to their strong overlap and interaction. Nevertheless, the IST theory enables manipulation of individual solitons by applying localized perturbations.

The soliton eigenvalue response is proportional to the RFs integral in the perturbation area; see the gray strip in Fig. 1(a). This visual interpretation allows one to choose the areas of the pulse, where perturbations contribute the most, so one can selectively manipulate soliton velocities. As an example Fig. 1(c) shows the pulse disintegration into three separate solitons moving along trajectory well fitted by the theoretical formula (8) together with Eq. (4). The perturbed pulse exhibits dynamics of strongly interacting overlapped solitons, which
<figure><img src="image_5.png" /><figcaption>FIG. 3: Experimental propagation of light pulses over 1200 km: (a, e) no phase perturbation, (b-d, f-h) different locations of the phase perturbation. Panels structure is similar to Fig. 2(c, d). Dashed lines show soliton trajectories obtained from theoretical expressions (8). In panels (b) and (d) we also added solid trajectories corresponding to the eigenvalues obtained from experimental IST of the perturbed pulse.</figcaption></figure>

asymptotically separate from each other. Note that, the perturbation is applied at a chosen position so that the two solitons with the smallest amplitudes receive nearly opposite significant velocities. In contrast, the soliton of the largest amplitude acquires a weakly positive velocity.

To verify the RFs approach we use the experimental setup represented schematically in Fig. 2(a). Here a continuous wave (CW) laser at 1550 nm (APEX-AP3350A) is modulated in both intensity and phase by two electro-optic modulators EOM$_{I-\phi}$ (iXblue MX-LN-20 and iXblue MPZ-LN-10). The setup enables the design of various initial conditions consisting of square-shaped pulses with superimposed localized phase perturbations. The phase perturbation of the field $\psi \rightarrow \psi e^{i\phi} \simeq \psi(1+i\phi)$, allows us to generate almost imaginary perturbations triggering the velocity of solitons [Eq. (8)]. The designed initial conditions are launched into a recirculating fiber loop consisting of an 8 km-long single-mode fiber (SMF) closed back on itself with a 90/10 fiber coupler. Light circulating in the loop is partially extracted at each round trip through the 10% arm of the fiber coupler and detected by a calibrated fast photodiode (PD) (Finisar - XPDV2120R-10-20-040) coupled to an 80 GSa/s sampling oscilloscope (Teledyne Lecroy - LabMaster 10-Zi 65 GHz) leading to a detection bandwidth around 36 GHz. Additionally, the losses accumulated over one circulation in the fiber loop are partially compensated by a counter-propagating Raman pump coupled in and out of the loop via wavelength division multiplexers (WDM), reducing the effective power decay rate of the circulating field to 0.0017 dB/km or $\alpha = 4 \times 10^{-7} \text{ m}^{-1}$, see Fig. 2(b). Our data analysis enables the reconstruction of the spatiotemporal dynamics of the intensity recorded in single-shot [11, 24, 25]. Additionally, the phase dynamics of the initial conditions is measured by using a 3 × 3 coupler-

based interferometric technique; see [38] and Supplemental Material [34]. The full characterization of the initial light field in terms of phase and amplitude allows us to compute the experimental IST spectrum.

In the experiments, the propagation of square-shaped pulses with $P_0 = 20 \text{ mW}$ and full width at half maximum (FWHM) of 220 ps perturbed by a Gaussian-shaped localized phase modulation (FWHM ~ 45ps) was recorded over a distance of 1200 km. Figs. 2(c,d) show 150 round trips inside the fiber loop, that corresponds to the propagation of a pulse with FWHM = 8.75 over 30 nonlinear lengths in normalized units. The generated pulse is not a sharp box due to the limited bandwidth of EOMs; however it is well approximated by a super-Gaussian function, which is close to our theoretical model, see Fig. 2(d,f) and Supplemental Material [34]. The spatiotemporal dynamics recorded in the experiment [Fig. 2(c)] is very similar to the one obtained in numerical simulations of the full NLSE model (1) [Fig. 3(e)] and demonstrates the disintegration of the pulse induced by the phase perturbation. The discrete IST spectra of the perturbed pulse measured in the experiment and its counterpart from simulations are plotted with circles in Figs. 2(g,h).

Theoretical description of the experiments starts from computing perturbations to the soliton eigenvalues of the initial pulse $\psi(\tau)$, using a box potential of the same power value $\int |\psi|^2 d\tau$, see Fig. 2(d). To account for the imperfections of the box-shaped light generation, we use experimentally measured soliton eigenvalues of the unperturbed pulse as the initial IST spectrum, to which we add the theoretical corrections in accordance with the applied phase perturbation; see Supplemental Material [34] for details of the data analysis procedures. Our theory predicts that the largest and smallest solitons ac-
quire similar negative velocities. In contrast, the middle-amplitude soliton acquires a positive velocity, so the perturbed pulse decays into two parts. Experimental results match very well our theoretical predictions; see Fig. 2(c) for the spatiotemporal diagram of the light propagation and the theoretical trajectories. In particular, the left detached pulse exhibits strong internal oscillations produced by the interaction of two solitons moving together. In addition, we investigate experimentally the influence of the perturbation location on the velocities of individual solitons [Fig. 3], which includes a reference example of an unperturbed pulse [Figs. 3(a)]. The latter demonstrates the standard oscillating dynamics similar to the theoretical one shown in Fig. 1(b), moderately affected by the power loss. Fig. 3(b) shows an analog of the experiment from Fig. 2(c) in the case when the perturbation is closer to the pulse center. Meanwhile, Fig. 3(d) is another scenario when the middle-amplitude soliton detaches to the left while the remaining two propagate together to the right. As expected, perturbation located at the pulse center does not disintegrate the light pulse since the functions $RF_n^{im}$ are anti-symmetric and thus, solitons do not acquire velocities, see Fig. 3(c). Note that in Figs. 3(b,d), the smallest amplitude solitons exhibit moderate discrepancy with theory, which can be mainly attributed to second-order corrections of the perturbation theory. We confirm this conjecture by presenting theoretical trajectories computed from experimental IST eigenvalues, see solid lines in Figs. 3(b,d), which considerably better fit the experimentally observed data.

Our approach highlights that the wave functions of the ZS system play not only an auxiliary role in the IST construction but also are directly linked to the physical properties of the nonlinear wave fields governed by the NLSE. Like the probability amplitude in quantum mechanics formalism [39], the constructed RFs tell us how perturbations address certain solitons even though the latter are completely delocalized within the pulse due to strong overlapping and interaction. Further studies of the wave function properties will, we believe, illuminate the description of complex scenarios of the NLSE dynamics [40–44], benefit to the perspective IST-based optical telecommunications approaches [45–47] and contribute to the rapidly developing field of integrable turbulence and dense soliton gases [25, 48–54]. The proposed here RFs method of examining nonlinear pulses sensitivity is very general and can be applied to other integrable systems, such as the Korteweg–de Vries or sine-Gordon equations [15], as well as other wave shapes, e.g., to the sech-pulse, when the wave functions are known [36]. When wave function solutions are unavailable, the RFs can be evaluated numerically and analyzed graphically, extending the capacity of the IST perturbation theory [21, 22, 55, 56] along with the recently advanced semi-analytical IST approaches [57]. Putting forward experimental techniques for control and manipulation of solitons in optical fibers [10–12], we emphasize that they can be also applied in nearly integrable hydrodynamical or Bose-Einstein condensate systems, see e.g. [52, 58, 59].

# SUPPLEMENTARY MATERIALS

## Eigenvalue Response Functions

Here we derive complete theoretical expressions for the eigenvalue response functions (RFs) of the box-shaped potentials. We employ the corresponding wave functions $\Phi(\tau, \lambda) = (\phi_1, \phi_2)^T$ of the Zakharov-Shabat (ZS) system. For convenience, we write the latter in the following form, which is equivalent to the Eq. (3) of the main Letter,

$$ \Phi_{\tau} = \begin{pmatrix} -i\lambda & \psi(\tau) \\ -\psi^*(\tau) & i\lambda \end{pmatrix} \Phi. \qquad (9) $$

For the box-shaped potential, $\psi_{\perp}(\tau)$ of length L and amplitude A, the wave function $\Phi(\tau)$ can be obtained by solving the scattering problem outside and inside the box and matching the solution at $\tau = \pm L/2$, see Ref. [26] of the main Letter, leading to the following result at an arbitrary spectral parameter $\lambda$:

$$ \phi_1 = \begin{cases} e^{-i\lambda\tau}, & \text{if } \tau < -L/2 \\ \{-i\lambda \sin(s)\chi^{-1} + \cos(s)\}e^{iL\lambda/2}, & \text{if } -L/2 < \tau < L/2 \\ a(\lambda)e^{-i\lambda\tau}, & \text{if } \tau > L/2 \end{cases} \qquad (10) $$

$$ \phi_2 = \begin{cases} 0, & \text{if } \tau < -L/2 \\ -A\chi^{-1}\sin(s)e^{iL\lambda/2}, & \text{if } -L/2 < \tau < L/2 \\ b(\lambda)e^{i\lambda\tau}, & \text{if } \tau > L/2 \end{cases} \qquad (11) $$

where

$$ \chi = \sqrt{A^2 + \lambda^2}, \quad s = (L/2 + \tau)\chi, \qquad (12) $$
while $a(\lambda)$ and $b(\lambda)$ are the so-called scattering coefficients, see Ref. [14-15] of the main Letter for details on the IST theory formalism,

$$a(\lambda) = e^{iL\lambda}\{\cos(\chi L) - i\lambda \sin(\chi L)\chi^{-1}\}, \qquad (13)$$

$$b(\lambda) = -A \sin(\chi L)\chi^{-1}. \qquad (14)$$

According to Eq. (8) of the main Letter, the RFs corresponding to real-valued and imaginary-valued perturbations of the box-shaped pulse evaluated at the discrete eigenvalue points $\lambda = \lambda_n$ are given by,

$$\mathrm{RF}_{n}^{\mathrm{re}}(\tau) = -\frac{i}{2} \left. \frac{\phi_{1}^{2} + \phi_{2}^{2}}{\int_{-\infty}^{\infty} \phi_{1}\phi_{2} d\tau} \right|_{\lambda=\lambda_{n}}, \qquad \mathrm{RF}_{n}^{\mathrm{im}}(\tau) = -\frac{1}{2} \left. \frac{\phi_{1}^{2} - \phi_{2}^{2}}{\int_{-\infty}^{\infty} \phi_{1}\phi_{2} d\tau} \right|_{\lambda=\lambda_{n}}. \qquad (15)$$

As one can see, from the structure of the wave function components $\phi_1(\tau)$ and $\phi_2(\tau)$ are provided above in Eqs. (10), (11) the functions $\mathrm{RF}_n^{\mathrm{re}}(\tau)$ and $\mathrm{RF}_n^{\mathrm{im}}(\tau)$ are always pure imaginary-valued and pure real-valued. Note that according to the IST theory, discrete eigenvalues are zeros for the scattering coefficient (13), i.e. $a(\lambda)|_{\lambda=\lambda_n} = 0$, that can be also seen from the eigenvalue relation (6) provided in the main Letter. In addition, the second scattering coefficient $b(\lambda)|_{\lambda=\lambda_n} = (-1)^n$. Evaluating integrals in expressions (15) at $\lambda_n = i\eta_n$, and after some transformations we obtain the following closed-form expressions for the box-shaped potential response functions,

$$\mathrm{RF}_{n}^{\mathrm{re}}(\tau) = -\frac{i}{\Delta_n} \begin{cases} e^{2\eta_n \tau}, & \text{if } \tau < -L/2 \\ \left\{ \left(\cos[s_n(\tau)] + \eta_n \chi_n^{-1} \sin[s_n(\tau)]\right)^2 + A^2 \chi_n^{-2} \sin^2[s_n(\tau)] \right\} e^{-\eta_n L}, & \text{if } -L/2 < \tau < L/2 \\ e^{-2\eta_n \tau}, & \text{if } \tau > L/2 \end{cases} \qquad (16)$$

$$\mathrm{RF}_{n}^{\mathrm{im}}(\tau) = -\frac{1}{\Delta_n} \begin{cases} e^{2\eta_n \tau}, & \text{if } \tau < -L/2 \\ \left\{ \left(\cos[s_n(\tau)] + \eta_n \chi_n^{-1} \sin[s_n(\tau)]\right)^2 - A^2 \chi_n^{-2} \sin^2[s_n(\tau)] \right\} e^{-\eta_n L}, & \text{if } -L/2 < \tau < L/2 \\ -e^{-2\eta_n \tau}, & \text{if } \tau > L/2 \end{cases} \qquad (17)$$

where $\chi_n = \chi(\lambda_n)$, $s_n(\tau) = (L/2 + \tau)\chi_n$, and the normalization constant coefficient $\Delta_n = 2 \int_{-\infty}^{\infty} \phi_1 \phi_2 d\tau|_{\lambda=\lambda_n}$ is,

$$\Delta_n = -A \chi_n^{-2} e^{-\eta_n L} \left[ \sin^2 L \chi_n - \frac{1}{2} \eta_n \chi_n^{-1} \sin 2L \chi_n + \eta_n L \right]. \qquad (18)$$

## Numerical methods

To locate soliton eigenvalues of the box-shaped potentials [Fig. 1(a) and Fig. 2(h) of the main Letter], we solve numerically the transcendental Eq. (6) of the main Letter at fixed parameters A and L. In the case of experimental wave fields [Fig. 2(g) of the main Letter] we use the standard Fourier collocation method, allowing one to solve the ZS eigenvalue problem numerically; see Ref. [22] of the main Letter for details on its implementation. Being fast and fairly accurate, the Fourier collocation approach is based on the Fourier decomposition of the wave field and transformation of the ZS differential equations [Eq. (9)] into a matrix algebraic equations. The latter provides a straightforward access to the ZS system eigenvalues. Alternatively, soliton eigenvalues can be found using more advanced IST algorithms, see Refs. [27-30] in the main Letter.

The conservative NLSE evolution presented in Fig. 1(b,c) of the main Letter was simulated using standard pseudospectral second-order split-step (Fourier) method. The pulse discretization step was $\Delta\tau = 0.04$ providing accurate description of the wavefield Fourier spectra, while the propagation step $\Delta z$ was chosen in order to satisfy the algorithm stability criteria $\Delta z < \Delta\tau^2/4$. For the NLSE evolution with losses, presented in Fig. 2(e) of the main Letter, the same method was used with the addition of an exponential decay term, $\exp(-\alpha\Delta z/2)$, affecting the wavefield as it numerically propagates, with $\alpha = 4 \times 10^{-7} \mathrm{m}^{-1}$ or $\alpha = 0.016$ in normalized units. The discretization and propagation steps for the pulse shown in Fig. 2(f) of the main Letter were chosen similar to the conservative NLSE case.

## Phase measurement using 3 × 3 coupler

The ideal 3 × 3 coupler distributes the injected intensity equally across each of the input ports to each of the output ports. We denote the electric field at the output or input port j as $E_j^{\mathrm{out}}$ and $E_j^{\mathrm{in}}$. Each output consists of its own
<figure><img src="image_6.png" /><figcaption>FIG. 4: (a) Schematic experimental recirculating fiber loop setup. (b) Experimental setup for the phase measurement of the initial condition. EOMI: Intensity electro-optic modulator. EOMφ: Phase electro-optic modulator. WDM: Wavelength-division multiplexing. SMF: Single mode fiber. PD: Photodiode. PC: Polarization controller.</figcaption></figure>

respectively input and the phase-shifted contribution of $2\pi/3$ from the other two inputs. For example, the output from port 1 is given by:

$$E_1^{\text{out}} = \frac{1}{\sqrt{3}} [E_1^{\text{in}} + E_2^{\text{in}} e^{i2\pi/3} + E_3^{\text{in}} e^{i2\pi/3}]. \quad (19)$$

In our setup, port 1 is used for the signal, port 2 for the reference, and port 3 is unused, simplifying the system of $E_j^{\text{out}}$. Considering the relative phase $\Delta\phi$ associated with the signal, we obtain the relations for the intensities:

$$\left\{
\begin{aligned}
I_1^{\text{out}} &= \frac{1}{3} \left[ I_1^{\text{in}} + I_2^{\text{in}} + 2\sqrt{I_1^{\text{in}} I_2^{\text{in}}} \cos(\Delta\phi - 2\pi/3) \right], \\
I_2^{\text{out}} &= \frac{1}{3} \left[ I_1^{\text{in}} + I_2^{\text{in}} + 2\sqrt{I_1^{\text{in}} I_2^{\text{in}}} \cos(\Delta\phi + 2\pi/3) \right], \\
I_3^{\text{out}} &= \frac{1}{3} \left[ I_1^{\text{in}} + I_2^{\text{in}} + 2\sqrt{I_1^{\text{in}} I_2^{\text{in}}} \cos(\Delta\phi) \right].
\end{aligned}
\right.
\quad (20)$$

Combining these expressions we obtain the following relation for the phase in the signal's frame of reference,

$$\tan(\Delta\phi) = \frac{\sqrt{3}}{2} \left[ \frac{I_1^{\text{out}} - I_2^{\text{out}}}{I_3^{\text{out}} - I_1^{\text{out}}/2 - I_2^{\text{out}}/2} \right]. \quad (21)$$

The phase term appears in the tangent function, which is periodic and has discontinuities at $\pm\pi/2$, leading to function jumps when the phase changes monotonically. Using Eq. (21) we find expressions for $\cos(\Delta\phi)$ and $\sin(\Delta\phi)$,

$$\sin(\Delta\phi) = \frac{\sqrt{3}}{2} \left[ \frac{I_1^{\text{out}} - I_2^{\text{out}}}{\sqrt{I_1^{\text{in}} I_2^{\text{in}}}} \right], \quad \cos(\Delta\phi) = \frac{I_3^{\text{out}} - I_1^{\text{out}}/2 - I_2^{\text{out}}/2}{\sqrt{I_1^{\text{in}} I_2^{\text{in}}}}. \quad (22)$$

Now, with the expressions for $\cos(\Delta\phi)$, $\sin(\Delta\phi)$ and $\tan(\Delta\phi)$, we impose conditions on the relative phase measurements and unfold the phase value across the entire range of the pulse amplitudes.

In Fig. 5 we represent the cosine, sine, and tangent functions over the interval $[-\pi; \pi]$ with three colored zones separated by the discontinuities of the tangent function at $\pm\pi/2$. In the central zone (white), the cosine always has a positive value, while in the left (blue) and right (red) zones, it is negative. The sine, on the other hand, is positive from $[0, \pi]$ and negative from $[-\pi, 0]$. This representation allows us to determine the phase value implementing the following procedure,

* Left (blue) zone: adding $\pi$ to the phase measurement.

* Central (white) zone: no phase measurement corrections.

* Right (red) zone: subtracting $\pi$ from the phase measurement.

Using this procedure, we avoid phase jumps, enabling an unambiguous measurement (up to a general constant factor) of the whole light pulse phase modulations.
<figure><img src="image_7.png" /><figcaption>FIG. 5: Behaviour of the cosine, sine, and tangent functions in the three colored zones, separated by the discontinuities of the tangent function associated with the phase unwrapping conditions (21). To enable unambiguous finding of the light pulse phase modulations we add $\pm \pi$ to the phase measurement in the left (blue) and right (red) zones.</figcaption></figure>

## Experimental data processing

Fig. 6 illustrates our data processing protocol applied to each initial experimental pulse. The protocol enables the characterization of the initial condition in terms of its discrete IST spectrum and the eigenvalue response functions (RFs). In the first step, see Fig. 6(a), we fit the signal with a super-Gaussian function $f(T)$, extracting the peak power of our pulses $P_0$, which we use for normalizing the experimental data,

$$f(T) = P_0 \exp \left( - \left( \frac{(T - T_0)^2}{2\sigma^2} \right)^P \right) \qquad (23)$$

In the second step, see Fig. 6(b); since the sensitivity functions are defined for a sharp box-shaped potential, we find an equivalent box-shaped version of our experimental super-Gaussian fit. We use the equality of the potential surfaces as a criterion because it preserves the number of solitons, $N = \text{Integer}[1/2 + AL/\pi]$. For the pulse shown in Fig. 6, we obtain the length of the equivalent box $L_{\text{box}} = 8.75$.

Then, we compute the IST spectrum of the super-Gaussian and box fits numerically and compare their soliton composition as shown in Fig. 6(c). Both potentials have three solitons with no relative velocity. The amplitudes of the super-Gaussian fit potential solitons are slightly lower than their box potential counterpart. This feature is more pronounced for the smallest soliton as it is more sensitive to perturbations. Finally, using the eigenvalues from the box-shaped potential, we retrieve the eigenvalue RFs, which we use for theoretical predictions of soliton trajectories in the main Letter, see Fig. 6(d).

## Soliton trajectories

In our experiments we perturb the initial field $\psi(\tau)$ by phase modulation $\phi(\tau)$. The modulation is characterized by a Gaussian shape with a full width at half maximum (FWHM) of 48 ps $\approx 1.65 \tau$ and an average amplitude of 0.55 rad. We perform the phase measurement independently from the recording of the spatio-temporal dynamics. Then we combine two separate recordings and reconstruct the initial complex field. Finally, we obtain the perturbation term $\delta\psi(\tau)$ as,

$$\delta\psi(\tau) = \psi(\tau)_{\text{pt}} - \psi(\tau) = \psi(\tau)e^{i\phi(\tau)} - \psi(\tau), \qquad (24)$$

and complete the data needed for the prediction of soliton trajectories.

The perturbation $\delta\psi(\tau)$ is primarily imaginary-valued for small-amplitude phase modulations $\phi(\tau)$; however, it also has a real-valued component that we account for when calculating the deviation of the eigenvalues. We obtain two types of the IST spectrum – theoretical and experimental, shown in Fig. 2(g,h) of the main Letter by dots and circles. The dots (theory) correspond to the following two-step procedure: i) first we compute initial IST spectrum
<figure><img src="image_8.png" /></figure>

FIG. 6: Experimental data processing protocol. (a) In red line: initial light pulse (at $z = 0$ km). In black dotted line: the super-Gaussian fit with (23). (b) In black dotted line: the super-Gaussian fit potential obtained from the fit with $P = 2.5$. In black line: the equivalent box-shaped potential. (c) IST spectrum: dots for the super-Gaussian potential (smooth box), crosses for the box-shaped potential. (d) Eigenvalue response functions (RFs) for imaginary perturbations associated to each soliton of the box-shaped potential. The color scheme is the same as used for the IST spectrum in (c).

for absolute-valued intensity profiles leading to eigenvalues alight along imaginary axis, ii) then with the eigenvalue response functions we predict theoretically the soliton eigenvalue deviations corresponding to the perturbation $\delta\psi(\tau)$ and add them to the initial IST spectrum. The circles (experiment) are experimental soliton eigenvalues, obtained directly from the complex-valued field $\psi(\tau)_{pt}$.

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

FIG. 7: Spatiotemporal diagrams of the perturbed pulse propagation from Fig. 2(c,e) of the main Letter. In addition to theoretical trajectories of soliton propagation (dashed lines), solid lines show the trajectories obtained from direct numerical computation of the IST spectrum.

Once the real parts of the perturbed soliton eigenvalues is found, we compute soliton velocities (see Eq. (4) of the
main Letter) and plot soliton trajectories. The incline of the trajectories corresponds to soliton velocities, while their
origin is chosen so to fit position of solitons at large distance where they well separated ($z \gtrsim 15$).

ACKNOWLEDGMENTS

This work of PS, FC, AM, SR has been partially supported by the Agence Nationale de la Recherche through the LABEX CEMPI project (ANR-11-LABX-0007), the SOGOOD (SOGOOD ANR-21-CE30-0061) and StormWave (StormWave ANR-21-CE30-0009) projects, 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 Wavetech). PS, FC, AM and SR thank the Centre d'Etudes et de Recherches Lasers et Applications CERLA for the technical help and and equipment. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for its support and hospitality during the program "Emergent phenomena in nonlinear dispersive waves," during which significant progress was made on this work. The work of RM on the response functions derivation was supported by RSF Grant No. 19-79-30075-II. The work of AG was funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101033047.

* Corresponding author : agelash@gmail.com; Current address: Institute of Physics, Swiss Federal Institute of Technology Lausanne (EPFL), CH-1015 Lausanne, Switzerland.

[1] L. F. Mollenauer, R. H. Stolen, and J. P. Gordon, Physical Review Letters **45**, 1095 (1980).

[2] L. Mollenauer and K. Smith, Optics letters **13**, 675 (1988).

[3] Y. S. Kivshar and G. P. Agrawal, *Optical solitons: from fibers to photonic crystals* (Academic press, 2003).

[4] A. Hasegawa, *Optical solitons in fibers* (Springer Science & Business Media, 2013).

[5] A. Hasegawa and F. Tappert, Applied Physics Letters **23**, 142 (1973).

[6] A. Hasegawa and F. Tappert, Applied Physics Letters **23**, 171 (1973).

[7] M. Remoissenet, *Waves called solitons: concepts and experiments* (Springer Science & Business Media, 2013).

[8] J. M. Dudley, G. Genty, and S. Coen, Reviews of modern physics **78**, 1135 (2006).

[9] D. V. Skryabin and A. V. Gorbach, Reviews of Modern Physics **82**, 1287 (2010).

[10] J. K. Jang, M. Erkintalo, S. Coen, and S. G. Murdoch, Nature communications **6**, 7370 (2015).

[11] P. Suret, M. Dufour, G. Roberti, G. El, F. Copie, and S. Randoux, Physical Review Research **5**, L042002 (2023).

[12] N. Englebert, C. Simon, C. M. Arabí, F. Leo, and S.-P. Gorza, arXiv preprint arXiv:2406.12848 (2024).

[13] V. E. Zakharov and A. B. Shabat, Sov. Phys. JETP **34**, 62 (1972).

[14] M. J. Ablowitz and H. Segur, *Solitons and the inverse scattering transform*, Vol. 4 (Siam, 1981).

[15] S. Novikov, S. Manakov, L. Pitaevskii, and V. Zakharov, *Theory of solitons: the inverse scattering method* (Springer Science & Business Media, 1984).

[16] A. Osborne, *Nonlinear ocean waves* (Academic Press, 2010).

[17] S. Wabnitz, Optics letters **18**, 601 (1993).

[18] V. Afanasjev and N. Akhmediev, Optics letters **20**, 1970 (1995).

[19] T. Okamawari, A. Hasegawa, and Y. Kodama, Physical Review A **51**, 3203 (1995).

[20] V. Gerdjikov, D. Kaup, I. Uzunov, and E. Evstatiev, Physical review letters **77**, 3943 (1996).

[21] Y. S. Kivshar and B. A. Malomed, Reviews of Modern Physics **61**, 763 (1989).

[22] J. Yang, *Nonlinear waves in integrable and nonintegrable systems*, Vol. 16 (SIAM, 2010).

[23] R. Jenkins and K. D. T.-R. McLaughlin, Communications on Pure and Applied Mathematics **67**, 246 (2014).

[24] A. E. Kraych, P. Suret, G. El, and S. Randoux, Physical review letters **122**, 054101 (2019).

[25] A. E. Kraych, D. Agafontsev, S. Randoux, and P. Suret, Physical Review Letters **123**, 093902 (2019).

[26] S. Manakov, Zh. Eksp. Teor. Fiz **65**, 10 (1973).

[27] G. Boffetta and A. R. Osborne, Journal of Computational Physics **102**, 252 (1992).

[28] S. Wahls and H. V. Poor, IEEE Transactions on Information Theory **61**, 6957 (2015).

[29] A. Vasylchenkova, J. E. Prilepsky, and S. K. Turitsyn, Optics letters **43**, 3690 (2018).

[30] R. Mullyadzhanov and A. Gelash, Optics Letters **44**, 5298 (2019).

[31] D. J. Kaup, SIAM Journal on Applied Mathematics **31**, 121 (1976).

[32] V. I. Karpman and E. M. Maslov, Soviet Physics JETP **46**, 281 (1977).

[33] R. Mullyadzhanov and A. Gelash, Physical Review Letters **126**, 234101 (2021).

[34] See *Supplemental Material* at http://link.aps.org/supplemental for details.

[35] M. Klaus and J. K. Shaw, Physical Review E **65**, 036607 (2002).

[36] J. Satsuma and N. Yajima, Progress of Theoretical Physics Supplement **55**, 284 (1974).
[37] J. E. Prilepsky and S. A. Derevyanko, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 75, 036616 (2007).

[38] B. Kelleher, D. Goulding, B. Baselga Pascual, S. P. Hegarty, and G. Huyet, The European Physical Journal D 58, 175 (2010).

[39] L. D. Landau and E. M. Lifshitz, *Quantum Mechanics: Non-relativistic Theory*. V. 3 of Course of Theoretical Physics (Pergamon Press, 1958).

[40] N. Akhmediev, J. M. Soto-Crespo, and A. Ankiewicz, Phys. Lett. A 373, 2137 (2009).

[41] M. Onorato, S. Residori, U. Bortolozzo, A. Montina, and F. T. Arecchi, Phys. Rep. 528, 47 (2013).

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

[43] M. Conforti, S. Li, G. Biondini, and S. Trillo, Optics letters 43, 5291 (2018).

[44] A. Bendahmane, G. Xu, M. Conforti, A. Kudlinski, A. Mussot, and S. Trillo, Nature Communications 13, 3137 (2022).

[45] S. T. Le, V. Aref, and H. Buelow, Nature Photonics 11, 570 (2017).

[46] S. K. Turitsyn, J. E. Prilepsky, S. T. Le, S. Wahls, L. L. Frumin, M. Kamalian, and S. A. Derevyanko, Optica 4, 307 (2017).

[47] L. L. Frumin, A. Gelash, and S. K. Turitsyn, Physical Review Letters 118, 223901 (2017).

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

[49] P. Walczak, S. Randoux, and P. Suret, Physical review letters 114, 143903 (2015).

[50] E. N. Pelinovsky, E. Shurgalina, A. Sergeeva, T. G. Talipova, G. El, and R. Grimshaw, Physics Letters A 377, 272 (2013).

[51] J. M. Soto-Crespo, N. Devine, and N. Akhmediev, Physical review letters 116, 103901 (2016).

[52] I. Redor, E. Barthélemy, H. Michallet, M. Onorato, and N. Mordant, Physical review letters 122, 214502 (2019).

[53] A. Slunyaev and T. Tarasova, Chaos: An Interdisciplinary Journal of Nonlinear Science 32 (2022).

[54] P. Suret, S. Randoux, A. Gelash, D. Agafontsev, B. Doyon, and G. El, Physical Review E 109, 061001 (2024).

[55] D. Kaup, SIAM Journal on Applied Mathematics 31, 121 (1976).

[56] V. Karpman and E. Maslov, Zh. Eksp. Teor. Fiz 73, 537 (1977).

[57] A. Gelash, S. Dremov, R. Mullyadzhanov, and D. Kachulin, Physical Review Letters 132, 133403 (2024).

[58] A. Chabchoub, N. Hoffmann, M. Onorato, G. Genty, J. M. Dudley, and N. Akhmediev, Physical review letters 111, 054104 (2013).

[59] 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).