# Chaotic Internal Dynamics of Dissipative Optical Soliton Molecules

Youjian Song, Defeng Zou, Omri Gat, Minglie Hu, Philippe Grelu

► To cite this version:

Youjian Song, Defeng Zou, Omri Gat, Minglie Hu, Philippe Grelu. Chaotic Internal Dynamics of Dissipative Optical Soliton Molecules. Laser and Photonics Reviews, 2023, 17 (8), 10.1002/lpor.202300066. hal-04875152

HAL Id: hal-04875152

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

Submitted on 20 Jan 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
Chaotic internal dynamics of dissipative optical soliton molecules

Youjian Song, Defeng Zou, Omri Gat,\* Minglie Hu,\* and Philippe Grelu\*

Y. Song, D. Zou

Ultrafast Laser Laboratory, Key Laboratory of Opto-electronic Information Science and Technology of Ministry of Education, School of Precision Instruments and Opto-electronics Engineering, Tianjin University, Tianjin 300072, China

O. Gat

Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel

E-mail: omrigat@mail.huji.ac.il

M. Hu

Ultrafast Laser Laboratory, Key Laboratory of Opto-electronic Information Science and Technology of Ministry of Education, School of Precision Instruments and Opto-electronics Engineering, Tianjin University, Tianjin 300072, China

E-mail: huminglie@tju.edu.cn

Ph. Grelu

Laboratoire ICB UMR 6303 CNRS, Université de Bourgogne, 9 avenue A. Savary, Dijon
21000, France

E-mail: philippe.grelu@u-bourgogne.fr

Keywords: ultrafast lasers, dissipative solitons, chaos, bifurcations, soliton molecules, ultrafast characterization, fiber lasers

**Abstract.** When a laser cavity supports the propagation of several ultrashort pulses, these pulses can form compact bound states called soliton molecules. Soliton molecules are fascinating objects of nonlinear science, presenting striking analogies with their matter molecules counterparts. The relative timing and phase between the copropagating pulses are the most salient internal degrees of freedom of the soliton molecule. The soliton pair, composed of two identical pulses, constitutes the chief soliton molecule of fundamental interest. Its two major internal degrees of freedom allow self-oscillating soliton molecules, which have been
recurrently observed. However, despite theoretical predictions, the low-dimensional chaotic dynamics of a soliton-pair molecule remains elusive. This article reports the observation of chaotic soliton-pair molecules within an ultrafast fiber laser, by means of a direct measurement of the relative optical pulse separation with sub-femtosecond precision in real time. Moreover, it demonstrates an all-optical control of the chaotic dynamics followed by the soliton molecule, by injection of a modulated optical signal that resynchronizes the internal periodic vibration of soliton molecule. The fast error-free switching between ordered and chaotic soliton molecules enabled by pump current sweeping and external injection highlights the potential prospects of all-optical logic gates and chaotic communication using soliton molecules.

# 1. Introduction

Chaotic dynamics are characterized by aperiodicity, high sensitivity to initial conditions and long-term unpredictability. These are intrinsic features and behaviors observable in most nonlinear systems. The manifestation of chaos in extended natural systems raises awe and fear, from the uncertain evolution of celestial bodies to the weather instability, respectively highlighted in the seminal works of Poincaré and Lorenz.[1,2] Therefore, from the past century, there has been a high impetus to investigate chaos within laboratory-confined systems over a wide range of scientific fields, from hydrodynamics and electromagnetics to quantum physics and optics.[3] Studying chaos in nonlinear optics is of particular interest since photonics exhibits handy systems such as the laser oscillator that can advance the fundamental understanding of chaotic dynamics and utilize chaos in technological innovations such as in secure optical communications and Lidar sensing in automotive environments.[4,5] Alternatively, understanding the early warnings signaling the onset of chaos helps circumvent its adverse manifestations in laser applications requiring high stability.[6-8]

The ultrafast laser constitutes an outstanding system for the investigation of chaos, being inherently wideband, highly nonlinear, and offering multiple parameter control schemes. Furthermore, the commonly accepted “dissipative soliton” paradigm explicates the great diversity of ultrashort laser pulse dynamics, which spans from stationary to chaotic pulse emission. Dissipative solitons are attracting states that result from the composite balance between dissipative and dispersive propagation effects and are potentially subjected to a great variety of bifurcations.[9] The standard mode-locked laser regime features a single laser pulse that retakes the same evolution along successive cavity roundtrips: this results from the existence of a stable focus attractor providing long term stability and robustness. Within the parameter space of the optical cavity, the laser can undergo abrupt bifurcations leading to
chaotic pulse dynamics.[10,11] Vivid illustrations of these chaotic dynamics include “soliton explosions”,[12] which manifest as intermittent chaos in the pulse evolution over numerous cavity roundtrips, and “noise-like” pulse dynamics, characterized by a persisting noisy optical waveform and a complete loss of mutual coherence between successive roundtrips.[13-14] Soliton explosions and noise-like pulsing belong to the category of “incoherent dissipative solitons”, which involves chaotic attractors to solve the apparent paradox of an average temporal localization for the pulse and its high instability.[15,16] We emphasize that these chaotic solitons evolve around attractors of extended dynamical systems featuring a high dimensionality, with similar possibilities within Kerr microresonators,[17] noting that the existence of low-dimensional chaotic soliton dynamics has been recently reported.[18,19]

Other abrupt bifurcations may result from the laser multi-pulsing instability (MPI).[20]
Typically, when the pump power is augmented, the MPI leads to an increased number of
intracavity pulses. Therefore, crossing the MPI with subsequent complex bifurcations leads to
a virtually unlimited landscape of complex and chaotic ultrafast multi-pulse dynamics.[21-23] The
interaction between dissipative solitons is also known to lead to the spontaneous formation of
compact bound states, also known as dissipative soliton molecules, since they display striking
analogies with their matter molecule counterparts. Akin to the hydrogen molecule, the soliton-
pair molecule is made of two identical constituents: it is the chief soliton molecule of
fundamental interest, and the focus of our present investigations. Once formed, a stable optical
soliton molecule can propagate indefinitely within the laser cavity. For a soliton molecule being
stationary in its moving frame, the relative temporal separation and phase between the soliton
constituents remain constant. Stationary soliton molecules are routinely generated in laser
oscillators and can be represented as point attractors in an infinite-dimensional phase space.[9]
However, upon a change in the laser parameters, a Hopf-type bifurcation can take place and
generate pulsating soliton molecules, represented as limit-cycle attractors, see **Figure 1c**.[24,25]
These self-excited oscillations with time-varying properties can manifest in various ways,
encompassing vibrating soliton molecules – where the soliton separation oscillates – as well as
phase-only oscillations.[26-32]

Recent laser experiments have implemented strategies of external periodic perturbation to probe and even induce transitions within soliton molecules, bringing in a close analogy to the approach of molecular spectroscopy. For instance, the resonant excitation of optical soliton molecules within a Kerr-lens mode-locked Ti:sapphire laser probed the intramolecular interaction to reveal its anharmonic feature.[33] The internal oscillatory dynamics of optical soliton molecules have been experimentally synchronized with a modulated signal injected into
an ultrafast fiber laser.[34] These experiments strengthen the hypothesis that the individual solitons composing a multisoliton waveform are often governed by low-dimensional dynamical systems. Since optical solitons break phase symmetry as well as time-translation symmetry, the internal dynamics of a soliton-pair molecule has at least two degrees of freedom, which is enough for self-oscillations but not for chaos. As a matter of fact, there has not been any observation to date of low-dimensional chaos concerning the internal motions within soliton molecules, either in theory or in experiment. The seminal theoretical predictions of chaotic soliton-pair molecules, based on the complex Ginzburg-Landau equation, needed either the addition of an external drive [35] or a soliton shape instability.[25]

