# Langevin's model for soliton molecules in ultrafast fiber ring laser cavity: Investigating experimentally the interplay between noise and inertia

Anastasiia Sheveleva, Aurélien Coillet, Christophe Finot, Pierre Colman

► To cite this version:

Anastasiia Sheveleva, Aurélien Coillet, Christophe Finot, Pierre Colman. Langevin's model for soliton molecules in ultrafast fiber ring laser cavity: Investigating experimentally the interplay between noise and inertia. Chaos, Solitons & Fractals, 2025, 197, pp.116472. 10.1016/j.chaos.2025.116472. hal-05060492

HAL Id: hal-05060492

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

Submitted on 8 May 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
<figure><img src="image_2.png" /></figure>

# Langevin's model for soliton molecules in ultrafast fiber ring laser cavity: Investigating experimentally the interplay between noise and inertia

Anastasiia Sheveleva ID*, Aurélien Coillet ID*, Christophe Finot ID*, Pierre Colman ID*

Laboratoire Interdisciplinaire Carnot de Bourgogne ICB UMR CNRS 6303, Université Bourgogne Europe, 9 avenue Alain Savary, Dijon, France

## ARTICLE INFO

**Keywords:**

Langevin model

SINDy regression

Soliton molecule

Fiber laser

Numerical twin

## ABSTRACT

The dynamics of soliton molecules in ultrafast fiber ring laser cavity is strongly influenced by noise. We show how a parsimonious Langevin model can be constructed from experimental data, resulting in a mathematical description that encompasses both the deterministic and stochastic properties of the evolution of the soliton molecules. In particular, we were able to probe the response dynamics of the soliton molecule to an external kick in a sub-critical approach, namely without the need to actually disturb the systems under investigation. Moreover, the noise experienced by the dissipative solitonic system, including its distribution and correlation, can now be also analyzed in details. Our strategy can be applied to any systems where the individual motion of its constitutive particles can be traced; the case of optical solitonic-system laser presented here serving as a proof-of-principle demonstration.

## 1. Introduction

The subtle balance between gain and loss mechanism, and dispersion and nonlinearity in an ultrafast fiber ring laser cavity can result in the spontaneous generation of very specific invariant pulses known as dissipative solitons [1]. A few of these co-propagating pulses can interact closely and form a bound state called the soliton molecule (SM). Like their atomic counterpart in chemistry, SMs can exhibit periodic vibrations of their constitutive parameters [2,3] (e.g. oscillations of the inter-soliton temporal spacing, of the optical phase difference, etc.). These SM systems have attracted much attention the last decade because they reveal the subtle dynamics at play in fiber lasers, and therefore call for more fundamental investigations with also possible applications in communication, and in all optical data processing and storage [4,5]. This will be of course only possible if they can be controlled to an acceptable extent [5]. From a fundamental prospect, SMs are a platform of choice to study the nonlinear dynamics of waves and their interaction. The existence of the SMs vibration is often ascribed to the formation of a limit-cycle attractor inside the laser cavity, which results in self-sustained oscillations even in the absence of any (operator-controlled) forced external excitation of the nonlinear oscillator [3]. In a sense, SMs are a very good example of nonlinear oscillators driven by the competition – or the balance – between dissipation and energy gain mechanisms. But fiber ring cavity lasers are extremely complex objects whose behaviors are set by multiple parameters [6–10]. If large efforts have been dedicated both to the

experimental measurements of the laser's photonic state [11,12] and to the control of the ring laser properties [10,13,14], many unknowns still remain. Therefore quantitative controls and prediction strategies rely most of the time on machine-learning [15,16]. This allows certain degree of control at the cost of intelligibility [17]. A direct consequence is that most advanced optimized control and stabilization techniques still cannot be implemented [18–20]. Another strong limitation is related to the intrinsic stochastic nature of the laser cavity. Because the observed dynamics is actually the result of the balance between various excitation and dissipative mechanisms, it is alas not possible to measure these two effects separately. This in turn creates a halo of uncertainty when comparing theoretical modeling against experiment. Such a limitation is actually found in any system whose evolution laws are not known, and that are driven by noise (namely quasi random events).

We demonstrate in this article that it is possible to create, based on experimental data, a Langevin model [21–23] for the vibration of the inter-soliton distance. Using the example of a fiber laser system, we namely derive a parsimonious set of equations [24–26] that reproduces the deterministic (i.e. latent) dynamics of the SM, which we then complement by a noise model that encompasses the stochastic variability observed in the SM evolution. We show in particular that the reconstructed Langevin model allows the subsequent independent study of the dissipation and driving mechanisms, and of the noise properties. Previous studies involving numerical fitting and other deep learning

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: pierre.colman@u-bourgogne.fr (P. Colman).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Experimental setup. PBS: Polarization Beam Splitter. PC: Polarization Controller. PD: Photodiode. HDF: Highly Dispersive Fiber (total dispersion of -50 ps/nm). (b) Evolution of time separation τ between the two dissipative solitons. Its exhibits a characteristic oscillation with a period T₀ = 117.6 RTs. 〈τ〉 = 3.48 ps. (c) Evolution of the optical phase difference ΔΦ between the two dissipative solitons. (d) Phase portrait in the plane {τ, ˙τ} corresponding to (b).</figcaption></figure>

