Formation and evolution dynamics of soliton molecules in spatiotemporal mode-locked fiber laser

Abstract

Spatiotemporal mode-locked (STML) fiber lasers have garnered significant interest in recent years for their unique ability to generate high-energy ultrashort pulses while simultaneously supporting complex transverse field structures. In particular, the formation of soliton molecules opens new avenues for ultrafast signal processing, high-bit-rate optical communications, and precision metrology. In this work, we presented a mode-resolved investigation into the birth, temporal evolution, and nonlinear dynamics of multiple soliton molecules in an all-fiber STML laser. Additionally, our mode-resolved analysis revealed that the transverse modes evolved synchronously during soliton-molecule formation while maintaining a relatively stable modal energy distribution. We further observed that both pulse spacing and relative phase gradually converged to constant values after a transient stage, indicating the formation of stable soliton molecules sustained by the balance among dispersion, nonlinearity, and gain-loss dynamics. Our study not only advances fundamental understanding of multi-soliton physics in complex laser cavities but also provides a numerical basis for engineering tailored spatiotemporal soliton molecules with controllable molecular order, spacing, and stability.

Graphical abstract

Keywords

fiber laser / spatiotemporal mode-locked / soliton molecules / nonlinear dynamics

1 Introduction

Mode-locked fiber lasers serve as a fundamental platform for generating ultrashort optical pulses and exploring nonlinear dynamics. Owing to their compact architecture, high environmental stability, and compatibility with all-fiber integration, they have been widely used in ultrafast spectroscopy, optical frequency combs, nonlinear optics, precision metrology, and high-speed optical communications [1–4]. In conventional single-mode fiber lasers, the pulse evolution is mainly governed by the balance among group-velocity dispersion, Kerr nonlinearity, and saturable absorption. Depending on the intracavity parameters, these systems can support diverse mode-locking regimes, including conventional solitons, self-similar pulses, dissipative solitons, noise-like pulses, and multi-soliton states [5, 6]. Among them, dissipative solitons and their bound states provide an important route for understanding the self-organization of ultrashort pulses in nonlinear optical cavities.

As the pump power or intracavity gain increases, a single soliton may become unstable and split into multiple pulses due to nonlinear phase accumulation [7]. Under suitable dissipative conditions, these pulses can form stable bound states with fixed temporal separations and locked relative phases, known as soliton molecules. Similar to chemical molecules composed of interacting atoms, optical soliton molecules consist of individual solitons coupled through dispersion, gain depletion and recovery, cross-phase modulation, and long-range nonlinear interactions [8]. In conventional single-mode fiber lasers, soliton molecules have been extensively investigated in terms of their formation, vibration, internal motion, phase evolution, and controllable switching dynamics [9, 10]. These studies have established soliton molecules as a representative nonlinear dynamical state and a promising carrier for multi-level optical encoding and ultrafast information processing [11–13].

Unlike conventional single-mode systems, spatiotemporal mode-locking (STML) fiber lasers support the generation of high-energy ultrashort pulses with rich transverse modal dynamics, offering unprecedented opportunities in fields such as high-capacity optical communication, nonlinear optics, and ultrafast spectroscopy [14–25]. While soliton molecules have been widely studied in single-mode systems, their behavior in multimode and spatiotemporal configurations remains largely unexplored. Multimode fibers, particularly graded-index multimode fibers (GIMFs), provide a natural platform for supporting multiple transverse modes and complex nonlinear interactions, which can dramatically enrich the dynamics of soliton molecules [26–30]. However, the interplay between modal dispersion, spatial filtering, saturable absorption, and nonlinear coupling presents significant challenges for both experimental observation and theoretical modeling.

In 2017, Wise et al. [31] first achieved spatiotemporal mode-locked output in a multimode fiber laser, opening a new era in the field of fiber lasers. In 2018, Qin et al. [32] first observed soliton molecules in a spatiotemporal mode-locked fiber laser. Different soliton molecules, including soliton pairs, soliton triplets, and soliton quartets, were obtained by adjusting the polarization state in the cavity, revealing the molecular formation mechanism and tunability in the context of spatiotemporal coupling. In 2021, Guo et al. [33] developed a single-shot multispeckle spectral-temporal measurement technique to reveal the real-time speckle-resolved spectral-temporal dynamics of spatiotemporal solitons in a 3D soliton laser. In 2023, Guo et al. [34] used a multi-speckle spatiotemporal measurement technique to achieve real-time, speckle-resolved observation of three-dimensional soliton molecules over long timescales, thereby capturing for the first time a variety of their dynamical behaviors. In 2025, Ning et al. [35] inserted a rotatable silica window into an STML fiber laser to introduce an extra spatial degree of freedom without changing key cavity parameters, thereby enabling rotation-controlled switching among multiple dynamic regimes via mode-energy redistribution and spatial filtering, with the switching mechanism further supported by misalignment simulations and field tests. Recently, Zhang et al. [36] investigated the transient spectral evolution of multi-soliton states in an all-fiber STML laser and observed the multimode double-pulse state in real time using time-stretch dispersive Fourier transform. Despite these advances, previous studies have mainly focused on the experimental observation and real-time characterization of spatiotemporal soliton molecules, whereas the mode-resolved physical mechanisms governing their birth, binding, and stabilization remain insufficiently clarified. In particular, how modal dispersion, spatial filtering, gain saturation, saturable absorption, and intermodal nonlinear coupling jointly drive the transition from a single spatiotemporal soliton to stable molecular states has not yet been systematically investigated. Therefore, a numerical framework capable of resolving the evolution of individual transverse modes is essential for revealing the physical origin of soliton-molecule formation in STML fiber lasers.