The motivation for discovering low-dimensional chaotic soliton molecules is two-fold. It will strengthen the analogy with the molecules of matter, which are indeed subjected to nonlinear and chaotic vibrations.[36] In addition, it will stimulate the development of effective low-dimensional interaction models applicable within a given range of system parameters inside otherwise highly multi-dimensional complex systems, improving the prospects of analysis and control of optical soliton molecules.

Hence, in this article, we demonstrate the excitation and all-optical control of a chaotic soliton-pair molecule within a passively mode locked fiber laser. Such experiment requires an extremely well-resolved recording of the relative soliton separation in real-time, which is achieved with the balanced optical cross-correlation (BOC) technique.[37] Indeed, BOC tracks the pulse separation variations within a vibrating soliton molecule with sub-femtosecond resolution, enabling the observation of its transitions to chaotic dynamics for the first time, to the best of our knowledge. The route through period doubling bifurcations to chaos is explicitly retrieved with a perfect reproducibility. The chaotic dynamics are qualitatively characterized by the direct BOC waveform, its radiofrequency (RF) spectrogram and phase portraits analysis. They are also quantitatively analyzed using the Lyapunov exponent and the correlation dimension analysis. Finally, we achieve an all-optical control of the chaotic dynamics followed by the soliton molecule, by means of a weak external signal injection, without altering the central laser parameters. The experimental observations are qualitatively well reproduced by numerical simulations.
## 2. Experimental results on chaotic soliton molecules

<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic diagram of the real-time characterization of the extension of dissipative optical soliton molecules. (a) Experimental setup of mode-locked fiber laser. LD, laser diode; WDM, wavelength division multiplexer; EDF, Er-doped fiber; OC, Output coupler; Col., collimator; PBS, polarization beam splitter; QWP, quarter wave plate; HWP, half wave plate; ISO, isolator. The soliton molecules with diverse internal dynamics are generated by the mode-locked fiber laser and sent into the BOC measurement. (b) An S-shaped curve is obtained when the delay between the two arms is scanned, which allows to define a linear measurement range (red-dotted line). (c) We could measure three intra-molecular dynamics: stationary states with a fixed separation, vibrating states following a regular oscillation, and chaotic states characterized by unpredictable separation dynamics. They respectively correspond to: fixed points, limit cycles and strange attractors of the effective dynamical system. OSC: oscilloscope. ESA: electric spectrum analyzer.</figcaption></figure>

### 2.1. Monitoring the extension of the soliton molecule in real time

Figure 1 displays the concept and implementation of the characterization of dissipative optical soliton molecules. These molecules are generated within a fiber ring laser cavity, which comprises an erbium-doped fiber for laser emission around the 1.5-micron telecom wavelength, and an effective and controllable ultrafast saturable absorber based on the nonlinear polarization evolution that takes place in the single-mode optical fibers (see the Method Section). Upon the adjustment of the pump power and the intracavity retarding waveplates, the laser can self-generate soliton-pair molecules in various dynamical states. The output soliton-pair molecules are beam-split into two branches and recombined with orthogonal polarization states using a Michelson interferometer, as sketched in Figure 1a. We then detect the real time dynamics of
the variations of the pulse separation within the soliton molecule by means of the BOC
measurement technique. When the overlap between the two arms is scanned, the BOC outputs
an S-shaped voltage signal (see **Figure 1b**). We set the zero crossing of the S-shaped curve as
the working point, around which the voltage excursions of the BOC output are proportional to
the variations of the pulse separation within the soliton molecule, whereas the laser intensity
noise is cancelled out by the balanced photodetection (more details about the BOC
measurement technique are provided in the Supporting Information, section A).

The recording of the BOC voltage signal allows true real time monitoring of the pulse
separation at the sub-femtosecond level without relying on any heavy data post-processing.
Supported with a spectral analysis, such monitoring facilitates a comprehensive exploration of
complex and low-amplitude internal soliton molecule dynamics. In the case of stationary soliton
molecules, the BOC output displays a flat curve with small random fluctuations, whereas for
vibrating molecules, periodic trajectories are observed. We anticipate that a breakdown of the
periodicity of the time-domain BOC signal can be, in case significant excursions of the BOC
signal remain, a signature of chaotic dynamics within the soliton molecule (see Figure 1c).

**2.2. Generating chaotic soliton molecules**

We have developed a controllable way to excite soliton molecules through a cascade of period-
doubling bifurcations leading to chaos. Figure 2a is an extensive recording of the evolution of
the RF spectrogram of the BOC output when the pump strength is gradually increased. As a
starting point, we set a fundamental period-one (P1) oscillating soliton molecule. An example
is shown in the lower part of Figure 2b, featuring quasi-sinusoidal oscillations of the pulse
separation. At a pump driving current I = 509 mA (optical pump power 305 mW), the oscillation
of the BOC output signal corresponds to an intramolecular vibration between the two solitons
having an amplitude of ~80 fs. The oscillation frequency f_s is measured from the corresponding
RF spectrum displayed in Figure 2c, featuring a sharp peak located at 2.36 MHz accompanied
with its 2^nd^ harmonic. The electronic noise floor is also displayed in Figure 2c (gray curve).
Accordingly, the oscillation period corresponds to about 19 cavity roundtrips. We note that the
oscillation period is not exactly an integer multiple of the roundtrip time, the latter phenomenon
being called subharmonic entrainment and requiring a careful search of specific system
parameters.[15,38,39]
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Generation of chaotic soliton molecules in an ultrafast fiber laser. (a) A full-record spectrogram shows the route of period doubling bifurcation to chaotic soliton molecules by increasing the pump strength. (b) Vibration waveform of the pulses separation for P1, P2, P4 and chaotic soliton molecules. (c) RF spectrum of the vibration waveforms displayed in (b). (d) RF spectrum of the laser output for P1, P2 and chaotic soliton molecules. (e) Average optical spectrum of a chaotic state recorded by an optical spectrum analyzer.</figcaption></figure>

When the pump strength is increased, the P1 branch first maintains its stability with a slightly varying fundamental frequency, until a new set of spectral peaks appear at odd multiples of $f_s/2$ in the RF spectrogram (Figure 2a), signaling a period-doubling bifurcation to a new family of period-doubled (P2) states. This bifurcation manifests in the vibration observed in the time domain (Figure 2b) as a breaking of the symmetry of the oscillation pattern, which now repeats every two oscillations. We display the RF spectrum of the P2 state at $I = 523$ mA in Figure 2c. The signal to noise ratio (SNR) of the $f_s/2$ frequency component is 42.6 dB, indicating a high stability. With a further increase of the pump driving current, the next bifurcation to the period-four (P4) family of steady states appears at $I \sim 525$ mA, see the vibration waveform in Figure 2b and its RF spectrum in Figure 2c. Finally, when the pump driving current increases above 526 mA, the sharp spectral peaks disappear and the RF spectrum becomes broadband, suggesting that the internal motion of soliton molecules has become chaotic. The two pulses, still bounded within the soliton molecule, experience irregular oscillations and their separation trajectory is aperiodic, in contrast with the periodic trajectories
of the P1, P2, and P4 families. Here, a clear P1-P2-P4-...-chaos sequence is demonstrated,
following the classical route of period-doubling bifurcation cascade to chaos. However, at a
pump driving current exceeding 532 mA, the soliton molecule becomes once again periodic,
first with limit cycle attractors belonging to a second P4 family, and then undergoes a reverse
period-doubling bifurcation to reach a second P2 family. Finally, another transition to chaos
starts for $I > 538$ mA. In **Figure 2d**, we plot the RF spectrum of the direct photodetection of the
laser output: it yields the cavity repetition frequency $f_c$, where both P1 oscillations (gray curve),
P2 oscillations (blue curve) and chaos (red curve) are considered for comparison. For the P1
state, the fundamental cavity repetition frequency $f_c$ is accompanied by sharp modulation
sidebands, whose separation from $f_c$ indeed matches with the BOC-measured oscillation
frequency of the soliton molecule.