techniques, particularly in the field of fiber lasers, usually considered the presence of noise as purely detrimental; and that it contains no valuable information. If noise is indeed problematic for the data-driven construction of mathematical models because of the numerical artifacts it could create [27]; it however does not influence much the final result if the numerical aspects are appropriately taken care of [28]. We show here that noise can actually improves the quality and the robustness of the model. For this first proof of principle, we apply our strategy to the vibration pattern of a 2-soliton molecule in an ultrafast fiber ring cavity laser. We show how the Langevin model can be constructed; and the information it provides. As complex nonlinear oscillators where noise plays an important role, SMs constitute an interesting platform to test our method. The conclusions are, however, general and the technique presented here can be readily applied to any other dynamical systems.

## 2. Experimental dataset and SINDy techniques

The experimental platform we consider in Fig. 1-(a) is an ultrafast fiber ring cavity laser composed of a 5-m long fiber ring laser cavity comprising 3 meters of Erbium doped fiber. The mode-locking mechanism comes from nonlinear polarization rotation in the single mode fiber (SMF) and a subsequent filtering by a polarizing beam splitter [29]. This laser emits pulses at about 44 MHz repetition rate, and also can support SMs for higher pump settings. The SMs dynamics is recorded thanks to the Dispersive Fourier Transform technique (DFT) [11,12,30]. For this, the laser output is dispersed in a -50 ps/nm Highly Dispersive Fiber (HDF), also known in the telecom industry as Dispersive Compensating Fiber — DCF. Finally, the DFT spectra are recorded by an 80 GS/s, 40 GHz electrical bandwidth, 8-bit depth oscilloscope. The photodiode has a 25 GHz bandwidth, and acts as a low pass filter for higher frequencies. The vibration pattern of the SM is shown in Fig. 1-(b-c): the parameters of interest are the time delay between the two pulses τ and their optical phase difference ΔΦ. Their evolution is observed over about 8500 round-trips, as limited by the oscilloscope memory. It exhibits a characteristic periodic oscillation that presents however some minor fluctuations [2,12,29]. They are more visible on the phase portrait in Fig. 1-(d) which reveals indeed that the SM does not strictly stick to its limit cycle, but experiences cycle to cycle fluctuations. This results in a thick donut-shaped ring, instead of a well-defined loop. We chose here to use the {τ, ˙τ} representation to point out the applicability of our method to other systems, instead of

the traditional representation {τ, φ} used for SMs. Indeed every other system may not allow the direct measure of two conjugated variables. In a more general context, the base of representation that is chosen may have a direct impact on the complexity of the subsequent mathematical model that is derived from the data, as discussed in more details in the next section.

## 3. The SINDy model

Following the Langevin's formalism, the objective here is to model the experimental data using the following equation:

$$ \frac{d}{dt} \begin{bmatrix} \tau \\ \dot{\tau} \end{bmatrix} = \begin{bmatrix} \dot{\tau} \\ F(\tau, \dot{\tau}) \end{bmatrix} + \begin{bmatrix} 0 \\ \sigma_L(\tau, \dot{\tau}, t) \end{bmatrix} \quad (1) $$

Where $\dot{\tau} = d\tau/dt$, and $\sigma_L$ is the generalized Langevin's force which is a null-centered random force. Note that its intensity could however depend on time and location {τ, ˙τ} in the phase portrait. $F(\tau, \dot{\tau})$ is a nonlinear function that describes the deterministic evolution of the oscillator. It can be well approximated using the Sparse Identification of Nonlinear Dynamic (SINDy) algorithm [24]. The latter consists in finding the minimal set of elementary functions $f_i$ whose combination would approximate F to the desired accuracy. It corresponds to finding the most compact vector $\Xi = [\xi_1, ..., \xi_p]$ that satisfies:

$$ \frac{d}{dt} \begin{bmatrix} \dot{\tau}_1 \\ \dot{\tau}_N \end{bmatrix} = \begin{bmatrix} f_1(\tau_1, \dot{\tau}_1) & \dots & f_p(\tau_1, \dot{\tau}_1) \\ \vdots & \ddots & \vdots \\ f_1(\tau_N, \dot{\tau}_N) & \dots & f_p(\tau_N, \dot{\tau}_N) \end{bmatrix} \Xi \quad (2) $$

The index notation $\tau_i$ represents the time records of $\tau(t_j)$. We define the normalized parameter $u = (\tau - \langle\tau\rangle)/\sigma_\tau$, where $\langle\tau\rangle = 3.48$ ps and $\sigma_\tau = 0.060$ ps. We chose the library of linear functions $f_i$ to be composed of 36 polynomial functions of u and ˙u that span all the possible combinations, up to the order 7. This choice is arbitrary; and its quality could be improved by relying on a physics informed strategy where could also be included specific functions, which are inferred by theory to be potential candidate of the SM evolution. This potential impact of arbitrary user-choices on the final result is completely analogue to the initial choice of the architecture for neural-networks. It is up to the user to ensure that the system used for the training indeed encompasses the function to be approximated. In this prospect, our approach to use only polynomial functions could be understood as performing a Taylor expansion of the actual SM's equation of evolution. We paid
<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) $\langle\tau\rangle = 3.48$ ps. Blue circles: Phase portrait of the SM oscillation in the $\{\tau: \dot{\tau}\}$ plane ($d\tau/dt = \dot{\tau}$). Black line: limit-cycle defined by Eq. (3). Red-dashed line: recovery dynamics for a noiseless model when the SM has been initially kicked from the limit-cycle to the yellow dot, which serves as starting point. (b) Details of the recovery of the red-dashed trajectory presented in (a). The geometrical distance from the limit cycle evolves as an exponential decaying law with a characteristic thermalization time $\tau_{th} = 98.2$ RTs. (c) Comparison between the experimental oscillations (solid blue) and the result of Eq. (3) (dashed red). (d) Evolution of the deviation between the SINDy prediction and the actual experimental evolution. The deviation is measured as the euclidean distance in the phase portrait. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