In this paper, we present a comprehensive numerical study of the formation, evolution, and stabilization of soliton molecules in an all-fiber STML laser cavity. Based on solving the generalized multimode nonlinear Schrödinger equation (GMMNLSE) and incorporating gain, saturable absorption, spatial filtering, and multimode coupling, our model enables mode-resolved tracking of the spatiotemporal pulse evolution, which is difficult to access directly in experiments. Different from previous studies that mainly focused on the experimental observation of spatiotemporal soliton molecules, this work provides a mode-resolved numerical interpretation of how modal dispersion, weak intermodal coupling, gain saturation, and intracavity dissipative effects jointly drive the transition from a single spatiotemporal soliton to stable molecular states. We systematically investigate the gain-dependent evolution from single-pulse operation to soliton pairs and soliton triplets, and reveal the stabilization dynamics of pulse separation and relative phase during molecular stabilization. Furthermore, the emergence of soliton triplets is identified as a higher-order spatiotemporal bound state enabled by enhanced intracavity energy and collective multimode synchronization, while the supported transverse modes maintain a relatively stable modal energy distribution. These results provide physical insight into the mode-resolved bonding mechanism of spatiotemporal soliton molecules and offer a numerical basis for tailoring their molecular order, spacing, and stability in multimode fiber lasers.

2 Simulation methods

The laser setup employed in the simulation is illustrated in Fig. 1. Each cavity element was represented by an appropriate mathematical model in the simulation to capture intracavity pulse dynamics. The pulses were first injected into a graded-index multimode fiber (GIMF), then sequentially traveled across a gain fiber, an output coupler, and a saturable absorber (SA). These components were connected by 0.2 m sections of GIMF, whose core diameter was 50 μm.

The numerical simulation model was based on the GMMNLSE, which can be written by Eq. (1):

(1)\begin{array}{rl}{\mathrm{\partial }}_z{A}_p(z,t)=& \mathrm{i}({\beta }_0^{(p)}-{\beta }_0^{(1)})A_p-({\beta }_1^{(p)}-{\beta }_1^{(1)}){\mathrm{\partial }}_t{A}_p\\ & +\sum _{j=2}^{N_d}\frac{i^{j+1}}{j!}{\beta }_j^{(p)}{\mathrm{\partial }}_t^j{A}_p\\ & +\mathrm{i}\frac{n_2{\omega }_0}{c}\sum _{l,m,n=1}^N{S}_{plmn}^K{A}_l{A}_m{A}_n^{\ast },\end{array}

where A_p denotes the optical field envelope of mode p, {\beta }_j^{(p)} represents the j-th order dispersion coefficient of mode p, n_2 is the nonlinear refractive index, which is set to 2.3\times {10}^{-20}{\mathrm{m}}^2/\mathrm{W} in the simulations, {\omega }_0 signifies the central angular frequency, S_{plmn}^K stands for the nonlinear mode-coupling tensors between modes p, l, m, and n for Kerr nonlinearity. In the present simulations, Raman scattering and self-steepening are neglected to reduce computational complexity. These higher-order effects are expected to play a secondary role for the pulse durations and energy levels considered here and thus do not qualitatively affect the main conclusions. The modal propagation and dispersion parameters used in the simulations are summarized in Table 1. Here, \mathrm{\Delta }{\beta }_0={\beta }_0^{(p)}-{\beta }_0^{(1)} and \mathrm{\Delta }{\beta }_1={\beta }_1^{(p)}-{\beta }_1^{(1)} are given relative to the fundamental {\mathrm{L}\mathrm{P}}_{01} mode, while {\beta }_2 – {\beta }_4 denote the higher-order dispersion coefficients of each mode.

In addition, the gain section is modeled using a simplified multimode gain model, in which six transverse modes are retained and identical gain parameters are applied to each mode. The propagation of each modal envelope A_p (p=1,\dots ,6) through the gain section is described by Eq. (2):

(2)\begin{array}{rl}\frac{\mathrm{\partial }{A}_p(z,\omega )}{\mathrm{\partial }z}=& \mathrm{i}\gamma {\left|A_p(z,\omega )\right|}^2{A}_p(z,\omega )\\ & +\frac{g_{0}f(\omega )}{2\left(1+\int {\left|A_p(z,\omega )\right|}^2\mathrm{d}\omega /E_{\mathrm{s}\mathrm{a}\mathrm{t}}\right)}{A}_p(z,\omega ),\end{array}

where \gamma =2.8{\mathrm{W}}^{-1}\cdot {\mathrm{k}\mathrm{m}}^{-1} is the nonlinear coefficient, E_{\mathrm{s}\mathrm{a}\mathrm{t}}=0.9\mathrm{n}\mathrm{J} is the gain saturation energy, g_0 is the small-signal gain, and f(\omega ) corresponds to a Gaussian-shaped gain spectrum with a full width at half maximum (FWHM) of 40 nm and a central wavelength of 1560 nm. The individual temporal and spectral profiles of all six transverse modes are retained during propagation. The gain section has a length of 0.4\mathrm{m}, and different soliton molecule states can be obtained by adjusting g_0. In the present model, the second-order dispersion of the gain fiber is neglected as a simplifying approximation, allowing the analysis to focus on gain saturation, nonlinear effects, and multimode dynamics associated with soliton molecule formation. This approximation may introduce moderate quantitative differences in the pulse width, pulse separation, spectral bandwidth, and spectral modulation spacing, while the principal formation dynamics and qualitative conclusions are not significantly altered within the parameter range considered here. Additional comparative simulations including the gain-fiber dispersion are presented in the Discussion section to further evaluate its influence on the soliton-molecule dynamics.

The linear mode coupling between the GIMF and the gain fiber is described by an effective 6\times 6 mode coupling matrix M. In principle, the coupling coefficients between transverse modes of two adjacent fiber sections can be evaluated from their normalized spatial overlap integrals, as given by Eq. (3):