Nevertheless, we emphasize that the relative amplitude of these sidebands remains small,
typically 35 dB or more below the main $f_c$ peak: this indicates a minor effect of the internal
motion of the soliton molecule on the total energy of the soliton molecule. This observation is
important to rule out any significant soliton shape instability contribution to the chaos observed
in our experiments.[25] For the P2 state, a pair of sub-sidebands appears, with a frequency
difference from $f_c$ being half of the initial P1 oscillating frequency. When reaching the chaotic
states, the sharp modulation peaks disappear while the floor around $f_c$ increases by more than
10 dB compared with the P1 state, further demonstrating the aperiodicity of the chaotic state.

As a highlighted feature of the internal soliton molecular motions reported here, from
periodic to chaotic, we emphasize again on the fact that the optical soliton molecule always
maintains its integrity, as a strongly bound soliton-pair system. This is linked to the fact that
the relative timing variations, typically in the range of ~100 fs, remain small compared to the
average soliton separation, which is here about 6.4 ps and matches with the spectral interfringe
of periodicity 1.28 nm shown in the average optical spectrum in **Figure 2e**. Finally, to support
the view that these results reveal widely represented soliton molecule dynamics and are not a
coincidence resulting from overly specific laser parameters, we have reproduced the bifurcation
route to chaos of the internal motion of soliton molecules from a fiber laser cavity featuring a
shorter SMF fiber length, see the Supporting Information, section B. We have also implemented
a numerical simulation to qualitatively account for the observed pulse dynamics, based on the
generalized nonlinear Schrödinger propagation equation, see the Supporting Information,
section C. The numerical results support our experimental findings of a clear route of period
doubling bifurcation to chaos for the internal motion within soliton-pair molecules. In spite of
these major bifurcations experienced by the soliton molecule, we verify that the associated
intracavity energy fluctuations remain small (~1% or below,) see the Supporting Information, section D, confirming the leading role of the internal dynamics on the observed phenomena.

## 2.3. Dynamical characterization of chaotic soliton molecules

Having established a deterministic route to chaotic soliton molecules, we next analyze the statistics of the chaotic trajectories recorded as in Figure 2b, an example of which is shown in **Figure 3a**. **Figure 3b** shows the chaotic orbit in a reconstructed three-dimensional phase space, which is obtained from the BOC output signal of Figure 3a and the delayed BOC signals with delay times $τ$ and $2τ$. The delay time $τ$ is chosen by the mutual information method, which is frequently used in the phase space reconstruction of one-dimensional time series.[40] The geometric structures can be seen from the projections of the phase trajectory on the x-y, x-z and y-z phase planes (gray curves). To quantitatively demonstrate that the reported interacting dynamics within the soliton molecule is chaotic, we extract the largest Lyapunov exponent $λ_L$ from the BOC signals. Lyapunov exponents quantify the divergence rate of nearby attractor trajectories in phase space and are widely used as a criterion for identifying chaos.[41] A necessary condition for chaos is the existence of at least one positive Lyapunov exponent $λ_L$. As shown in **Figure 3c**, the indicator of chaos is reflected by the average initial logarithmic divergence curve over all pairs of neighbors, where an exponential divergence is reasonably linear with the Lyapunov exponent $λ_L ≈ 6.5/μs$. We note that the divergence trajectory features saturation since the chaotic system is bounded in the phase space so that the average divergence cannot exceed the boundary of the strange attractor.

We also analyze the correlation dimension $D_2$ by using the Grassberger-Procaccia (GP) algorithm.[42] The correlation dimension is an efficient diagnostic to distinguish between chaos and random motion. For a random motion, the measured $D_2$ grows linearly with the increasing embedding dimension $m$, whereas in the situation of chaos, $D_2$ saturates to an asymptotic value. **Figure 3d** displays the logarithmic plot of the correlation function $C_m(r)$ as a function of the radius $r$ for increasing $m$ and the corresponding correlation integral slope versus the radius $r$ is shown in **Figure 3e**. The curve manifests a clear saturation of the correlation integral slope, towards an estimated correlation dimension $D_2 ≈ 2.28$. A finite value for the correlation dimension supports the hypothesis of chaotic internal soliton-pair dynamics and gives an estimate of the degree of complexity of the motion of the interacting solitons. A detailed description of the analyzing method involving the delay time $τ$, the embedding dimension $m$, the Lyapunov exponent $λ_L$, and the correlation dimension $D_2$ is presented in the Supporting Information, section E.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Statistical analysis of a chaotic internal soliton molecule motion. (a) BOC output signal. (b) A three-dimensional phase space orbit reconstructed from the BOC output signal and the delayed BOC output signals with delay times $τ$ and $2τ$. The gray curves are projections of the phase portrait on three phase planes. (c) Lyapunov exponent analysis of the BOC output signal. Inset: enlarged picture shows the exponent fitting result. (d) Logarithmic plot of the correlation function $C_m(r)$ versus the radius $r$ for increasing embedded dimension $m$. (e) Slope of the correlation function versus the radius $r$. A clear plateau is observed, with the black dashed curve giving an estimation of the correlation dimension.</figcaption></figure>

We verify the reproducibility of chaotic soliton molecule generation, by using a pump current modulated at 100 Hz with a wide modulation amplitude such that the soliton molecules experience a complete circular sequence of bifurcations (P1 - chaos - P1), as observed in Figure 2a for a driving current ranging from 500 to 530 mA. The resulting real-time BOC output signals of the pulse separation are shown in Figure 4a, with Fourier spectral analysis in Figure 4b: the recording demonstrates a perfect reproducibility of the bifurcation sequences leading to or departing from chaos. During the first half modulation period, a switching occurs from stable P1 to P2 states. A transient P4 state is then observed before the onset of chaos. Subsequently, the system returns to the P1 state through the reverse sequence of bifurcations, initiated by the reduction of the intracavity energy. This dynamical switching between regular and chaotic oscillations is also illustrated in the Supporting Information videos 1&2, where we videotaped this reversible process, showing the real-time evolution of the waveforms and the RF
spectrograms after setting a longer modulation period suitable for a clear view in real time. Whereas the bifurcation sequence is repeatedly observed, the details of intra-molecular chaotic dynamics can vary substantially over consecutive experimental runs, with another example mirroring the analysis of Figure 3 shown in the Supporting Information figure S5.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Real-time BOC output signals of the variations of pulse separation within a soliton molecule during the pump current modulation at 100 Hz. (b) The two-dimensional spectra extracted from windowed Fourier transforms of the BOC output signals in (a).</figcaption></figure>

**2.4. All optical control of chaotic soliton molecules**