attention that the typical terms that one could expect to find in a simple anharmonic oscillator subject to friction loss (in particular the terms $u^3$ and $\dot{u}$) were indeed included in the library of test-functions. They were however not selected by the regression model. This indicates that the oscillator we are observing is actually of a different nature from the simple pendulum. That said, if a systems' evolution exhibits characteristics symmetry or scaling law, it is still a wise strategy to include in the library some functions that respect that specific symmetry or scaling law. We like to stress that the SINDy technique does not forcefully impose any specific model, hence avoiding any subsequent miss-interpretation of the system's dynamics. Besides, this technique involves the use of non-greedy algorithm with a $L_1$ penalty. This type of penalty is essential to result in parsimonious models because it tends to zero any irrelevant coefficients. Regression is made here using the Least Angle Regression Algorithm (LARS [31]). Other popular regression algorithms are the LASSO [32] and the Elastic-net [33], to cite a few. It is possible to evaluate the quality of the fit by checking which functions $f_j$ are selected and how their respective coefficients vary with the free fitting parameter. In the present case in Fig. 3, the coefficients are varying smoothly with the tuning free parameter $\lambda$, and the functions with the highest polynomial order are also never selected. We see that for 7 non-zero parameters, the idealized dynamics of the SM can be reproduced with about 3.5% accuracy. Considering that the SINDy algorithm is basically a least square error fitting, moderate zero-average random fluctuations have little impact on the final result, as it has been demonstrated [28]. That said, an increased level of noise would greatly impact the precision of the regression coefficients of lesser magnitude. Performing the regression on isolated segments is essential to check that the SM does no undergo any notable change in its dynamics, as it would be characterized by large changes in the coefficients [34]. For constant dynamics however, the best strategy is nevertheless to perform the regression on the whole dataset so as to reduce to the minimum the impact of noise. The SM under study here can be modeled using the following seven-parameter nonlinear equation:

$$ \frac{d^2 u}{dt^2} = \xi_1 u + \xi_2 u^2 + \xi_3 u^5 + \dot{u}(\xi_4 u + \xi_5 u^2) + \dot{u}^2(\xi_6 u^2 + \xi_7 u^3) \quad (3) $$

$\Xi = (-0.758, -0.073, -0.01, 0.412, -0.171, -0.078, 0.011)$. Solution of Eq. (3) is shown graphically in Fig. 2-(c). The equation indeed mimics well the experimental data, minus its local fluctuations. The periodicity is especially well reproduced. Finer analysis of the

divergence between the experimental evolution and the prediction from the SINDy model are presented in Fig. 2-(d). We see that the error remains about constant for about the first 2000 RTs, hence about 17 oscillations cycles; then it increases by steps. This indicates that the divergence is not caused by a drift, hence the incapacity of SINDy model to reproduce the correct periodicity, but can be ascribed to localized phase-shifts caused by the noise that offsets the experimental cycle. It takes about 8000RTs to reach about a complete cycle opposition. It is interesting to note that the parameter $\xi_1$ in Eq. (3) would correspond to a pure-sine linear oscillation of 112.9 RTs periodicity. Therefore we can conclude, based on the actual average $T_0 = 117.6$ RTs periodicity which is observed over 72 full oscillations (about 8500 Round-Trips), that the nonlinearity contributes notably to the oscillator's dynamics. It is also possible to derive from Eq. (3) the shape of the limit cycle attractor (black line in Fig. 2-(a)). When compared to the actual experimental data, it clearly shows the impact of random fluctuations (noise) that forces the SM to explore a larger area around its limit-cycle. Therefore the fit is not only strictly valid for reproducing the ideal oscillator dynamics (i.e. the limit cycle), but also retains some validity to describe the evolution of the SM when it is driven out of its idealized trajectory. The zone of confidence would correspond to the blue area in Fig. 2-(a) that is actually explored by the SM. Outside of this area the dynamics is inferred by analytical continuation, hence it becomes purely speculative for large deviations.

## 4. SINDy: latent dynamics

A key point in the present analysis is that the noise forces the oscillator out of its equilibrium position. Therefore the oscillator can be observed over a larger area than its limit-cycle without having to actually drive forcefully the oscillator out of its nominal trajectory. SMs are indeed fragile soliton complexes that could be easily destroyed. Aside from resulting to a fit with a larger zone of validity, this subcritical approach has also other advantages. Each noise event can indeed be interpreted as a random kick that excites the system, causing it to deviate from its stable recurrent dynamics. In case these random kicks can be tracked individually, it is then possible to observe the systems' response to a perturbation without having to actually perturb it. As a result this approach is non-invasive, hence it does not risk destroying any fragile dynamics equilibrium that may exist, nor
<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) Evolution of the $\xi_i$ coefficients as function of $\lambda$, the free parameter controlling the LARS regression. (b) Number of non-zero coefficients. (c) Corresponding reconstruction error. $\lambda = 3.7$ is chosen for the present work, corresponding to 7 non-zero parameters.</figcaption></figure>