(3)M_{ij}={\displaystyle \frac{{\displaystyle \iint F_i^{(2)\ast }(x,y)F_j^{(1)}(x,y)\mathrm{d}x\mathrm{d}y}}{\sqrt{{\displaystyle \iint {\left|F_i^{(2)}(x,y)\right|}^2\mathrm{d}x\mathrm{d}y}}\sqrt{{\displaystyle \iint {\left|F_j^{(1)}(x,y)\right|}^2\mathrm{d}x\mathrm{d}y}}}},

where F_j^{(1)} and F_i^{(2)} refer to the transverse mode profiles of the two adjacent fiber sections [37]. The overlap integral formulation provides the physical basis for constructing the mode coupling matrix. In the present simulation, the specific elements of \mathit{M} are empirically selected through numerical parameter adjustment to represent the effective mode coupling, mode-dependent loss, and spatial filtering at the fiber interface, with the diagonal and off-diagonal elements describing modal retention and intermodal coupling, respectively [31, 38]. The mode fields are transformed according to A_{\mathrm{o}\mathrm{u}\mathrm{t}}=\mathit{M}{A}_{\mathrm{i}\mathrm{n}}, where A_{\mathrm{o}\mathrm{u}\mathrm{t}} and A_{\mathrm{i}\mathrm{n}} denote the output and input mode fields, respectively. Since \mathit{M} represents an effective lossy coupling process rather than a lossless transformation between complete modal bases, it is not required to satisfy the unitary condition {\mathit{M}}^{†}\mathit{M}=\mathit{I} [37]. The maximum singular value of the matrix is approximately 0.478, confirming that the coupling is passive and does not introduce additional optical gain. The corresponding power reduction represents effective coupling and spatial filtering loss outside the retained six-mode subspace [18]. In this paper, \mathit{M} is set as

(4)\mathit{M}=\left(\begin{array}{cccccc}0.3& 0.1& 0& 0& 0& 0.1\\ 0& 0.25& 0.15& 0& 0& 0\\ 0.1& 0.2& 0.2& 0.2& 0& 0\\ 0& 0& 0.1& 0.1& 0& 0\\ 0& 0& 0& 0.1& 0.15& 0.1\\ 0& 0& 0& 0& 0.1& 0.2\end{array}\right).

For SA, the spatiotemporal nonlinear transmittance T is described by Eq. (5):

(5)T(z,t)=1-{\displaystyle \frac{{\alpha }_s}{1+|A(z,t){|}^2/I_{\mathrm{s}}}}-{\alpha }_{\mathrm{n}\mathrm{s}},

where the non-saturable absorption coefficient {\alpha }_{ns} = 0.5, {\alpha }_s = 0.2 is the modulation depth, governing the pulse-shaping dynamics in mode-locked lasers, I_s = 10 W is the saturation intensity, which determines the threshold for intensity-dependent absorption saturation. The modulation depth {\alpha }_s = 0.2 is chosen based on experimentally reported saturable absorbers operating around the 1.5 μm wavelength region, where modulation depths of approximately 20% have been demonstrated [39]. The non-saturable absorption coefficient {\alpha }_{ns} = 0.5 is treated as an effective loss parameter, including the intrinsic non-saturable absorption, insertion loss, and spatial filtering loss in the simulated STML cavity.

In the present simulation, six transverse modes of the GIMF are retained, including {\mathrm{L}\mathrm{P}}_{01}, {\mathrm{L}\mathrm{P}}_{11\mathrm{a}}, {\mathrm{L}\mathrm{P}}_{11\mathrm{b}}, {\mathrm{L}\mathrm{P}}_{21\mathrm{a}}, {\mathrm{L}\mathrm{P}}_{21\mathrm{b}}, and {\mathrm{L}\mathrm{P}}_{02}. In the following, we denote these modes as mode 1−mode 6 corresponding to {\mathrm{L}\mathrm{P}}_{01}, {\mathrm{L}\mathrm{P}}_{11\mathrm{a}}, {\mathrm{L}\mathrm{P}}_{11\mathrm{b}}, {\mathrm{L}\mathrm{P}}_{21\mathrm{a}}, {\mathrm{L}\mathrm{P}}_{21\mathrm{b}}, and {\mathrm{L}\mathrm{P}}_{02}, respectively. The central wavelength is set to 1550 nm. Each mode is initialized with a Gaussian-shaped pulse with a FWHM of 10 ps. The initial total pulse energy is 50 pJ. The Gaussian pulse is used as the initial condition to facilitate the analysis of the subsequent pulse splitting and molecule formation dynamics. While the final molecule state is primarily determined by the cavity parameters, the detailed transient evolution path can depend on the initial field, particularly during the early buildup stage. Accordingly, quantitative features such as the pulse width, pulse separation, and temporal drift may show slight variations for different initial conditions. Additional simulations using a noise floor as the initial condition are presented in the Discussion section to further evaluate the influence of the initial field on the formation and stabilization of the soliton molecules.

3 Results and discussion

3.1 Single pulse

We first set g_0 to 11.6 {\mathrm{m}}^{-1}, under which the laser cavity supported the formation of a stable single spatiotemporal soliton. Figure 2 illustrated the temporal profile and spectrum of the stabilized pulse after 500 cavity roundtrips. Figure 2(a) presented the temporal profile, which exhibited a clean and symmetric shape without side lobes, indicating the absence of multi-pulse splitting. As illustrated in the inset, the beam profile confirmed that the pulse supported multiple modes. The pulse width of each mode was on the order of hundreds of femtoseconds. The corresponding spectra, shown in Fig. 2(b), displayed a smooth {\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}}^2-like envelope, with the central wavelengths of all modes located around 1552 nm. These temporal and spectral characteristics further confirm the stable single-pulse operation. The stability of the pulse was further examined by monitoring its evolution over multiple cavity roundtrips. The intracavity energy evolution, presented in Fig. 2(c), stabilized after around 150 roundtrips without exhibiting oscillatory behavior, which highlighted the robustness of the single soliton state under moderate gain conditions.