Chaotic dynamics is prone to exhibit considerable reaction to additional perturbations. In chaos theory, it was shown that chaos can be suppressed to a limited cycle attractor by non-feedback methods such as the application of external periodic forces.[43] The elimination of chaos is most effective when the drive frequency (frequency of the external force) is close to the eigenfrequency of self-excited oscillations that exist in the vicinity (in the parameter space) of the chaotic motion. In our experiment, the free-oscillation frequency of the vibrating soliton molecule varies within the range of 1-5 MHz, so we choose a drive frequency $f_d = 2$ MHz. We here demonstrate that the chaotic dynamics involving bound solitons in ultrafast fiber lasers can be suppressed by external cavity injection, without affecting the original laser parameters and the structure of the individual solitons. First, we experimentally verify that the soliton molecule remains a robust entity with respect to the injection signal, including within the chaotic region. Therefore, the soliton molecule presents an opportunity of external control. As shown in **Figure**
**5a**, we implement the external control of chaos via the injection of a modulated continuous wave (CW), which is an efficient method for pulse dynamics modulation. In the optical spectral domain, the CW is selected at 1530 nm, far enough from the soliton molecule spectrum to avoid a deleterious influence on the mode locking regime within a suitable range of injection power. The CW is amplified by an Er-doped fiber amplifier (EDFA), whose output power is adjustable. An electro optic modulator with a bandwidth exceeding 1 GHz driven by a function generator imprints a sinusoidal intensity modulation on the injection signal. The latter is injected into the cavity via a 10% coupler in the counterpropagating lasing direction, thus directly influencing the gain saturation level of the fiber laser with a bandwidth up to several MHz.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** All-optical control of chaotic soliton molecules. (a) Conceptual illustration of the all-optical control of chaotic soliton molecules. (b) RF spectrum of BOC output signals at three injection strengths. (c) A full-record spectrogram showing the control of chaotic soliton molecules, which retrieve a single-periodic oscillation (P1) upon the increase of the injection signal. (d) Vibration waveform of the inter-soliton separation corresponding to the cases in (b).</figcaption></figure>

We experimentally prepare a chaotic soliton molecule whose dynamics of inter-pulse separation is manifestly broadband, as attested by the RF spectrum of the BOC output signal (black curve in **Figure 5b**) at 550 mA pump strength. Then, we control this state by increasing the injected CW laser intensity, measured as the driving power of the external EDFA. At a driving power of 2.8 mW, the chaotic motion becomes intermittent, interspersed by P2-like
oscillations (blue curve in **Figure 5d**) with corresponding spectral peaks at half-integer multiples of $f_s$ located at about 2.5 MHz on top of a broad continuum (blue curve in Figure 5b). Robust P1 oscillations are obtained at the current drive of 9.2 mW, as attested by both the RF spectrum and the regular vibration waveform (red curves in Figure 5b & 5d). The complete and deterministic route of chaos control through period-doubling reversals via external cavity injection is shown in **Figure 5c**. Note that excessive injection power ($> 9.8$ mW) will result in the disruption of the soliton molecule. The frequency components of the injection signal and beat notes are also evidenced by the RF spectrum in Figure 5b, confirming that the external injection is efficiently coupled to the chaotic motion. Overall, the route of external injection for chaos control evolves in a highly reproducible and reversible manner, so that the internal motion of soliton molecules reverts to chaos if we weaken the injection strength. For a clearer visual presentation of our results, we also videotaped this fully efficient switching of chaos control, see Supporting Information section S3 and video3.

### 3. Discussion and conclusion

Whereas oscillating optical soliton molecules have been the subject of several seminal investigations, this article presents the first experimental demonstration of the route for excitation and all optical control of a chaotic soliton molecule in optical resonators and lasers. We emphasize on the fact that the chaos investigated here refers to the internal dynamics of the soliton molecule and is assessed to be low-dimensional, leaving the soliton molecule energy practically unchanged. This can be understood physically by the chosen situation where the variations of the inter-pulse separation are of the order of one percent of that separation, therefore perturbing only slightly the overall shape of the soliton molecule. Therefore, to observe such internal chaos requires a characterization method endowed with an exquisite sub-femtosecond sensitivity and true real-time monitoring. The BOC method fully qualifies, doing so without relying on any heavy post-processing of data. This allows to record experimentally the vibration waveform with the necessary high precision to observe and analyze the experimental internal chaos. We therefore demonstrated that the chaotic dynamics of bound pulses within a soliton molecule can be accessed and manipulated in an unambiguous way. The BOC method appears particularly suited in the investigation of the complex nonlinear dynamics of soliton molecules in optical oscillators, as much as it is widely used for timing jitter characterization in ultrafast lasers and ring resonators.[44,45] We anticipate that it will be a useful asset in precise investigation of nonlinear dynamics effects such as stochastic resonance[46] and subharmonic entrainment [39] of solitons molecules.
Following the present investigation, the intramolecular dynamics of optical soliton molecules is gaining diversity as it further strengthens the analogy with matter-like molecules. Indeed, our study constitutes a strong confirmation of the hypothesis that a wide range of soliton molecule oscillations can be viewed as low-dimensional dynamical systems that feature limit cycle attractors as well as strange attractors. This hypothesis stays reasonably valid if the amplitude of the intramolecular motion remains sufficiently small compared to the average separation between solitons.

The fast error-free switching between ordered and chaotic soliton molecules enabled by pump current sweeping and external injection highlights the potential prospects of all-optical logic gates and chaotic communication using soliton molecules. Beyond the fundamental case of soliton-pair molecules, there is a large variety of molecular complexes, macro-molecules and soliton crystals that could be similarly investigated. It is reasonable to assume that the chaotic dynamics observed in our experiment could be extended to three-soliton interactions in lasers,[47] and even more complicated multi-body dynamics such as interactions within soliton molecular complexes and supramolecular structures in parallel optical-soliton reactors.[48] More generally, the findings of chaotic interactions among dissipative solitons will be of significant interest in the context of analogous nonlinear systems, such as Bose-Einstein condensates and hydrodynamics.

**4. Experimental methods**

*Laser setup.* The chaotic soliton molecules under investigation are generated from an Er-doped fiber laser, mode locked by the nonlinear polarization evolution (NPE) technique. NPE provides a virtual quasi-instantaneous saturable absorber effect, whose transfer function can be widely adjusted by altering the orientations of the intracavity wave plates, leading to multifarious mode-locked states. The cavity incorporates two types of fibers: a 0.5-m long EDF (Liekke Er 110-4/125) serves as the gain medium and the others are standard single-mode fiber (SMF 28). By changing the length of the SMF, we obtained two mode-locking states with different fundamental repetition frequency and net chromatic dispersion. The total length of the laser cavity of 4.6 (4.1) m yields a fundamental repetition frequency of 44.8 (50.6) MHz, corresponding to a roundtrip time of 22.3 (19.76) ns. The laser has a net anomalous dispersion ~ -0.07(-0.06) ps² at 1.55-µm wavelength. In both states, the mode locking threshold is at a pump driving current of about 380 mA (optical pump power ~ 230 mW) and the threshold of multi-pulsing (MPI) is at about 455 mA (optical power ~ 270 mW). Soliton molecular chaotic
dynamics are observed in both mode-locking states, proving that the observation is not a coincidence resulting from specific laser parameters.

*Balanced optical cross-correlator.* A single-crystal balanced optical cross-correlation (BOC) scheme has been implemented to detect the internal soliton molecular motion. The BOC setup consists of a 4-mm type-II phase-matched periodically poled KTiOPO₄ (PPKTP) crystal (facet 1 coating: HT@ 1550 nm & 775 nm, facet 2 coating: HT@ 775 nm & HR@ 1550 nm), a dichroic mirror (HT@1550 nm & HR@ 775 nm), a focusing lens ($f=30$ mm) and a balanced photodetector (Thorlabs PDB420A). The output of the BOC is detected by an oscilloscope (Agilent infiniium) for monitoring of intramolecular pulse separation in real time and a RF spectrum analyzer (Rigol DSA815) for spectrum analysis of chaotic dynamics.