requires a mean to actually excite the system. That said because the SM is under constant perturbation from noise, it is not possible to observe the complete relaxation dynamics. This would require being able to switch off completely the noise for an extended time, which is not possible experimentally. But using the numerical model we just derived as a numerical twin, it is now possible to investigate how the SM would have relaxed in absence of noise from a position away from the limit-cycle (yellow circle in Fig. 2-(a)): following the noiseless model defined in Eq. (3), the SM will then dissipate the excess of energy and spiral down back to its limit-cycle.

The (normalized) geometrical distance from the limit-cycle is shown in Fig. 2-(b). The SM recovers to its nominal steady dynamics with a relaxation time $\tau_{th} = 98.23$ RTs. This set the typical time-range of excitations which the SM would be able to respond to. Note that this characteristic time is about equals to the actual oscillator periodicity (117 RTs). Faster excitations would leave no time for the oscillator to adapt to them; they could destroy the SM or cause its transfer to another limit-cycle (provided the latter exists) [35]. Similarly, Eq. (3) can be further analyzed to determine the loci of the driving and of the dissipation mechanisms, the strength of the potential confining the SM, etc. Interestingly, for the present case, the sole equilibrium position would correspond to the static molecule ($u = 0 \Leftrightarrow \tau = \langle \tau \rangle = C^{onst}$).

That said, the SINDy analysis still indicates that the SM does not spontaneously oscillate around its static equilibrium position: the sustainability of the current SM system is indeed due to the velocity contribution $u$, which is always repulsive for the small excursion $|\tau - \langle \tau \rangle| < 0.14$ ps experienced here by the SM. As a consequence, we can infer that the existence of the oscillation is actually intrinsic to the existence of an initial drift velocity between the two pulses during the formation of the molecule. An important level of noise can also prevent the molecule to remain static. This rises the question regarding the existence of specific classes of SMs depending on whether the oscillation occurs around a static equilibrium position, or if it is intrinsically linked to the stability of the SM. Answering such a question would require to perform the current analysis on numerous SM records at different but controlled laser's parameters. Note that the physical parameters (e.g. type and parameters of the gain medium, etc.) influence the laser dynamics is a non-trivial manner; and their impact is not known yet quantitatively. This would then require more complex fiber laser involving fine stabilization; and exceed the scope of the present proof-of-principle demonstration.

## 5. Impact of noise on the dynamics: the Langevin model

The SINDy analysis [24] is an important tool in order to infer latent (i.e. deterministic) dynamical features that cannot be measured directly because of the presence of noise. It can also be used to extract stochastic properties.

According to Eq. (1), the Langevin force $\sigma_L$ causes the SM to experience a dynamics different from the latent dynamics, approximated by the SINDy regression, as discussed in the previous session. The real latent dynamics differs of course slightly from the SINDy ones, so that Eq. (1) can be rewritten as:

$$ \frac{d}{dt} \begin{bmatrix} \tau \\ \dot{\tau} \end{bmatrix} = \begin{bmatrix} \dot{\tau} \\ F_{SINDy}(\tau, \dot{\tau}) \end{bmatrix} + \begin{bmatrix} \theta_{\tau}(\tau, \dot{\tau}) \\ \theta_{\dot{\tau}}(\tau, \dot{\tau}) + \sigma_L(\tau, \dot{\tau}, t) \end{bmatrix} \quad (4) $$

Where $\theta(\tau, \dot{\tau})$ is the error between the SINDy regression and the real latent dynamics. Comparing for each round-trip the divergence of evolution ($\delta\tau, \delta\dot{\tau}$) between the actual SM's evolution and the SINDy projection results for small divergences into the following linearized equation:

$$ \frac{d}{dt} \begin{bmatrix} \delta\tau \\ \dot{\delta\tau} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ \partial_{\tau} F_{SINDy} & \partial_{\dot{\tau}} F_{SINDy} \end{bmatrix} \begin{bmatrix} \delta\tau \\ \dot{\delta\tau} \end{bmatrix} + \begin{bmatrix} \theta_{\tau}(\tau, \dot{\tau}) \\ \theta_{\dot{\tau}}(\tau, \dot{\tau}) + \sigma_L(\tau, \dot{\tau}, t) \end{bmatrix} \quad (5) $$

With the observation of the divergence between the experimental data and the SINDy model, Eq. (5) can be solved for $\theta$ and $\sigma_L$. The results are presented graphically in Fig. 4-(a$_1$, b$_1$), where the histograms of $\theta + \sigma_L$ are presented depending on the angular position along the limit cycle (trigonometric orientation, zero angle at East direction). By definition, the Langevin's force is random, and therefore the average of the histogram, hence the average drift, can be attributed to $\theta$ (pink line in Fig. 4-(a$_1$, b$_1$)). The resulting $\sigma_L$ is presented in Fig. 4-(b$_2$): its amplitude and distribution does not depend on the loci. Moreover it exhibits a Gaussian distribution (Fig. 4-(c)), uniform over the whole phase portrait, hence $\sigma_L(\tau, \dot{\tau}, t) = \sigma_L(t)$. In the canonical Langevin model, the remaining distribution Fig. 4-(a$_2$) would be zero. It serves here as an estimate of the numerical artifacts caused by the noise and the subsequent numerical processing. In particular the computation of second order time derivative on noisy signals is known to be a source of problems. These numerical artifacts are however by three times weaker than the Langevin perturbation, so we can conclude that both our numerical treatment and the presented analysis is legitimate. The magnitude of the Langevin's force can be normalized to the dissipation capability of the limit-cycle. Its lesser magnitude, combined with the absence of extreme rogue events, indicates that the current SM is operated in a stable regime, away from any bifurcation point.