It should be noted that the simulation center wavelength is set to 1550 nm, whereas the gain spectrum has a broad bandwidth of 40 nm. Therefore, although the gain center is located at 1560 nm, the spectral region around 1550-1552 nm still experiences sufficient gain, and the stabilized pulse spectrum is not necessarily centered at the maximum of the gain profile. To examine this effect, we compared different wavelength settings, as shown in Fig. 3. When both the simulation center wavelength and the gain center are set to 1550 nm, the system remains in the same stable single-pulse state and the spectrum is still centered at approximately 1552 nm, as shown in Figs. 3(a, c). In contrast, Figs. 3(b, d) show that when both wavelengths are shifted to 1560 nm, the cavity evolves into a stable soliton pair. These results indicate that the apparent offset between the gain center and the stabilized spectrum mainly reflects the different wavelength settings of the gain profile and the simulated intracavity field, while the final state is determined by the overall balance among dispersion, nonlinearity, gain, saturable absorption, modal coupling, and cavity loss.

3.2 Soliton pairs

When g_0 was increased to 12.6 {\mathrm{m}}^{-1}, stable soliton pairs were obtained. This was due to the higher gain leading to excessive nonlinear phase accumulation in the solitons, causing a single soliton to split into multiple solitons, which formed soliton molecules under the combined effects of intracavity nonlinearity and dispersion. The characteristics of soliton pairs were illustrated in Fig. 4. Figures 4(a) and (b) presented the stabilized pulse profiles and the spectrum of soliton molecules after 500 resonator roundtrips. The pulse separation \mathrm{\Delta }\tau was 1.5 ps and the modulation period \mathrm{\Delta }\lambda was 4.9 nm. It can be verified that \mathrm{\Delta }\tau \cdot c\cdot \mathrm{\Delta }\lambda /{\lambda }^2=1 approximately where the central wavelength \lambda was 1552 nm and c is the speed of light, demonstrating the output corresponds to a typical soliton molecule.

In order to study the formation and evolution of soliton pairs, the evolution dynamics of the soliton pair were analyzed. As the intensities of different modes varied significantly, each mode was independently normalized to more clearly reveal its evolution process. Figures 5(a) and (b) illustrated the temporal and spectral evolution of six modes in the spatiotemporal mode-locked fiber laser. The temporal and spectral evolution of the different modes exhibited similar trends, suggesting collective multimode evolution under the combined effects of dispersion, nonlinearity, modal coupling, and gain-loss dynamics. Such collective dynamics enabled the establishment of stable phase relationships and fixed pulse separations, which were essential for the formation of soliton molecules [41].

Although the intensity levels of different transverse modes are different, all six modes exhibit nearly synchronized temporal splitting and spectral modulation during the formation of the soliton pair. This behavior indicates that the constituent pulses are not independently formed in separate modes; instead, they evolve as a collectively locked spatiotemporal wave packet. The physical origin of the synchronized multimode evolution can be understood from the joint action of modal walk-off, intermodal coupling, nonlinear interaction, and cavity dissipation. The relatively small group-delay mismatch among the supported transverse modes preserves their temporal overlap during propagation, providing favorable conditions for sustained intermodal interactions. Under this weak walk-off condition, Kerr nonlinear interactions and mode coupling correlate the temporal and spectral evolution of the different modes, while spatial filtering imposes a common round-trip constraint on the multimode field. In addition, gain saturation and saturable absorption regulate the intracavity pulse energy and stabilize the pulse evolution, thereby providing the dissipative conditions required for the formation of a stable collectively locked state. Consequently, the supported transverse modes evolve collectively toward a stable spatiotemporal bound state rather than behaving as independent pulse trains. The equilibrium pulse separation is manifested as a common temporal separation shared by all supported modes rather than as mode-dependent separations. This differs from conventional single-mode soliton molecules, whose bonding characteristics are determined solely by scalar temporal effects such as dispersion, Kerr nonlinearity, gain saturation, and saturable absorption. In the present STML cavity, the bonding process additionally depends on the preservation of modal overlap and collective transverse-mode evolution under weak intermodal walk-off.

Figure 6 illustrated the energy evolution of the spatiotemporal mode-locked laser. Figure 6(a) was the energy evolution of the entire process. Low-order modes, especially mode 3, dominated the cavity. In the initial stage, the intracavity energy increased rapidly owing to amplification in the gain fiber. After a period of time, the growth rate slowed down significantly. This was because the gain fiber reached saturation, which limited the further increase of the energy in the cavity. After about 50 roundtrips, the energy in the cavity began to oscillate, indicating that the cavity could not maintain a stable single pulse state. After the 150th roundtrip, the energy in the cavity stabilized, indicating the formation of a stable soliton pair. Figure 6(b) presented the energy evolution of the formation process of soliton pairs. During the formation of soliton pairs, the changing trends of the energies of various modes were basically the same. It should be noted that no pronounced modal energy redistribution was observed during the formation of the soliton pair. Instead, the energies of different modes followed similar evolution trends and finally reached stable values. This indicates that in this GIMF-based cavity, the formation of the soliton molecule is not driven by strong energy transfer among transverse modes. Rather, the multimode field maintains a relatively stable modal energy distribution while the temporal pulse structure evolves from a single-pulse state into a bound molecular state. Therefore, the intermodal interaction in this regime is mainly reflected in collective modal synchronization rather than strong modal energy exchange. Such mode-resolved energy information provides a complementary perspective to experimental observations, where the measured field generally contains the superposition of multiple transverse modes and the energy contribution of each mode is difficult to isolate.