*Pump current modulation.* The pump laser diode (LD) is driven by a commercial laser diode controller (Thorlabs CLD1015), which has an external RF input port for pump current modulation. A 100 Hz sinusoidal modulation generated by a function generator (Rigol DG5102) is used to modulate the laser. The modulation amplitude can be adjusted by varying RF voltage from the function generator output. The modulation bandwidth can be as high as 10 MHz, while, in Er-fiber lasers, it is limited to few kilohertz due to the low pass filtering effect caused by the damped gain relaxation oscillations. Ref. 51 used this scheme for intramolecular pulse separation manipulation. Here, this technique allows us to reproduce the bifurcation sequences leading to or departing from chaos.

*All-optical control of chaotic soliton molecules.* The external optical injection is based on an all-fiber configuration. The 1530 nm CW laser is generated by a commercial low noise tunable laser (Santec TSL-550). The output is power amplified by a home-made all polarization maintaining EDFA. A fiber-coupled high speed electro optic modulator (Conquer KG-AM-15-2.5G) is used to sinusoidally modulate the intensity of the CW and the MHz modulation signal is from a function generator (Rigol DG5102).

**Supporting Information**

Supporting Information is available from the Wiley Online Library or from the author.

**Acknowledgements**
This work was supported by the National Natural Science Foundation of China (Grant 61975144, 61827821); Ph.G. acknowledges support from PIA3-ISITE-BFC (ANR-15-IDEX-0003); O. Gat acknowledges support from the Israel Science Foundation (ISF) (Grant No. 2403/20). Y.S and D.Z. contributed equally to this work.

Received: ((will be filled in by the editorial staff))

Revised: ((will be filled in by the editorial staff))

Published online: ((will be filled in by the editorial staff))

References

[1] Poincaré, Henri. Sur le problème des trois corps et les équations de la dynamique. Divergence des séries de M. Lindstedt. *Acta Mathematica*. **13**, A3–A270 (1890).

[2] Lorenz, E. Deterministic non-periodic flow. *Journal of the Atmospheric Sciences* **20**, 130–141 (1963).

[3] Strogatz, S. H. *Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering*. CRC press, 2018.

[4] Argyris, A. et al. Chaos-based communications at high bit rates using commercial fibre-optic links. *Nature* **438**, 343–346 (2005).

[5] Spitz, O. et al. Private communication with quantum cascade laser photonic chaos. *Nat. Commun.* **12**, 3327 (2021).

[6] Sciamanna, M. & Shore, K. A. Physics and applications of laser diode chaos. *Nat. Photonics* **9**, 151–162 (2015).

[7] Jumpertz, L., Schires, K., Carras, M., Sciamanna, M. & Grillot, F. Chaotic light at mid-infrared wavelength. *Light Sci. Appl.* **5**, e16088 (2016).

[8] Deng, Y. et al. Mid-infrared hyperchaos of interband cascade lasers. *Light Sci. Appl.* **11**, 1–10 (2022).

[9] Grelu, P. & Akhmediev, N. Dissipative solitons for mode-locked lasers. *Nat. Photonics* **6**, 84–92 (2012).

[10] Sucha, G., Bolton, S. R., Weiss, S. & Chemla, D. S. Period doubling and quasi-periodicity in additive-pulse mode-locked lasers. *Opt. Lett.* **20**, 1794–1976 (1995).

[11] Soto-Crespo, J.M., Grapinet, M., Grelu, P. & Akhmediev, N. Bifurcations and multiple-period soliton pulsations in a passively mode-locked fiber laser. *Phys. Rev. E* **70**, 066612 (2004).
[12] Cundiff, S. T., Soto-Crespo, J. M. & Akhmediev, N. Experimental evidence for soliton explosions. *Phys. Rev. Lett.* **88**, 073903 (2002).

[13] Horowitz, M., Barad, Y. & Silberberg, Y. Noiselike pulses with a broadband spectrum generated from an erbium-doped fiber laser. *Opt. Lett.* **22**, 799-801 (1997).

[14] Runge, A. F., Aguergaray, C., Broderick, N. G. & Erkintalo, M. Coherence and shot-to-shot spectral fluctuations in noise-like ultrafast fiber lasers. *Opt. Lett.* **38**, 4327-4330 (2013).

[15] Soto-Crespo J. M. & Akhmediev, N. Soliton as strange attractor: Nonlinear synchronization and chaos. *Phys. Rev. Lett.* **95**, 024101 (2005).

[16] Wang, Z. Q., Nithyanandan, K., Coillet, A., Tchofo-Dinda, P., & Grelu, P. Buildup of incoherent dissipative solitons in ultrafast lasers. *Phys. Rev. Research 2*, 013101 (2020).

[17] Lucas, E., Karpov, M., Guo, H., Gorodetsky, M. & Kippenberg, T. Breathing dissipative solitons in optical microresonators. *Nat. Commun.* **8**, 736 (2017).

[18] L. B. A. Mélo, G. F. R. Palacios, P. V. Carelli, L. H. Acioli, J. R. Rios Leite, and M. H. G. de Miranda. Deterministic chaos in an ytterbium-doped mode-locked fiber laser. *Opt. Express* **26**, 13686–13692 (2018).

[19] Wang, W., Xian, T., Zhang, M. & Zhan, L. Real-time revealing phase transition of bifurcation evolution in modelocked lasers. *Opt. Lett.* **47**, 1234-1237 (2022).

[20] Shlizerman, E., Ding, E., Williams, M. O. & Kutz, J.N. Characterizing and suppressing multi-pulsing instabilities in mode-locked lasers. *Proc. SPIE* **7933**, 79331L (2011).

[21] Lecaplain, C., Grelu, P., Soto-Crespo, J. M. & Akhmediev, N. Dissipative rogue waves generated by chaotic pulse bunching in a mode-locked laser. *Phys. Rev. Lett.* **108**, 233901 (2012).

[22] Chouli, S. & Grelu, P. Soliton rains in a fiber laser: an experimental study. *Phys. Rev. A* **81**, 063829 (2010).

[23] Meng, F., Lapre, C., Billet, C. et al. Intracavity incoherent supercontinuum dynamics and rogue waves in a broadband dissipative soliton laser. *Nat. Commun.* **12**, 5567 (2021).

[24] Grapinet, M., & Grelu, P. Vibrating soliton pairs in a mode-locked laser cavity. *Opt. Lett.* **31**, 2115-2117 (2006).

[25] Soto-Crespo, J. M., Grelu, P., Akhmediev, N. & Devine, N. Soliton complexes in dissipative systems: Vibrating, shaking, and mixed soliton pairs. *Phys. Rev. E* **75**, 016613 (2007).

[26] Ortaç, B. et al. Observation of soliton molecules with independently evolving phase in a mode-locked fiber laser. *Opt. Lett.* **35**, 1578-1580 (2010).
[27] Krupa, K., Nithyanandan, K., Andral, U., Tchofo-Dinda, P. & Grelu, P. Real-time observation of internal motion within ultrafast dissipative optical soliton molecules. *Phys. Rev. Lett.* **118**, 243901 (2017).

[28] Herink, G., Kurtz, F., Jalali, B., Solli, D. & Ropers, C. Real-time spectral interferometry probes the internal dynamics of femtosecond soliton molecules. *Science* **356**, 50–54 (2017).

[29] Zhou, Y., Ren, Y.-X., Shi, J., Mao, H. & Wong, K. K. Y. Buildup and dissociation dynamics of dissipative optical soliton molecules. *Optica* **7**, 965–972 (2020).