## 6. Discussion

We would like to stress that the method we propose here never assumes any peculiar shapes regarding the noise distribution. This is a very practical advantage over other stochastic optimization methods like the Fokker-Planck optimization [36,37]. That said a first limitation of the naive approach used here is the requirement to track the evolution of individual particles. This could be problematic for other systems composed of numerous indiscernible particles. In order to showcase in more details this advantage provided by our method, we generated synthetic data sets by solving the Eq. (3) in presence of controlled noise sources. These datasets were then processed using the same technique as the experimental data. The extracted SINDy coefficients varied slightly depending on the data set, but remained within the natural variability observed for the coefficients of the experimental data. Furthermore no notable changes where observed in the dynamics, namely all SINDy equations demonstrate about the same limit cycle. The noise distributions extracted from the synthetic data are presented in Fig. 5.

The nominal noise distributions used for the data generation are well retrieved, but with more pronounced deformations for the uniform
<figure><img src="image_6.png" /><figcaption>Fig. 4. ($a_1, b_1$) Results of Eq. (5). Cyan line indicates the average drift force, attributed to $\theta_z$ (resp. $\theta_x$). ($a_2, b_2$) Resulting distribution of the Langevin's Force $\sigma_L$, according to Eq. (1). It is obtained after subtraction of the drift $\theta$ in ($a_1, b_1$). By definition $F(\tau, \hat{\tau}) = F_{SINDy}(\tau, \hat{\tau}) + \theta(\tau, \hat{\tau})$. According to the Langevin's formalism, ($a_2$) must be null, and therefore is ascribed to numerical artifact. (c) Blue: Histogram of the Langevin's force shown in ($b_2$). Gray: histogram of the distribution of the numerical errors, constructed from ($a_2$). (d) Correlation $\langle \sigma_L(t+t_0)\sigma_L(t) \rangle_t / \langle \sigma_L(t)^2 \rangle_t$ as a function of the time offset $t_0$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Noise distribution retrieved from synthetic data sets generated using different noise source. Shape of the nominal noise distribution is shown in dashed-magenta. (a) Uniform distribution. (b) Exponential probability distribution. (c) Gaussian distribution. (d) Distribution with a double hump and zero occupation probability at zero Langevin's force. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

(a) or the doubled peaked (d) distribution. We see indeed that the sharp edges tend to be smoothed out, and that the distribution is repopulated around the null noise, respectively. These alteration of the shape are directly related to the numerical treatment undergone by the data. The nominal SINDy process involves numerical derivations, hence the combination of different data-points and their respective noise. As per the central limit theorem, this blending of different noise events would gradually result in an effective normal distribution. Further smoothing of data, as sometimes required for the accurate derivation of noisy data, would eventually worsen further the situation by averaging the noise events over a larger time windows. This potential problem can however be mitigated by more advanced SINDy implementations [38]. Moreover, noise usually presents a characteristic cut-off frequency, directly related to the nature of the process which has generated it. Therefore as long as the time windows required for the numerical processing remains shorter than the noise cut-off frequency, the numerical processing has lesser impact. Synthetic data were generated here in an unfavorable scenario considering true white noise. In the present article, data are recorded at about 44MS/s (i.e. 1S/RT), a pace much faster than the noise dynamics as expected and observed for fiber lasers.

Additional information can be obtained from the noise time correlation, as shown in Fig. 4-(d). In details, we first see at small evolution time a strong positive correlation followed by a strong negative anti-correlation. This provides an estimate of the noise cut-off frequency, which must be around 12 RTs. On longer timescale, the noise retains about 3% correlation/anti-correlation oscillations with about 1.1 MHz periodicity. We ascribe this behavior to the intrinsic dynamic of the gain medium where the evolution of the gain level is inversely correlated to the variation of the optical power circulating within the cavity. Any power drop in the cavity would result into less power taken at each round-trip from the gain medium; in turn the latter can build up extra energy and gain, which results later into an increase of the optical intra-cavity power. Same reasoning works in case of a transient increase of optical power. Negative delayed feedback is known to create oscillations, with typical µs-timescale for ion-doped fibers [39–41]. For cavity laser with active control, this residual oscillation could eventually be decreased further if the right control strategy is employed [42]. This could involve for example modulating the cavity loss with the right frequency and phase offset to equalize the intra-cavity power. And the eventual instability created by the feedback control itself could also
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** (a) Nominal Langevin Model as defined by Eq. (4). (b) Parametric study, at nominal noise level, of the impact of the noise correlation's period and strength. Metric is the probability of the SM to dissociate for a 1000 RTs propagation. (c) red-dash: Result of Langevin simulation where the noise amplitude (i.e. std($\sigma_L$)) has been multiplied by two. The SM breaks apart quickly. (d) red-dash: Result of Langevin simulation where the frequency of the noise correlation has been set to the oscillation period. The SM breaks apart eventually, after a rather long evolution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