Figure 7 illustrated the temporal evolution of the relative phase and pulse separation between the two constituent pulses of the soliton pairs. In the initial stage, the intracavity pulse built up from the initial condition, with its energy gradually accumulating and the waveform being temporally compressed under the combined action of the gain, saturable absorption, Kerr nonlinearity, and spatial filtering. Around the 100th cavity roundtrip, the initially formed STML single pulse became unstable and underwent a splitting event, initiating the formation of a two-pulse bound state. In the spatiotemporal mode-locking framework, this transition was driven by the coupled evolution of the temporal waveform and the transverse modal field: intermodal coupling, together with the effective spatial filtering, modified the dissipative balance of the cavity and favored a multi-pulse attractor over the single-pulse one. Subsequently, the two pulses separated in time while their relative phase evolved under the combined action of Kerr nonlinearity, dispersion, and gain-loss modulation. After approximately 200 cavity roundtrips, the phase difference and pulse spacing converged to constant values. The relative phase remained locked, and the pulse separation stabilized around 1.5 ps, consistent with the modulation period observed in the optical spectrum. The stabilization of pulse separation and relative phase further confirms the establishment of a stable soliton molecule state. In conventional single-mode fiber lasers, such stabilization is mainly governed by scalar temporal interactions, including dispersion, Kerr nonlinearity, gain saturation, and saturable absorption. In contrast, in the STML cavity considered here, the temporal binding process occurs in the presence of multiple transverse modes. Owing to the relatively small modal group-delay mismatch in the GIMF, the supported modes remain temporally overlapped during propagation and participate in molecular formation in a synchronized manner. As a result, the final soliton molecule can be regarded as a weak walk-off, collectively locked spatiotemporal bound state rather than a purely temporal bound state, as observed in conventional single-mode fiber lasers.

3.3 Soliton triplets

We also obtained soliton triplets by adjusting g_0 to 12.9 {\mathrm{m}}^{-1}. Figure 8 presented the temporal and spectral characteristics, as well as the evolution dynamics, of soliton triplets. As shown in Fig. 8(a), the temporal profile exhibited three well-separated pulses with nearly identical peak intensities, forming a well-defined three-pulse temporal structure. The corresponding optical spectrum in Fig. 8(b) displayed pronounced modulation, confirming the multi-pulse interference typical of higher-order soliton molecules. Figure 8(c) showed the intracavity energy evolution, which stabilized after approximately 150 roundtrips, confirming the formation of a robust soliton triplet driven by strong nonlinearity and mode interactions. It should be noted that the energy shown in Fig. 8(c) represents a stable but non-optimized operating condition. The parameters used here were chosen to clearly track the formation and stabilization of the soliton triplet, rather than to maximize the output energy. Figures 8(d) and (e) illustrated the temporal and spectral evolution of the dominant transverse mode (mode 3), respectively, during the formation of the soliton triplet. A stable three-pulse structure emerged after a transient growth and splitting phase.

Compared with the soliton pair state, the formation of soliton triplets requires a higher intracavity gain, which provides sufficient energy for further pulse splitting while the cavity maintains a stable dissipative balance. Although the six transverse modes have different propagation constants and group velocities, the temporal positions of the three pulses remain consistent among the different modes, indicating that the triplet is not independently formed in a single dominant mode. Instead, the constituent pulses are established through synchronized multimode evolution. The stable modal energy evolution in Fig. 8(c) further indicates that no pronounced modal energy redistribution occurs during the stabilization stage; rather, the multimode field reaches a stable modal energy distribution as the temporal three-pulse bound state forms. The temporal and spectral evolution of mode 3 in Figs. 8(d) and (e) clearly shows the transition from an initially broadened pulse to a stable three-pulse molecular state, accompanied by the emergence of regular spectral modulation. This behavior confirms that the higher-order soliton molecule is sustained by the joint balance among gain saturation, Kerr nonlinearity, saturable absorption, spatial filtering, and weak intermodal coupling. Therefore, the soliton triplet observed here should be regarded as a higher-order spatiotemporal bound state with collective modal locking rather than a simple extension of conventional single-mode soliton triplets.

3.4 Discussion

Figure 9 depicted the relationship between the number of soliton pulses within a soliton molecule and the small-signal gain g_0. When g_0 was below the mode-locking threshold, the intracavity field could not build up to form a stable pulse. Conversely, excessively large g_0 drove the system into an unstable regime due to the breakdown of the dispersion-nonlinearity balance. As the gain increased, the accumulated intracavity energy surpassed the threshold for multi-pulse formation, leading to the emergence of higher-order soliton molecules. The graph demonstrated a step-like increase in the number of pulses with discrete gain intervals, suggesting that the number of pulses was quantized and dictated by the balance between gain saturation, nonlinear interaction, and spatial filtering. This result indicated that by adjusting g_0, one can effectively control the order and structure of the soliton molecule, offering a potential mechanism for tunable ultrafast pulse generation in multimode laser systems.

Although the quantitative characteristics of the obtained soliton molecules depend on the specific coupling matrix and fiber parameters, the underlying collective multimode mechanism is expected to persist over a finite parameter range. Variations in the coupling matrix alter the coupling and filtering of different transverse modes in the cavity, which can modify their relative energy distribution and may produce a different dominant mode in the final stable state. However, collective multimode evolution can still occur provided that several transverse modes remain active, are sufficiently coupled, and are not completely suppressed by spatial filtering. The modal group-delay mismatch is also a critical parameter. When the accumulated intermodal walk-off \mathrm{\Delta }{T}_{pq}=|\mathrm{\Delta }{\beta }_{1,pq}|L, where L denotes the length of the GIMF, remains comparable to or smaller than the characteristic pulse duration, the transverse modes retain sufficient temporal overlap for sustained interactions. Moderate variations in {\beta }_2, the nonlinear coefficients, and the modal overlap coefficients are expected primarily to modify quantitative features, including the pulse duration, spectral bandwidth, molecular separation, and gain thresholds for the formation of soliton pairs and triplets, rather than altering the fundamental collective synchronization mechanism.