[30] J. Peng, Z. Zhao, S. Boscolo, C. Finot, S. Sugavanam, D. V. Churkin, and H. Zeng, Breather Molecular Complexes in a Passively Mode-Locked Fiber Laser. *Laser & Photonics Rev.* **15**, 2000132 (2021).

[31] Nimmesgern, L., Beckh, C., Kempf, H., Leitenstorfer, A. & Herink, G. Soliton molecules in femtosecond fiber lasers: universal binding mechanism and direct electronic control. *Optica* **8**, 1334–1339 (2021).

[32] Hamdi, S., Coillet, A., Cluzel, B., Grelu, P. & Colman P. Superlocalization reveals long-range synchronization of vibrating soliton molecules. *Phys. Rev. Lett.* **128**, 213902 (2022).

[33] Kurtz, F., Ropers, C. & Herink, G. Resonant excitation and all-optical switching of femtosecond soliton molecules. *Nat. Photonics* **14**, 9-13 (2020).

[34] Zou, D., Song, Y., Hu, M., Gat., O. & Grelu P. Synchronization of the internal dynamics of optical soliton molecules. *Optica* **9**, 1307 (2022).

[35] Turaev, D., Vladimirov, A.G. & Zelik S. Chaotic bound state of localized structures in the complex Ginzburg-Landau equation. *Phys. Rev. E* **75**, 045601 (2007).

[36] Wu, G. *Nonlinearity and Chaos in Molecular Vibrations*. Elsevier Science (2005).

[37] Kim, J., Chen, J., Cox, J., & Kärtner, F. X. Attosecond-resolution timing jitter characterization of free-running mode-locked lasers. *Opt. Lett.* **32**, 3519–3521 (2007).

[38] Cole, D. C., & Papp, S. B. Subharmonic entrainment of Kerr breather solitons. *Phys. Rev. Lett.* **123**, 173904 (2009).

[39] Xian, T., Zhan, L., Wang, W., & Zhang, W. Subharmonic entrainment breather solitons in ultrafast lasers. *Phys. Rev. Lett.* **125**, 163901 (2020).

[40] Fraser, A. M., & Swinney, H. L. Independent coordinates for strange attractors from mutual information. *Phys. Rev. A* **33**, 1134 (1986).

[41] Rosenstein, M. T., Collins, J. J. & De Luca, C. J. A practical method for calculating largest Lyapunov exponents from small data sets. *Physica D: Nonlinear Phenomena* **65**, 117–134 (1993).
[42] Grassberger, P. & Procaccia, I. Measuring the strangeness of strange attractors. *The theory of chaotic attractors*. (pp. 170-189) Springer, New York, (2004).

[43] Boccaletti, S., Grebogi, C., Lai, Y. C., Mancini, H. & Maza, D. The control of chaos: theory and applications. *Phys. Rep.* **329**, 103–197 (2000).

[44] Benedick, A. J., Fujimoto, J. G. & Kärtner, F. X. Optical flywheels with attosecond jitter. *Nat. Photon.* **6**, 97–100 (2012).

[45] Bao, C. et al. Quantum diffusion of microcavity solitons. *Nat. Phys.* **17**, 1–5 (2021).

[46] Monifi, F. et al. Optomechanically induced stochastic resonance and chaos transfer between optical fields. *Nat. Photon.* **10**, 399–405 (2016).

[47] Xin, F., et al. Evidence of Chaotic Dynamics in Three-Soliton Collisions. *Phys. Rev. Lett.* **127**, 133901 (2021).

[48] He, W. et al. Synthesis and dissociation of soliton molecules in parallel optical-soliton reactors. *Light Sci. Appl.* **10**, 1–15 (2021).
Table of contents entry:

Can the analogy between optical soliton molecules and matter molecules be further developed?
For instance, do optical soliton molecules allow low-dimensional chaos? An ultrafast fiber laser
experiment with advanced real-time characterization yields a positive answer. It demonstrates
sequences of bifurcations leading to chaotic relative motions of the bound optical solitons.
Furthermore, through signal injection, these internal motions within the soliton molecule are
resynchronized to the external oscillator, opening new possibilities in terms of modeling,
control, and application of these sophisticated optical pulse waveforms.

Youjian Song, Defeng Zou, Omri Gat,\* Minglie Hu,\* and Philippe Grelu\*

Chaotic internal dynamics of dissipative optical soliton molecules

ToC tentative Figure:

<figure><img src="image_7.png" /></figure>
# Supporting Information

## Chaotic internal dynamics of dissipative optical soliton molecules

Youjian Song, Defeng Zou, Omri Gat,\* Minglie Hu,\* and Philippe Grelu\*

In the following, we provide the technical details concerning: (A) the balanced optical cross-correlation (BOC) method utilized in the experiment; (B) an additional experimental result demonstrating perfect reproducibility of the chaotic soliton molecules by fast pump current sweeping; (C) detailed descriptions of the criteria for chaos and another exampleed statistics characteristic of the BOC signal; (D) an additional experimental result demonstrating the universality of the route to chaotic soliton molecules; (E) experimental results of the pulse energy evolution during the route to chaos; (F) the ultrafast fiber laser model used to numerically confirm the existence of the chaotic attractor of interacting solitons that we observed in experiments.

### A. Calibration of the BOC signal

As shown in **Figure S1a**, the fiber laser used here is mode-locked by nonlinear polarization evolution (NPE), which is facilitated through the inclusion of a combination of three wave plates and a free space isolator. A short length of erbium-doped fiber (0.5-m Liekki Er 110-4/125) serves as the gain medium and the other fibers are standard single mode fiber (SMF). The laser yields a net anomalous dispersion ~ -0.07(-0.06) ps² at 1.55-µm wavelength and stretch pulses of conventional solitons.

The output soliton molecule is then split into two branches by an unequal-arm Michelson interferometer (MI). This allows partially overlapping the leading and the trailing solitons with orthogonal polarizations. As shown in **Figure S1b**, we implement a single crystal BOC configuration, by using a 4-mm type-II phase-matched periodically poled KTiOPO₄ (PPKTP) crystal for sum frequency generation in forward/reverse direction and two dichroic mirrors (DM1 transmits 1550 nm laser and reflects the 775 nm SFG signal, DM2 transmits 775 nm SFG signal and reflects 1550 nm laser) for separating the two SFG signals in two directions. We use a balanced photodetection to detect the SFG signals. By sweeping M1 in the MI, a S-shaped output voltage versus the relative pulse timing is observed on an oscilloscope [see **Figure S1c**]. The zero crossing (represented by a green point) is found when the averaged pulse
separation at the input of the BOC equals exactly $t_g$, which represents the group delay between the orthogonal polarizations in the PPKTP crystal. A linear fit around the zero-crossing region gives a slope about 2.21 mV/fs, which serves as the timing error discriminator.

<figure><img src="image_8.png" /><figcaption>**Figure S1.** Principle to detect variations of intra-molecular pulse separation based on BOC method. (a) Experimental setup of mode-locked fiber laser. LD, laser diode; WDM, wavelength division multiplexer; EDF, Er-doped fiber; OC, Output coupler; Col., collimator; PBS, polarization beam splitter; QWP, quarter wave plate; HWP, half wave plate; ISO, isolator. (b) Principle of optical cross-correlation method. M: silver mirror; DM1,2: dichroic mirrors; PPKTP, periodically poled KTiOPO₄. (c) The S curve represents the BOC output signal obtained by sweeping Mirror in the Michelson interferometer.</figcaption></figure>

## B. Reproducibility of chaotic soliton molecules generation through fast pump current modulation