be identified. In a nutshell, the information provided by the Langevin regression method proposed here serves as an advanced diagnosis tool, so as to select the best working point, or to identify potential issues. Actual control is possible provided that the SINDy regression is done in function of a few controllable laser parameters. The SM dynamics and the noise features would then be expressed as functions of controllable parameters; and therefore an optimization and control is possible. The possibility to characterize experimentally the strength and distribution of fluctuations would allow to get the SM progressively arbitrarily close (or far away, depending on the application) from any transition. We present graphically in Fig. 6-(b) the result of the full Langevin model Eq. (1) that combines the latent dynamics and the noise model as defined in Fig. 4-(c,d). This can be directly compared with the full experimental data shown in the left panel Fig. 6-(a). We see that the main features are reproduced correctly. The differences come from the over-simplification of the oscillator model through the SINDy regression, where parts of the confinement potential cannot be approximated analytically with enough accuracy. The asymmetry of the confining potential is directly reflected on the skewness of the SM probability of occupation around its limit cycle. This problem can be corrected through the further addition of extra small non-analytical corrective terms into $\theta(\tau, t)$ in Eq. (3). This is subject to the fact that there are enough sampled points at each location to determine the latter with precision. This increase in numerical precision would however come at the cost of the pure analytical description as provided by the SINDy algorithm. This is actually an inherent weakness of the SINDy method where good knowledge of the system is required in order to select the right set of functions for the optimization step. That said the numerical twin which we just created brings also a strong advantage in the sense that it is now possible to experience virtually different configurations and interplay scenario between the noise and the latent dynamics, as presented in Fig. 6-(b-d). Firstly, increasing the amplitude of the Langevin force, hence noise, by a factor two would result in the destabilization of the SM which would break apart 80% of the time within 8500 RTs (Fig. 6-(c)). No breaking is observed at nominal noise level. This indicates that while the SM is stable, a slight increase of the noise would have dramatic impact. Conversely, a low frequency noise, hence constant over a few RTs, would enlarge notably the parameters

space explored by the SM (Not shown here). We then investigated in Fig. 6-(d) the impact of the frequency of the residual noise temporal correlation. At the nominal noise and correlation levels, but with a correlation frequency synchronized with the SM oscillation period, the SM is gradually exploring a larger and larger area. Eventually the SM breaks also apart. Quantitative data can be extracted from these graphics and we present in Fig. 6-(b) the probability rate of the SM to dissociate depending on the Noise correlation's period and strength. We observe a clear cut-off for correlation periods faster than half of an oscillation cycle where the noise correlation strength has little impact. For slower correlation, any increase of the correlation strength degrades the stability of the SM. This impact is especially striking for correlation periods about matching the oscillation ones (synchronous noise driving). We finally note the presence of dark stripes at specific correlations periods indicating that the noise correlation has little impact on the SM for these discrete set of frequencies. If the noise's features are by definition hard to control, such investigation provides still valuables information. It indicates also clearly that by injecting some noise at specific frequency, it is also possible to prevent the binding of the SM; and hence extends the working range of the fiber laser in harmonic-locking. We can conclude from this that any remaining noise correlation, albeit weak, seems to impact greatly the SM dynamics.

## 7. Conclusion

To conclude, we constructed a Langevin model for the accurate description of the vibration pattern of a soliton molecule (SM). Interestingly, the presence of fluctuations does not hinder the regression process, but on the contrary makes it stronger by exploring part of the parameters space outside the limit-cycle of the nonlinear oscillator. The key advantage of our approach is that, in a subsequent analysis after the mathematical model is constructed, it is possible to differentiate between the latent evolution of the nonlinear oscillator and its stochastic properties. Therefore, characteristic features such as the strength of noise or the intrinsic oscillator's relaxation dynamics can be inferred, albeit they cannot be measured independently. Moreover the SM properties are inferred in a sub-critical approach that does not require to actively perturb or probe the laser cavity. But thanks to
the numerical twin that is constructed, it is still possible to explore further the impact of noise, including the onset at which it could alter notably the SM dynamics. These general features are intrinsically linked to the global fitting procedure that we propose, which encompasses both features of the latent dynamics and of the noise. Outside the field of fiber lasers, we believe that the construction of a Langevin model, based on a SINDy regression, could be applied to any system composed of trackable particles. More specific to systems for which the final objective is to control them, performing a parametric Langevin analysis against a few of the parameters that can be controlled would allow a deterministic control. Relying on subsequent optimal control techniques would be a competitive approach to the use of deep-learning techniques. The use of the latter is trendy for the control of black-box systems, but their training requires large data set which is not always possible to obtain. Moreover the reliability of the resulting control cannot be assessed; and any deviation from the system's points of operation tested during the training could lead to failure. Here the langevin's analysis provides estimate for the stability zones; and ensures the absence of uncontrolled transition states while changing dynamically the point of operation. We believe that it would allow better control and comprehension of complex laser systems.

For this first demonstration, we applied our analysis on the {$\tau, \dot{\tau}$} parameter space to illustrate the generality of our method to any situation where a time signal is recorded. Note that in specific situations where conjugated parameters could be directly recorded, it would then remove the necessity to compute a numerical time derivative, which is always a problematic matter in presence of noise, hence improve the overall quality of the analysis.