In order to further evaluate the influence of gain fiber dispersion, we performed additional simulations by including a typical second-order dispersion coefficient of {\beta }_2=16.5\mathrm{p}{\mathrm{s}}^2/\mathrm{k}\mathrm{m} in the 0.4 m gain fiber, while keeping the other cavity parameters unchanged. For g_0=12.6{\mathrm{m}}^{-1}, the system still evolves into a stable soliton pair, as shown in Fig. 10. During the initial stage, the intracavity energy increases rapidly and the pulse undergoes temporal reshaping and splitting. The two constituent pulses then gradually establish a stable temporal separation, while the corresponding spectral interference fringes become increasingly regular. The six transverse modes exhibit nearly synchronized temporal and spectral evolution and finally converge to a stable collectively locked state. Compared with the case without gain-fiber dispersion, the final pulse separation changes to approximately 1.1 ps and the spectral modulation period to approximately 6.5 nm, together with moderate changes in the pulse duration and spectral bandwidth.

For g_0=12.9{\mathrm{m}}^{-1}, the system similarly evolves through an initial energy buildup and pulse splitting stage before forming a stable soliton triplet, as shown in Fig. 11. Three localized pulses gradually emerge and establish a stable bound state, while the spectral modulation structure becomes stable after the transient evolution. The synchronized evolution among the six transverse modes is also preserved. The final adjacent pulse separations are approximately 1.15 ps and 1.1 ps, with additional changes in the pulse duration, spectral bandwidth, and modulation period compared with the case without gain-fiber dispersion.

These quantitative differences arise because the additional dispersion in the gain fiber modifies the pulse chirp and temporal width during each cavity roundtrip, which consequently changes the peak power and nonlinear phase accumulation. Through the repeated interplay of dispersion, Kerr nonlinearity, gain saturation, and saturable absorption, the equilibrium temporal and spectral parameters of the soliton molecules are therefore shifted. Nevertheless, under the operating conditions considered here, the gain fiber dispersion does not change the soliton molecule order, the principal formation pathway, or the synchronized multimode evolution. Therefore, its omission in the original model mainly affects quantitative pulse characteristics rather than the qualitative conclusions of this work.

In addition, we performed additional simulations using a noise floor while keeping all other cavity parameters unchanged to further examine the influence of the initial condition. For g_0=12.6{\mathrm{m}}^{-1}, the intracavity field first builds up from random fluctuations, followed by pulse localization, splitting, and binding, and finally evolves into a stable soliton pair, as shown in Fig. 12. The final pulse separation is approximately 1.63 ps and the spectral modulation period is approximately 4.5 nm. Compared with the Gaussian initial condition, the overall molecular formation process and final two-pulse state are preserved, although moderate differences appear in the pulse width and equilibrium separation.

For g_0=12.9{\mathrm{m}}^{-1}, the noise-seeded field similarly undergoes an initial buildup and self-organization stage before forming a stable soliton triplet, as shown in Fig. 13. The final adjacent pulse separations are approximately 1.74 ps and 2.36 ps. Compared with the Gaussian initial condition, the three-pulse molecular order is preserved, while the pulse separations, pulse widths, and detailed temporal configuration show more noticeable variations. This stronger sensitivity can be attributed to the more complex interactions among three constituent pulses, which involve multiple pulse separations and relative phases.

These comparisons indicate that the initial condition mainly affects the transient formation pathway and the quantitative characteristics of the final molecular state. A Gaussian initial pulse already provides a localized field, whereas a noise floor requires an additional buildup and self-organization process governed by gain, saturable absorption, dispersion, Kerr nonlinearity, and spatial filtering. Consequently, different initial fields can lead to different transient gain saturation, nonlinear phase accumulation, and relative phases during pulse formation. Nevertheless, under the operating conditions considered here, both initial conditions produce the same molecular order, confirming that the stable soliton pair and triplet states are not imposed by the pre-shaped Gaussian input.

Compared with conventional single-mode soliton molecules, the soliton molecules in STML fiber lasers possess additional transverse degrees of freedom. In single-mode systems, the molecular binding is described by the evolution of a scalar pulse field, and the key observables are mainly pulse separation, relative phase, and spectral modulation. In the present multimode system, these temporal molecular features are accompanied by stable modal participation and simultaneous evolution of multiple transverse modes. Therefore, the soliton molecule obtained here is not simply a temporal bound state, but a spatiotemporal bound state sustained by the joint balance of temporal dissipative dynamics and transverse mode-locking.

Our results are consistent with previous experimental observations of spatiotemporal soliton molecules, such as the real-time measurements reported by Guo et al. [34]. Guo et al. [34] experimentally investigated the real-time spatiotemporal dynamics of three-dimensional soliton molecules over long timescales, whereas the present numerical framework provides direct access to the evolution of individual transverse modes. By independently tracking the temporal profiles, spectra, and energies of the supported modes, our simulations reveal how the transverse modal components collectively evolve during molecular formation and eventually establish a stable spatiotemporal bound state. This mode-resolved perspective therefore complements real-time experimental measurements by providing information on the modal-level dynamics that is difficult to isolate directly in experiments.

It should be noted that the gain section employed in the present simulation is a simplified multimode model in which identical gain parameters are applied to all six transverse modes. This treatment preserves the mode-resolved temporal and spectral evolution throughout the cavity, but does not fully account for the mode-dependent properties of a realistic multimode active fiber. In practice, different transverse modes may experience different effective gain coefficients and saturation characteristics because of their different spatial overlaps with the doped region [42, 43]. Mode-dependent dispersion and nonlinear coefficients may also introduce additional temporal walk-off, nonlinear phase accumulation, and spectral differences among the modes [44]. In addition, gain competition may modify the relative modal amplification and steady-state energy distribution [45]. These effects could influence the relative modal energies, dominant transverse mode, degree of synchronization, gain thresholds for different molecule states, equilibrium pulse separation, and corresponding stability ranges. Therefore, the specific operating points and quantitative values obtained here should be interpreted within the framework of the present simplified gain model. Nevertheless, the model retains the transverse modal degrees of freedom throughout the cavity, and the qualitative picture of collective multimode evolution is expected to remain applicable when multiple modes maintain sufficient temporal overlap and coupling.