In this section, we use fast pump power modulation to verify the reproducibility of chaotic soliton molecules. This electronic control method is similar to that, which implement a laser perturbation via current modulation of the semiconductor diode-laser driver in Reference S4. The difference is that we do so to realize the reliable precipitate and recurrence of the chaotic internal dynamics of optical soliton molecules. The bandwidth of the semiconductor diode-laser
is about several kHz, beyond which the modulation will be filtered by the gain relaxation oscillations. We implement the sinusoidal modulation to force the pump current into the range of period doubling bifurcation to chaos exactly. The modulation period is 10 ms and the resultant real-time BOC output signals of the pulse separation are shown in **Figure S2**.

<figure><img src="image_9.png" /><figcaption>**Figure S2.** (a) Real-time BOC output signals of the variations of pulse separation within a soliton molecule during the fast pump power modulation. The modulation period is 10 ms. (b) The two-dimensional spectra extracted from windowed Fourier transforms of the BOC output signals in (a).</figcaption></figure>

As shown in **Figure S2a**, multiple consecutive cycles of BOC signals (~ 40 thousand oscillations of the intra-molecular pulses separation and ~2 million laser cavity roundtrips) show perfect reproducibility during the fast pump sweeping. Further insight into the detailed bound-state response is provided by windowed Fourier transforms of BOC output signals. Specifically, the set of Fourier spectra along the frequency sweep forms a two-dimensional map (**Figure S2b**), representing the complete fast response of the soliton molecules. During the first half modulation period, switching occur from stable P1 state to P2 state. A transient P4 state is also observed before the chaos. Then, the chaos return back to P1 state through the reverse progress, initiated by the global change of intracavity energy in the laser cavity. Notably, the smooth processes instead of overnight processes are observed for all bifurcation processes. This dynamically switching between regular and chaotic oscillation is also apparent from the Supporting Information Video 1, in which we videotaped this reversible process by real-time
evolution of RF spectrograms. The modulation period in the video was adjusted to be 2 seconds for clear viewing.

### C. Phase space reconstruction, Lyapunov exponents and Correlation dimension

*Phase space reconstruction.* For experimental signals with noise, delay time $τ$ and embedding dimension $m$ for phase space reconstruction should be carefully selected, avoiding that reconstructed trajectories in the phase space are indistinguishable or independent.[⁵⁵] In our manuscript, the optimal value of $τ$ is recommended as the value of the first local minima of mutual information, as described in Reference S6. The recommended value of embedding dimension $m$ is estimated as the value where the percentage of false nearest neighbors falls nearly close to zero (below 1/1000 here).[⁵⁷] The parameters are selected as $τ = 70$ ns, $m = 3$ with time span of 5 µs and sampling rate of 1 ns for Fig. 3a in the main text.

*Lyapunov exponents.* Lyapunov exponents quantify the presence of chaos in dynamical systems, by examining the property of sensitivity to initial conditions. At least one positive Lyapunov exponents exists for the chaotic systems. Practically, the chaotic component could be confirmed by the positive largest Lyapunov exponents $λ_L$. There are several algorithms developed for the calculation of the $λ_L$ from experimental scalar time series (such as BOC output signals in our experiment), during which algorithm of Rosenstein et al. (R-algorithm) have been widely used.[⁵⁸] For the reconstructed phase space of vector times series, the basic idea of R-algorithm is that the $j^{th}$ pair of nearest neighbors diverge approximately at a rate given by $λ_L$

$$d_j(i) \approx d_0 \cdot e^{\lambda_L(i\Delta t)} \qquad (S1)$$

where $d_j(i)$ represents the mean Euclidean distance between neighboring trajectories in reconstructed phase space after $i$ discrete-time steps, i.e., $i\Delta t$ seconds. $d_0$ is the initial distance. Taking the logarithm of both sides of Equation S1, one obtains

$$\ln d_j(i) \approx \ln d_0 + \lambda_L i \Delta t \qquad (S2)$$

Therefore, the largest Lyapunov exponents $λ_L$ can be calculated by

$$y(i) = \frac{1}{\Delta t} \langle \ln d_j(i) \rangle \qquad (S3)$$

where $\langle \rangle$ denotes the average over all values of $j$. One could calculate the value of $λ_L$ by fitting the slope of the linear region of the average logarithmic divergence, as done in Figure 3c of the main article text. We verify the accuracy of the used program by Lorentz attractor using the
same parameters in Ref. S8. The $\lambda_L$ value is calculated to be 1.496, consistent with the ideal value of 1.5.

*Correlation dimension.* Correlation dimension ($D_2$) method is a popular method to evaluate the ‘strangeness’ of the strange attractor of chaotic system. According to Grassberger-Procaccia algorithm,[S9] the correlation integral $C_m(r)$ is defined as the probability of two points in the reconstructed phase space separated by a distance $\rho$ less than or equal to $r$:

$$C_m(r) = \lim_{N \to \infty} \frac{1}{N^2} \sum_{i=1}^{N} \sum_{j=1}^{N} \Theta(r - \rho) \quad (S4)$$

where $\Theta$ is the Heaviside step function, $N$ represents the number of points available, $i$ and $j$ are coordinates of the points. $C_m(r)$ behaves as a power of $r$ for small $r$:

$$C_m(r) \propto r^{D_2(m)} \quad (S5)$$

Taking the logarithm of both sides of Eq. S5, one obtains:

$$D_2(m) \approx \lim_{r \to 0} \frac{\ln C_m(r)}{\ln(r)} \quad (S6)$$

In practice, this value is estimated by the slope of the regression region of the logarithmic plot. The value tend to saturate on the finite value of $D_2$ with the increasing embedded dimension $m$ for the chaotic signals and grows linearly for random motion. We calculate $D_2$ of the Lorentz attractor to be 2.06 (Figure S3a) and infinite value for white noise (Figure S3b), consistent with Reference S9.

<figure><img src="image_10.png" /><figcaption>**Figure S3.** (a) Logarithmic plot of $\ln C(r)$ and $\ln(r)$ of Lorentz attractor; (b) Logarithmic plot of $\ln C(r)$ and $\ln(r)$ of white noise.</figcaption></figure>

## D. Other period doubling bifurcation to chaos observed after laser cavity alteration
Another evidence of period doubling bifurcation to chaos of the chaotic soliton molecules obtained at cavity repetition frequency ~51 MHz is shown in **Figure S4**. The adjustment of repetition rate is realized by increasing or decreasing the length of the single-mode fiber.

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

**Figure S4.** (a) Vibration waveform of the BOC output signals of P1, P2, P4 and chaotic soliton molecules. (b) RF spectrum of the vibration waveform as in (a). (c) RF spectrum of the cavity fundamental repetition frequency. (d) Optical spectrum of chaotic states recorded by optical spectrum analyzer (Left: logarithmic coordinates; Right: linear coordinates).

E. Influence of the internal soliton motions on the soliton molecule energy

To investigate the pulse energy evolution during the period doubling bifurcation to chaos, we record the output pulse sequences with a photodiode that is connected to a digital oscilloscope with 1 GHz analog bandwidth. The changes in the intensity of the pulse sequences reflect the pulse energy change during each cavity roundtrip. For soliton molecules with P1 oscillating internal dynamics, periodic modulation of the amplitude of the pulse sequences is observed in **Figure S5a**. A 5 MHz low-pass filter is used to only extract the envelope of pulse train and the result is shown by the black curve. Quasi-sinusoidal oscillation with period of ~19 roundtrips is clearly observed, which matches the results of BOC signals in Figure 2b in the main text. For P2 oscillation state, the period doubling bifurcation of the pulse energy is observed. Finally, the periodic modulation turns into strong and irregular fluctuations, indicating the existence of the strange attractor of the pulse energy of chaotic soliton molecules.
<figure><img src="image_12.png" /></figure>