### CRediT authorship contribution statement

**Anastasiia Sheveleva:** Writing – review & editing, Visualization, Validation, Software, Investigation, Formal analysis, Data curation.
**Aurélien Coillet:** Writing – review & editing, Resources, Data curation.
**Christophe Finot:** Writing – review & editing, Funding acquisition, Conceptualization.
**Pierre Colman:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.

### Funding

This work was supported by the French ANR funded CoMuSim and NAC-NIP (ANR-20-CE24-0007) and OPTIMAL projects, France (ANR-20-CE30-0004).

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

P.C. thanks Omri Gat for fruitful discussion.

### Data availability

Data will be made available on request.

### References

[1] Grelu P, Akhmediev N. Dissipative solitons for mode-locked lasers. Nat Photonics 2012;6:84. http://dx.doi.org/10.1038/nphoton.2011.345.

[2] Igbonacho J, Nithyanandan K, Krupa K, Dinda PT, Grelu P, Moubissi AB. Dynamics of distorted and undistorted soliton molecules in a mode-locked fiber laser. Phys Rev A 2019;99:063824. http://dx.doi.org/10.1103/PhysRevA.99.063824.

[3] Zavyalov A, Iliew R, Egorov O, Lederer F. Discrete family of dissipative soliton pairs in mode-locked fiber lasers. Phys Rev A 2009-05;79:053841. http://dx.doi.org/10.1103/PhysRevA.79.053841.

[4] Pang M, He W, Jiang X, Russell PSJ. All-optical bit storage in a fibre laser by optomechanically bound states of solitons. Nat Photonics 2016-05;10(7):454-8. http://dx.doi.org/10.1038/nphoton.2016.102.

[5] Liu Y, Huang S, Li Z, Liu H, Sun Y, Xia R, et al. Phase-tailored assembly and encoding of dissipative soliton molecules. Light: Sci Appl 2023-05;12(1). http://dx.doi.org/10.1038/s41377-023-01170-x.

[6] Zavyalov A, Iliew R, Egorov OA, Lederer F. Dissipative soliton molecules with independently evolving or flipping phases in mode-locked fiber lasers. Phys Rev A: At Mol Opt Phys 2009;80. http://dx.doi.org/10.1103/PhysRevA.80.043829.

[7] Kutz JN, Collings BC, Bergman K, Knox WH. Stabilized pulse spacing in soliton lasers due to gain depletion and recovery. IEEE J Quantum Electron 1998;34(9):1749-57. http://dx.doi.org/10.1109/3.709592.

[8] Xian T, Wang W, Zhan L. Real time revealing relaxation dynamics of ultrafast mode-locked lasers. Phys Rev Res 2022-03;4(1):013202. http://dx.doi.org/10.1103/physrevresearch.4.013202.

[9] Soto-Crespo JM, Akhmediev N, Grelu P, Belhache F. Quantized separations of phase-locked soliton pairs in fiber lasers. Opt Lett 2003;28(19):1757-9. http://dx.doi.org/10.1364/OL.28.001757.

[10] Jang JK, Erkintalo M, Murdoch SG, Coen S. Ultraweak long-range interactions of solitons observed over astronomical distances. Nat Photonics 2013-07;7(8):657-63. http://dx.doi.org/10.1038/nphoton.2013.157.

[11] Goda K, Jalali B. Dispersive Fourier transformation for fast continuous single shot measurements. Nat Photonics 2013;7:102. http://dx.doi.org/10.1038/nphoton.2012.359.

[12] Herink G, Kurtz F, Jalali B, Solli DR, Ropers C. Real-time spectral interferometry probes the internal dynamics of femtosecond soliton molecules. Science 2017;356(6333):50. http://dx.doi.org/10.1126/science.aal5326.

[13] Nimmesgern L, Beckh C, Kempf H, Leitenstorfer A, Herink G. Soliton molecules in femtosecond fiber lasers: universal binding mechanism and direct electronic control. Optica 2021-10;8(10):1334. http://dx.doi.org/10.1364/optica.439905.

[14] He W, Pang M, Yeh DH, Huang J, Menyuk CR, Russell PSJ. Formation of optical supramolecular structures in a fibre laser by tailoring long-range soliton interactions. Nat Commun 2019-12;10(1). http://dx.doi.org/10.1038/s41467-019-13746-6.

[15] Girardot J, Billard F, Coillet A, Hertz E, Grelu P. Autosetting mode-locked laser using an evolutionary algorithm and time-stretch spectral characterization. IEEE J Sel Top Quantum Electron 2020-09;26(5):1-8. http://dx.doi.org/10.1109/jstqe.2020.2985297.

[16] Woodward RI, Kelleher EJR. Towards 'smart lasers': self-optimisation of an ultrafast pulse source using a genetic algorithm. Sci Rep 2016-11;6(1). http://dx.doi.org/10.1038/srep37616.

[17] Haluszczynski A, Räth C. Controlling nonlinear dynamical systems into arbitrary states using machine learning. Sci Rep 2021-06;11(1). http://dx.doi.org/10.1038/s41598-021-92244-6.

[18] Berkovitz LD, Medhin NG. Nonlinear optimal control theory. Chapman and Hall/CRC; 2012-08. http://dx.doi.org/10.1201/b12739.

[19] Boscain UV, Piccoli B. An Introduction to Optimal Control. 1967. URL https://apisemanticscholar.org/CorpusID:13140527.