4 Conclusion

In this work, we numerically investigated the formation and dynamic evolution of soliton molecules in a spatiotemporal mode-locked fiber laser based on graded-index multimode fibers. Using the GMMNLSE together with models for gain, saturable absorption, and spatial coupling, we successfully simulated the generation of stable soliton molecules including soliton pairs and soliton triplets under different small-signal gain conditions. Our work provides access to mode-resolved quantities, including individual modal energy evolution and modal synchronization, which are difficult to directly resolve in experiments and are essential for understanding the bonding mechanism of STML soliton molecules. Our results highlight the potential of spatiotemporal mode-locked lasers as a versatile platform for exploring complex nonlinear dynamics and developing next-generation ultrafast photonic systems.

References

[1] J. Wang , X. Wang , J. Lei , M. Ma , C. Wang , Y. Ge , and Z. Wei , Recent advances in mode-locked fiber lasers based on two-dimensional materials, Nanophotonics 9, 2315 (2020)

[2] B. Fu , J. Sun , C. Wang , C. Shang , L. Xu , J. Li , H. Zhang , MXenes: Synthesis , and optical properties , and applications in ultrafast photonics, Small 17, 2006054 (2021)

[3] Y. Han , Y. Guo , B. Gao , C. Ma , R. Zhang , H. Zhang , Generation , and optimization , and application of ultrashort femtosecond pulse in mode-locked fiber lasers, Prog. Quantum Electron. 71, 100264 (2020)

[4] C. Zhang , L. Cong , J. Zhang , H. Meng , S. Peng , S. Li , H. Zhang , Y. Wang , B. Fu , J. Wang , and M. Qiu , Advanced metasurfaces: precise modulations of light in space and time domains, Adv. Photonics 8, 014003 (2026)

[5] C. K. Nielsen and S. R. Keiding , All-fiber mode-locked fiber laser, Opt. Lett. 32, 1474 (2007)

[6] U. Keller , Recent developments in compact ultrafast lasers, Nature 424, 831 (2003)

[7] J. Peng and H. Zeng , Build-up of dissipative optical soliton molecules via diverse soliton interactions, Laser Photon. Rev. 12, 1800009 (2018)

[8] Z. Q. Wang , K. Nithyanandan , A. Coillet , P. Tchofo-Dinda , and P. Grelu , Optical soliton molecular complexes in a passively mode-locked fibre laser, Nat. Commun. 10, 830 (2019)

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

[10] C. Zhang , Z. Ouyang , Y. Lyu , B. Fu , L. Zhang , and M. Qiu , Modulation of ultrafast soliton molecules enabled by plasmonic metafibers, PhotoniX 6, 43 (2025)

[11] Y. Liu , S. Huang , Z. Li , H. Liu , Y. Sun , R. Xia , L. Yan , Y. Luo , H. Liu , G. Xu , Q. Sun , X. Tang , and P. P. Shum , Phase-tailored assembly and encoding of dissipative soliton molecules, Light Sci. Appl. 12, 123 (2023)

[12] S. Liu , Y. Cui , E. Karimi , and B. A. Malomed , On-demand harnessing of photonic soliton molecules, Optica 9, 240 (2022)

[13] K. Krupa , K. Nithyanandan , U. Andral , P. Tchofo-Dinda , and P. Grelu , Real-time observation of internal motion within ultrafast dissipative optical soliton molecules, Phys. Rev. Lett. 118, 243901 (2017)

[14] B. Cao , C. Gao , K. Liu , X. Xiao , C. Yang , and C. Bao , Spatiotemporal mode-locking and dissipative solitons in multimode fiber lasers, Light Sci. Appl. 12, 260 (2023)

[15] J. Sun , Y. Wang , C. Zhang , L. Xu , and B. Fu , Spatiotemporal nonlinear dynamics in multimode fiber laser based on carbon nanotubes, Front. Phys. (Beijing) 19, 52201 (2024)

[16] L. G. Wright , D. N. Christodoulides , and F. W. Wise , Controllable spatiotemporal nonlinear effects in multimode fibres, Nat. Photonics 9, 306 (2015)

[17] W. Wang , J. Duan , S. Ding , B. Yang , X. Cui , X. Meng , M. Shi , H. Huang , S. Zhao , F. He , J. Wang , X. Xu , Y. Yao , Z. Wang , T. Wang , J. Zhang , and F. Grillot , Tunable amplitude and frequency modulated optical frequency combs in external-cavity-locked quantum dot lasers, Adv. Photonics 7, 066011 (2025)

[18] L. G. Wright , P. Sidorenko , H. Pourbeyram , Z. M. Ziegler , A. Isichenko , B. A. Malomed , C. R. Menyuk , D. N. Christodoulides , and F. W. Wise , Mechanisms of spatiotemporal mode-locking, Nat. Phys. 16, 565 (2020)

[19] Y. Ding , X. Xiao , K. Liu , S. Fan , X. Zhang , and C. Yang , Spatiotemporal mode-locking in lasers with large modal dispersion, Phys. Rev. Lett. 126, 093901 (2021)

[20] G. Fu , Y. Li , Q. Xiao , D. Li , M. Gong , Q. Liu , S. Wabnitz , and P. Yan , Spatiotemporal dissipation dynamics: a route to high-beam-quality and high-peak-power spatiotemporal mode-locked fiber lasers, Optica 11, 1503 (2024)

[21] J. Sun , G. Wang , J. Chao , X. Wang , H. Yang , and B. Fu , Buildup of multiple spatiotemporal nonlinear dynamics in an all-fiber multimode laser, Opt. Lett. 48, 6019 (2023)

[22] C. Zhang , L. Zhang , H. Zhang , B. Fu , J. Wang , and M. Qiu , Pulsed polarized vortex beam enabled by metafiber lasers, PhotoniX 5, 36 (2024)