**Figure S5.** Pulse energy and output power during the bifurcation process. (a-c) pulse sequences (blue curve) and filtered energy fluctuation (black curve) of the P1, P2 and chaotic oscillations, respectively. (d) Output power of chaotic soliton molecules in 60 minutes, inset is the corresponding histogram of the amplitude of 3600 pulses.

We emphasize that although internal interactions exhibit complex and unpredictable behavior within chaotic soliton molecules, they show great stability as an ensemble circulating in the cavity because the oscillation amplitude remains negligible compared with the amplitude of the pulses train, as shown in **Figure S5 a,b &c**. The robustness of the chaotic soliton molecules could also be witnessed by the other three aspects. At first, a record of the 60-minutes output power of chaotic states is shown in **Figure S5d**, in which the corresponding rms power fluctuation is calculated as low as 0.094%. The histogram of the amplitude of 3600 pulses is shown in inset in Figure S5d, a normal distribution in a narrow range of the pulse intensities is obtained. Besides, the stability of the mode-locked state of chaotic soliton molecules could be confirmed by the high SNR (>> 70 dB) of the repetition frequency $f_c$ in Figure 2d in the main text. Thirdly, the robustness of the chaotic soliton molecules is also evidenced by multiple consecutive cycles tracing of the oscillation waveform of the pulse separation by fast pump power modulation shown in Section B, Supplementary materials. The chaotic state could be
reliably translated into different stability regimes with perfect reproducibility and be maintained for a long time without a special protection device.

## F. Numerical simulation of chaotic soliton molecules

In this section, we conduct numerical simulations of the ultrafast laser dynamics to qualitatively confirm our experimental observations and ascertain the general manifestation of chaotic soliton molecules. These simulations are based on a generalized nonlinear Schrödinger equation (GNLSE), which includes dissipative (gain/loss) terms. The model also includes a saturable absorber (SA) to promote mode locking. This model has been widely used to reproduce complex pulsed laser dynamics, such as soliton pulsations, soliton molecule vibrations and incoherent dissipative solitons:

$$ \frac{\partial A}{\partial z} + \frac{i}{2} \left( \beta_2 + ig \frac{1}{\Omega_g^2} \right) \frac{\partial^2 A}{\partial \tau^2} = \frac{g}{2} A + \frac{\beta_3}{6} \frac{\partial^3 A}{\partial \tau^3} + i\gamma |A|^2 A \quad (S7) $$

where A is the field envelope, Z the propagation coordinate and τ the time in the frame co-moving with the pulses. The fiber dispersion and effective nonlinearity parameters are selected from Corning SMF-28 and Liekki Er 110-4 for simulation. g = g₀/(1+∫|A|²dτ/Eₐ₉₅ₐ₉₅+(ω-ω₀)/Δω₂) is the EDF gain, where g₀ is the small-signal gain, Eₐ₉₅ₐ₉₅ is the saturation energy, ω is the optical angular frequency, ω₀ is the central angular frequency. The saturable gain g has a FWHM bandwidth Δω of 50 nm. The SA is modeled by a power dependent transmittance function: T = 1 - αₘ₅-q₀(1+P/Pₐ₉₅ₐ₉₅)⁻¹, where αₘ₅ is the non-saturable loss, P = |A|² is instantaneous power, Pₐ₉₅ₐ₉₅ is saturation power. The SA is modeled with 8% non-saturable loss, 30% modulation depth and 50 W saturation power. Guided by the experiment, we adjust the energy saturation coefficient Eₐ₉₅ₐ₉₅ of the modeled gain fiber to seek irregular dynamics of interacting solitons. Basically, Eₐ₉₅ₐ₉₅ will be proportional to the pumping power well above the laser threshold.

**Figure S6** illustrates typical examples of soliton molecule dynamics found numerically, revealing the route from period-doubling bifurcations to chaotic soliton molecules. **Figure S6 a & b** show a P1 oscillation dynamics at Eₐ₉₅ₐ₉₅ = 0.63 nJ. From the Fourier-transform of the optical spectra displayed over consecutive cavity roundtrips in **Figure S6a**, an oscillating soliton temporal separation of period ~25 round trips is obtained and displayed in **Figure S6b**. The observed amplitude of the separation oscillation is about 0.05 ps. The soliton molecule energy is also plotted in Figure S6a, showing the same periodicity. Starting from this P1 state, increasing the Eₐ₉₅ₐ₉₅ to 0.69 nJ results in the transition to P2 states. The period doubling characteristics are observed in the optical spectrum, the soliton molecule energy, and the pulse separation, shown in **Figure S6c&d**. At Eₐ₉₅ₐ₉₅ = 0.71 nJ, the chaotic dynamics dominates the
interaction, see **Figure S6e**. We plot the phase portrait of the chaotic state, a three-dimensional plane spanned by pulse separations in **Figure S6f**. The complex geometrical structure of the strange attractor is obtained, revealing a strong comparison with limit cycles attractors (P1 and P2 states). The bifurcation diagrams of both the soliton molecule energy local maxima and pulse separation local maxima are shown in **Figure S6g**. It clearly shows the route of P1-P2-P4-...-chaos.

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

**Figure S6.** Numerical simulations. (a-b) Dynamics of a P1 oscillation state at $E_{sat} = 0.63$ nJ. (a) Spectral evolution and pulse energy (inset) over 150 consecutive round trips. (b) Corresponding evolution of the inter-soliton temporal separation. (c-d) Dynamics of a P2 oscillation state at $E_{sat} = 0.69$ nJ. (c) Spectral evolution and pulse energy (inset) over 150 consecutive round trips. (d) Corresponding evolution of the inter-soliton separation. (e) Spectral evolution of the chaotic soliton molecule and pulse energy (inset) over 1500 consecutive round trips at $E_{sat} = 0.71$ nJ. (f) Phase portraits of chaotic states. (g) Bifurcation diagram of the soliton molecule energy local maxima and pulse separation local maxima.

**Supplementary references**

S1. Schibli T R, Kim J, Kuzucu O, et al. Attosecond active synchronization of passively mode-locked lasers by balanced cross correlation. Optics Letters, 2003, 28, 947-949.

S2. Kim J, Kaertner F X. Attosecond-precision ultrafast photonics. Laser & Photonics Reviews, 2010, 4(3): 432-456.

S3. Song Y, Zhou F, Tian H, et al. Attosecond timing jitter within a temporal soliton molecule. Optica, 2020, 7(11): 1531-1534.
S4. S1-Nimmesgern, L., Beckh, C., Kempf, H., Leitenstorfer, A., & Herink, G. (2021). Soliton molecules in femtosecond fiber lasers: universal binding mechanism and direct electronic control. Optica, 8(10), 1334-1339.

S5. F. Takens, Detecting strange attractors in turbulence, Lecture Notes in Mathematics, Vol. 898 (1981) p. 366.

S6. Fraser, A. M., & Swinney, H. L. (1986). Independent coordinates for strange attractors from mutual information. Physical Review A, 33(2), 1134.

S7. Kennel, M. B., Brown, R., & Abarbanel, H. D. (1992). Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical Review A, 45(6), 3403.

S8. S2-Rosenstein, M. T., Collins, J. J., & De Luca, C. J. (1993). A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 65(1-2), 117-134.

S9. Grassberger, P., & Procaccia, I. (1983). Characterization of strange attractors. Physical Review Letters, 50(5), 346.