[20] Sugny D, Kontz C, Jauslin HR. Time-optimal control of a two-level dissipative quantum system. Phys Rev A 2007-08;76(2):023419. http://dx.doi.org/10.1103/physreva.76.023419.

[21] Zwanzig R. Nonlinear generalized langevin equations. J Stat Phys 1973-11;9(3):215-20. http://dx.doi.org/10.1007/bf01008729.

[22] Darve E, Solomon J, Kia A. Computing generalized langevin equations and generalized Fokker-Planck equations. Proc Natl Acad Sci 2009-07;106(27):10884-9. http://dx.doi.org/10.1073/pnas.0902633106.

[23] Brückner DB, Ronceray P, Broedersz CP. Inferring the dynamics of underdamped stochastic systems. Phys Rev Lett 2020-07;125(5):058103. http://dx.doi.org/10.1103/physrevlett.125.058103.

[24] Brunton SL, Proctor JL, Kutz JN. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc Natl Acad Sci 2016-03;113(15):3932-7. http://dx.doi.org/10.1073/pnas.1517384113.

[25] Floryan D, Graham MD. Data-driven discovery of intrinsic dynamics. Nat Mach Intell 2022-12;4(12):1113-20. http://dx.doi.org/10.1038/s42256-022-00575-4.
[26] Chen B, Huang K, Raghupathi S, Chandratreya I, Du Q, Lipson H. Automated discovery of fundamental variables hidden in experimental data. Nat Comput Sci 2022-07;2(7):433-42. http://dx.doi.org/10.1038/s43588-022-00281-6.

[27] Prokop B, Gelens L. From biological data to oscillator models using SINDy. IScience 2024-04;27(4):109316. http://dx.doi.org/10.1016/j.isci.2024.109316.

[28] Ermolaev AV, Sheveleva A, Genty G, Finot C, Dudley JM. Data-driven model discovery of ideal four-wave mixing in nonlinear fibre optics. Sci Rep 2022-12;12(1). http://dx.doi.org/10.1038/s41598-022-25523-5.

[29] Hamdi S, Coillet A, Grelu P. Real-time characterization of optical soliton molecule dynamics in an ultrafast thulium fiber laser. Opt Lett 2018;43(20):4965-8. http://dx.doi.org/10.1364/OL.43.004965.

[30] Godin T, Sader L, Kashi AK, Hanzard P-H, Hideur A, Moss DJ, et al. Recent advances on time-stretch dispersive Fourier transform and its applications. Adv Phys: X 2022-05;7(1):2067487. http://dx.doi.org/10.1080/23746149.2022.2067487.

[31] Efron B, Hastie T, Johnstone I, Tibshirani R. Least angle regression. Ann Statist 2004-04;32(2). http://dx.doi.org/10.1214/009053604000000067.

[32] Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B Stat Methodol 1996-01;58(1):267-88. http://dx.doi.org/10.1111/j.2517-6161.1996.tb02080.x.

[33] Zou H. The adaptive lasso and its oracle properties. J Amer Statist Assoc 2006-12;101(476):1418-29. http://dx.doi.org/10.1198/0162145060000000735.

[34] Lejarza F, Baldea M. Data-driven discovery of the governing equations of dynamical systems via moving horizon optimization. Sci Rep 2022-07;12(1). http://dx.doi.org/10.1038/s41598-022-13644-w.

[35] Kurtz F, Ropers C, Herink G. Resonant excitation and all-optical switching of femtosecond soliton molecules. Nat Photonics 2019-10;14(1):9-13. http://dx.doi.org/10.1038/s41566-019-0530-3.

[36] Risken H. The Fokker-Planck equation. Springer Berlin Heidelberg; 1996, http://dx.doi.org/10.1007/978-3-642-61544-3.

[37] Cheng XZ, Jalil MBA, Lee HK, Okabe Y. Mapping the Monte Carlo scheme to langevin dynamics: A Fokker-Planck approach. Phys Rev Lett 2006-02;96(6):067208. http://dx.doi.org/10.1103/physrevlett.96.067208.

[38] Kaheman K, Kutz JN, Brunton SL. SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proc R Soc A: Math Phys Eng Sci 2020-10;476(2242). http://dx.doi.org/10.1098/rspa.2020.0279.

[39] Desurvire E. Analysis of transient gain saturation and recovery in erbium-doped fiber amplifiers. IEEE Photonics Technol Lett 1989-08;1(8):196-9. http://dx.doi.org/10.1109/68.36040.

[40] Kuroda K, Shibata K, Yoshikuni Y. Wavelength-dependent transition time of gain saturation in an erbium-doped fiber amplifier. Appl Phys B 2015-04;120(1):111-5. http://dx.doi.org/10.1007/s00340-015-6109-x.

[41] Khamis M. Dynamic investigation of gain-clamped thulium and erbium-doped fiber amplifiers for reconfigurable WDM systems: comparative study. J Opt 2019-11;49(1):76-82. http://dx.doi.org/10.1007/s12596-019-00581-7.

[42] Fasel U, Kaiser E, Kutz JN, Brunton BW, Brunton SL. SINDy with control: A tutorial. In: 2021 60th IEEE conference on decision and control. CDC, IEEE; 2021-12, p. 16-21. http://dx.doi.org/10.1109/cdc45484.2021.9683120.