[23] J. Sun , G. Wang , C. Zhang , W. Lyu , Y. Lyu , and B. Fu , Buildup and real-time control of spatiotemporal nonlinear dynamics, Ultrafast Sci. 5, 0105 (2025)

[24] G. Wang , B. Fu , J. Sun , C. Zhang , C. Zhang , J. Wu , C. Brès , and L. Xu , Bridging the mode-locked and Q-switched regimes in fiber lasers, Laser Photon. Rev. 19, 2402171 (2025)

[25] J. Zou , H. He , L. Song , L. Huang , J. Lan , L. Lan , T. Chen , C. Kong , Z. Sun , and Z. Luo , Sub-100-fs deep-red mode-locked fiber laser for multicolor two-photon microscopy, Adv. Photonics 7, 046009 (2025)

[26] X. Lin , Y. Gao , J. Long , J. Wu , W. Hong , H. Cui , Z. Luo , W. Xu , and A. Luo , All few-mode fiber spatiotemporal mode-locked figure-eight laser, J. Lightwave Technol. 39, 5611 (2021)

[27] J. Wu , G. Liu , Y. Gao , X. Lin , H. Cui , Z. Luo , W. Xu , M. E. Likhachev , S. S. Aleshkina , V. M. Mashinsky , M. V. Yashkov , and A. Luo , Switchable femtosecond and picosecond spatiotemporal mode-locked fiber laser based on NALM and multimode interference filtering effects, Opt. Laser Technol. 155, 108414 (2022)

[28] D. Mao , Z. Yuan , K. Dai , Y. Chen , H. Ma , Q. Ling , J. Zheng , Y. Zhang , D. Chen , Y. Cui , Z. Sun , and B. A. Malomed , Temporal and spatiotemporal soliton molecules in ultrafast fibre lasers, Nanophotonics 14, 677 (2025)

[29] C. Zhang , W. Lyu , L. Cong , Z. Gu , Z. Zhang , G. Wang , Y. Lyu , B. Yang , I. Ahmad , and B. Fu , Tunable soliton molecules enabled by a nanotubes-based mode-locked fiber laser, Front. Phys. (Beijing) 20, 052201 (2025)

[30] Y. Ding , X. Xiao , P. Wang , and C. Yang , Multiple-soliton in spatiotemporal mode-locked multimode fiber lasers, Opt. Express 27, 11435 (2019)

[31] L. G. Wright , D. N. Christodoulides , and F. W. Wise , Spatiotemporal mode-locking in multimode fiber lasers, Science 358, 94 (2017)

[32] H. Qin , X. Xiao , P. Wang , and C. Yang , Observation of soliton molecules in a spatiotemporal mode-locked multimode fiber laser, Opt. Lett. 43, 1982 (2018)

[33] Y. Guo , X. Wen , W. Lin , W. Wang , X. Wei , and Z. Yang , Real-time multispeckle spectral-temporal measurement unveils the complexity of spatiotemporal solitons, Nat. Commun. 12, 67 (2021)

[34] Y. Guo , W. Lin , W. Wang , R. Zhang , T. Liu , Y. Xu , X. Wei , and Z. Yang , Unveiling the complexity of spatiotemporal soliton molecules in real time, Nat. Commun. 14, 2029 (2023)

[35] Y. Ning , J. Li , J. He , S. Wang , M. Xiao , Z. Li , P. Wang , Y. Liu , and Z. Wang , Controllable switching of soliton dynamics in spatiotemporal mode-locked fiber laser, Opt. Laser Technol. 192, 113926 (2025)

[36] W. Zhang , X. Xiao , H. Zhang , J. Peng , K. Chen , M. Sun , F. Yan , and Y. Zheng , Transient dynamics of multi-soliton states in a spatiotemporal mode-locked fiber laser, Opt. Laser Technol. 196, 114675 (2026)

[37] J. Vuong , P. Ramantanis , Y. Frignac , M. Salsi , P. Genevaux , D. F. Bendimerad , and G. Charlet , Mode coupling at connectors in mode-division multiplexed transmission over few-mode fiber, Opt. Express 23, 1438 (2015)

[38] T. Matsui , K. Nakajima , and T. Sakamoto , Effective mode-field diameter for few-mode fibers for considering splice loss characteristics, Appl. Opt. 56, 7484 (2017)

[39] Y. Song , S. Chen , Q. Zhang , L. Li , L. Zhao , H. Zhang , and D. Tang , Vector soliton fiber laser passively mode locked by few layer black phosphorus-based optical saturable absorber, Opt. Express 24, 25933 (2016)

[40] L. Yun and W. Zhao , PbS quantum dots saturable absorber for dual-wavelength solitons generation, Nanomaterials 11, 2561 (2021)

[41] D. Zou , R. Liu , Y. Shi , A. Zhang , J. Li , G. J. Chen , H. Dang , Y. Song , M. Hu , and P. P. Shum , Resonantly driven nonlinear dynamics of soliton molecules in ultrafast fiber lasers, Adv. Photonics 7, 016005 (2025)

[42] N. Bai , E. Ip , T. Wang , and G. Li , Multimode fiber amplifier with tunable modal gain using a reconfigurable multimode pump, Opt. Express 19, 16601 (2011)

[43] Z. Zhang , C. Guo , L. Cui , Y. Zhang , C. Du , and X. Li , All-fiber few-mode erbium-doped fiber amplifier supporting six spatial modes, Chin. Opt. Lett. 17, 100604 (2019)

[44] M. Zitelli , F. Mangini , M. Ferraro , O. Sidelnikov , and S. Wabnitz , Conditions for walk-off soliton generation in a multimode fiber, Commun. Phys. 4, 182 (2021)

[45] M. Gong , Y. Yuan , C. Li , P. Yan , H. Zhang , and S. Liao , Numerical modeling of transverse mode competition in strongly pumped multimode fiber lasers and amplifiers, Opt. Express 15, 3236 (2007)

Rights & permissions

Higher Education Press