Compiled

# On-demand harnessing of photonic soliton molecules

Shilong Liu  State Key Laboratory of Modern Optical Instrumentation, College
of Optical Science and Engineering, Zhejiang University,Hangzhou, Zhejiang
310027, China  Department of Physics, University of Ottawa, 25 Templeton
Street, Ottawa, Ontario, K1N 6N5 Canada  Corresponding author:
<EMAIL_ADDRESS>Yudong Cui  State Key Laboratory of Modern Optical
Instrumentation, College of Optical Science and Engineering, Zhejiang
University,Hangzhou, Zhejiang 310027, China  Corresponding author:
<EMAIL_ADDRESS>Ebrahim Karimi  Department of Physics, University of Ottawa,
25 Templeton Street, Ottawa, Ontario, K1N 6N5 Canada  Boris A. Malomed
Department of Physical Electronics, Faculty of Engineering, and Center for
Light-Matter Interaction, Tel Aviv University, Tel Aviv 69978, Israel
Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica,
Chile

###### Abstract

Soliton molecules (SMs) are fundamentally important modes in nonlinear optical
systems. It is a challenge to experimentally produce SMs with a required
temporal separation in mode-locked fiber lasers. Here, we propose and realize
an experimental scenario for harnessing SM dynamics in a laser setup. In
particular, we tailor SMs in a mode-locked laser controlled by second-order
group-velocity dispersion and dispersion losses: the real part of dispersion
maintains the balance between the dispersion and nonlinearity, while the
dispersion loss determines the balance of gain and losses. The experimental
results demonstrate that the dispersion loss makes it possible to select
desired values of the temporal separation (TS) in bound pairs of SMs in the
system. Tunability of the SM’s central wavelength and the corresponding
hysteresis are addressed too. The demonstrated regime allows us to create
multiple SMs with preselected values of the TS and central wavelength, which
shows the potential of our setup for the design of optical data-processing
schemes.

###### doi:

http://dx.doi.org/10.1364/optica.XX.XXXXXX

††journal: optica

##  1 Introduction

Photonic “soliton molecules" (SMs) are bound states of two or more fundamental
solitons, whose existence and stability are determined by the separation and
phase shift between the constituents [1, 2, 3, 4]. For a stable SM, the
relative phase of the bound solitons is fixed in the course of the
transmission. Nevertheless, the variation of the phase makes it possible to
observe complex dynamical effects, such as vibrations of the “molecule" [5],
stepwise evolution [6], and phase sliding [7, 8]. Direct studies of these
effects have been recently made possible with the help of the single-shot
measurement technique, enabled by the time-stretch dispersive Fourier
transform (TSDFT) [9, 10, 11, 12, 13]. Besides the relative phase, the
temporal separation (TS) between the bound solitons in the SM is another basic
degree of freedom. It directly affects the force of the interaction between
the solitons, which is strong for TS taking values comparable to the single-
soliton’s width, and weak for essentially larger separations [2, 14, 15, 8,
16]. Much interest has been drawn to the TS degree of freedom due to the fact
that it affects numerous phenomena, including nonlinear evolution [7, 6, 17],
short- and long-range interactions [3, 16], SM complexes [8], and the Hopf
bifurcation which excites robust intrinsic oscillations in SMs [18, 5]. The TS
also plays the basic role in the development of various SM-based optical
applications, such as “multialphabetic" encoding for communications [19],
optical switching [20], and all-optical storage [21]. These studies make it
necessary to develop means for “on-demand" manipulations of bound states of
solitons in nonlinear dissipative systems. In particular, a machine-learning
optimization algorithm, which helps to create a bound state of two solitons
with a required TS in the fiber laser, was reported very recently [22].

Passively mode-locked fiber lasers offer a highly efficient platform for the
generation of the bound states. Many experiments performed in such systems
have produced stable SMs with the fixed relative phase [4, 23, 15], and TS
corresponding to tightly or loosely bound “molecules" [4, 23, 6, 24, 16]. Such
SMs were predicted by theoretical models based on the nonlinear
Schrödinger/Ginzburg-Landau equations [25, 2, 3]. However, controllable
manipulations of SMs in a nonlinear dissipative system is a challenge to the
experiment. In particular, it is not easy to produce SMs with a predetermined
value of TS in a mode-locked fiber laser. Although several SMs with different
separations can be created by varying parameters of the fiber cavity, viz.,
the polarization or pump power, it is difficult to exactly adjust the TS, or
switch SMs between states with certain TSs “at will" [23, 6, 15, 7]. In
particular, the ability to produce bound states with predetermined (“on-
demand") values of the separation is important for SM-based optical data-
processing schemes, communications [19], optical switching [20], storage [26,
21], and soliton trapping [27]. Recently, the technique of active pump
modulation was adopted to efficiently manipulate SMs in mode-locked fiber
lasers [20, 28]. In particular, all-optical switching of two SMs was realized
by dint of the frequency-swept pump modulation [20]. Actually, the pump
modulation is a powerful technique to control the SMs, while pump
perturbations may lead to unwanted effects, such as generation of higher-order
harmonics [20].

In this paper, we propose and experimentally realize a scenario allowing one
to perform on-demand manipulations of photonic SMs in a passively mode-locked
laser. Different from the external pump modulation, our method relies on
directly varying parameters of the nonlinear passive system, viz., the group-
velocity dispersion (GVD) and dispersion losses, to switch the SM from one
state to another, with predetermined values of the TS. The method works via
controlling the complex dispersion: its real part determines the balance
between the GVD and nonlinearity in the mode-locked laser, while the imaginary
(dispersion-loss) part maintains the compensation of the frequency-dependent
losses by the gain. Studies of manipulations of the dispersion had a long
history in the course of the development of the dispersion management for
mode-locked lasers [29, 30, 31, 32, 33, 34]. In particular, a fourth-order
soliton was recently created by engineering both second- and fourth-order GVD
in a mode-locked laser system [32]. The result shows the imposed dispersion in
the cavity largely affects the nature of the soliton pulse, such as the actual
temporal profile shape.  However, the dispersion losses [3, 35, 18] have drawn
less attention in studies of nonlinear dissipative systems, especially as
concerns soliton pairs. Effects akin to the dispersion losses were used for
shaping temporal pulses in linear systems [36, 37]. Here, we demonstrate that
additional contributions to the dispersion losses play an important role in
the formation of the photonic SMs. In particular, in a specific range of
parameters, a large contribution to the imaginary part of the complex
dispersion increases the TS between the bound solitons. Experimentally, we
employ a $4f$ pulse shaper placed in the cavity which uses a hologram, to
tailor both GVD and the dispersion losses via encoding phase and amplitude
distributions in the hologram. The setup can realize producing multiple SMs
on-demand in the range of TS from $3.014$ to $5.478$ ps. By employing the
TSDFT technique, we observe the SM’s spectral dynamics and demonstrate the
realization of quaternary encoding and optical switching between four SMs with
different temporal sizes. Further, by slightly moving the hologram in the
spatial light modulator, the tunability of the wavelength is realized, and
essential hysteretic phenomena are revealed. Thus, the hologram with the
adjustable structure and position is the main tool that is used in this work
to realize efficient control of SMs in the fiber laser. The results not only
offer a robust way to manipulate SMs in the mode-locked fiber system, but also
supports a new tunable degree of freedom for the experimental control of
nonlinear dissipative systems. The multiple SMs reported in this work may find
essential applications to the design of “multi-alphabetic" communications,
optical switching, and storage also in molecular spectroscopy.

##  2 Methods for on-demand generation of SMs

Figure 1: The scheme for the generation of SMs (soliton molecules). (a): The
layout of the mode-locked fiber laser. Here, Pump stands for the pump source
at $980$ nm; WDM+ISO is the wavelength-division multiplexer for $980$ nm and
$1560$ nm combined with the isolator for $1560$ nm; SA is the saturable
absorber; PC1(2) are polarization controllers that includes one half-wave and
two quarter-wave plates; and OC is the output coupler. The pulse shaper system
includes two infrared optical gratings G1,2 (the grating period is 940 gr
mm-1), input and output lenses $F_{1,2}$(the focal length is 125 mm), and an
infrared spatial light modulator (SLM). It has the resolution of $1920\times
1080$ and the pixel pitch of 6.4 $\mu m$. (b): The output SMs built of two
solitons with the $\mathrm{sech}(t)$ profile, corresponding to different
values of the spectral filtering,
$\mathrm{Im}\left(\tilde{\beta}_{2,\mathrm{holo}}\right)$. Here,
$\ket{m}_{A,B}$ are symbols for SMs with different values of TS (temporal
separation between the bound solitons), while subscripts $A$ and $B$ represent
SMs with different central wavelengths. (c): Phase and amplitude distributions
of the encoding factor, given by Eq. (3) with $\tilde{n}=10\left(1+i\right)$
and $\delta\lambda=0$, in the case when only the second-order GVD and
dispersion losses are presented. The distributions are displayed versus the
offset of the operational wavelength from the central point,
$\Delta\lambda=\lambda-1554$ nm. (d): Same as in (c), but for the setup with
the central wavelength shifted by $\delta\lambda=2.5$ nm.

Figure 1 illustrates the principle of manipulating the photonic SMs determined
by the effects of GVD and dispersive losses. Panel 1(a) shows the layout of
the passively mode-locked laser, which includes fiber-guided and free-space
sections. The former one consists of a $3.3$ m long erbium-doped fiber
providing the necessary gain, a homemade saturable absorber providing the
mode-locking in the cavity, and the output coupler extracting $10\%$ of the
power from the fiber cavity. In the free-space section, we employ a $4f$ pulse
shaper to morph the amplitude and phase of the spectrum. The total length of
the cavity is $\approx 19.88$ m, which includes $1.1$ m in the free space for
spectrum shaping. The measured cavity’s fundamental repetition rate is
$10.458$ MHz. The net dispersion of the fiber cavity is $-0.2672$ ps2, i.e.,
the mode-locked laser operates in the anomalous GVD region.

The generation of SMs in the mode-locked fiber laser can be adequately modeled
by the complex Ginzburg-Landau equations (CGLEs) [2, 3, 38]. The appropriate
CGLE for amplitude $A\left(z,t\right)$ of the output optical pulse is [39, 32]

$\frac{{\partial A}}{{\partial
z}}=i\sum\limits_{k}{\frac{{{{\tilde{\beta}}_{k}}}}{{k!}}}{\left({i\frac{\partial}{{\partial
T}}}\right)^{k}}A+\frac{{g-\alpha}}{2}A+i\gamma|A{|^{2}}A+\mathrm{H.O.E}.$ (1)

Here, $z$ and $t$ are the propagation distance and local time; $\gamma$ is the
Kerr parameter; $g$ and $\alpha$ represent the linear gain and loss,
respectively; and $\mathrm{H.O.E}$ stands for terms representing higher-order
nonlinear effects, which includes the saturation of the nonlinear gain,
saturation of the nonlinear refractive index, etc. The higher-order terms are
essential for the normal operation of the system [18, 3]. In this model,
$\tilde{\beta}_{k}=\mathrm{Re}(\tilde{\beta}_{k})+i\,\mathrm{Im}(\tilde{\beta}_{k})$
includes GVD, $\mathrm{Re}(\tilde{\beta}_{k})$, and the dispersion losses,
$\mathrm{Im}(\tilde{\beta}_{k})$, of orders $k=2,3,4,...$. In this paper,
$\tilde{\beta}_{k}$ are referred to as $k$-th order complex dispersion
coefficients. In particular, $\tilde{\beta}_{2}$ represents the complex
second-order dispersion.  Its GVD part, $\mathrm{Re}(\tilde{\beta}_{2})$,
plays the central role in the operation of the dispersion management in mode-
locked lasers [29, 30, 31, 32]. However, the dispersion-loss part,
$\mathrm{Im}(\tilde{\beta}_{2})$, being a necessary ingredient for the
creation of two- and multi- soliton bound states in dissipative nonlinear
systems [2, 3, 40, 35, 18], was less explored experimentally, especially in
the context of SMs. In the traditional mode-locked fiber laser system without
the addition filter, it is related to the gain bandwidth $\Omega_{g}$, viz.,
$\mathrm{Im}(\tilde{\beta}_{2})=g/2\Omega_{g}^{2}$ [38]. In our system,
$\tilde{\beta}_{k}$ can be precisely adjusted by means of the $4f$ pulse
shaper, that transforms the temporal pulse into the frequency domain, and then
back into the temporal one. Using the spatial light modulator (SLM) installed
at the center of the pulse shaper, it is possible to perform flexible
management of the complex dispersion coefficient of an arbitrary order
(details of the pulse shaper could be found in Section 1 of the Supplementary
Material). In this paper, the second-order term is addressed in the main text,
the higher-order terms being considered in Section A of the Appendix and
Section 2 of the Supplementary Material. As mentioned above,
$\mathrm{Re}\left({\tilde{\beta}_{2}}\right)$ provides the dispersion which
balances the nonlinearity in the setup of the mode-locked laser, which is the
ordinary condition necessary for the formation of optical solitons [41, 42,
32], while the loss part, $\mathrm{Im}\left({\tilde{\beta}_{2}}\right)$, is
responsible for maintaining the balance of the gain and losses in the system
[43]. To distinguish the complex dispersion coefficients in different sections
of the system, we define $\tilde{\beta}_{2,\mathrm{fiber}}$ as the coefficient
in the fiber section, and $\tilde{\beta}_{2,\mathrm{holo}}$ as one imprinted
by the hologram in the $4f$ pulse shaper.

Our experimental findings and simulations (see Section A in Appendix)
demonstrate that the dynamics of the SMs can be tailored by adjusting the real
and imaginary parts of ${\tilde{\beta}_{2,\mathrm{holo}}}$. In this context,
Fig. 1(b) shows that TS between the bound solitons gradually increases,
following the increase of $\mathrm{Im}{(\tilde{\beta}_{2,\mathrm{holo}})}$ in
a certain range. In particular, only the conventional single soliton,
symbolized by $\ket{0}$, is possible at
$\mathrm{Im}(\tilde{\beta}_{2,\mathrm{holo}})=0$. With the increase of
$\mathrm{Im}(\tilde{\beta}_{2,\mathrm{holo}})$, the system can create SMs on-
demand with different values of TS, $\ket{m}(m=1,2,3...)$ (see the
illustration in Fig. 1(b)). In addition, the setup could change the color
(central wavelength) of SMs, while keeping the TS constant. The latter option
is shown in the bottom line of Fig. 1(b) by means of labels $\ket{4}_{A,B}$,
where $A$ and $B$ represent the central wavelength of two SMs. Experimentally,
both the GVD and dispersion-loss constituents of
$\tilde{\beta}_{2,\mathrm{holo}}$ can be manipulated through the phase and
amplitude distribution of the hologram in the SLM. The hologram’s
characteristic is defined as
$H(x)=\phi_{R}(\mathrm{Re}(\tilde{\beta}_{2,\mathrm{holo}}),x)\cdot\phi_{I}(\mathrm{Im}(\tilde{\beta}_{2,\mathrm{holo}}),x)\cdot\phi_{G}(d,x)$,
which includes three factors used for the implementation of the precise
manipulations. Here,
$\phi_{R}(\mathrm{Re}(\tilde{\beta}_{2,\mathrm{holo}}),x)$ and
$\phi_{I}(\mathrm{Im}(\tilde{\beta}_{2,\mathrm{holo}}),x)$ represent real and
imaginary terms (corresponding to the phase and amplitude encoding,
respectively), while $\phi_{G}(d,x)$ represents an optical grating that helps
to avoid crossover with the zero-order diffraction [44, 45]. For the
description of the experiment, we define the dimensionless complex coefficient
$\tilde{n}$ to encode the hologram by means of its real and imaginary parts:

$\displaystyle\mathrm{Re}\left({\tilde{\beta}_{2,\mathrm{holo}}}\right)$
$\displaystyle={k_{R}}\cdot\mathrm{Re}\left({\tilde{n}}\right)$ (2a)
$\displaystyle\mathrm{Im}\left({\tilde{\beta}_{2,\mathrm{holo}}}\right)$
$\displaystyle={k_{I}}\cdot\mathrm{Im}\left({\tilde{n}}\right).$ (2b) Here,
$k_{R}=-0.2170\times 10^{-3}$ ps${}^{2}/$m and $k_{I}=0.9309\times
10^{-3}\leavevmode\nobreak\ $ps${}^{2}/$m are two appropriate constants that
are selected empirically in accordance with the experiment; the values are
subject to the constraint imposed by the capacity of the pulse shaper (see the
details in the Section 1 of the Supplementary Material).
$\mathrm{Re}\left({\tilde{n}}\right)$ represents the positive or negative
number associated with GVD, while $\mathrm{Im}\left({\tilde{n}}\right)$
connects to the dispersion losses, hence it may be only positive. Following
these definitions and ignoring the higher-order dispersion (which is addressed
in Section 2 of the Supplementary Material), the encoding part of the hologram
is written as
${\phi}\left({{{\tilde{\beta}}_{2,\mathrm{holo}}}}\right)\simeq\exp\left[{i{{\tilde{\beta}}_{2,\mathrm{holo}}}(\tilde{n}){{L_{\mathrm{cavity}}(\omega-\delta\omega)}^{2}}/2}\right],$
(3)

where $L{{{{}_{\mathrm{cavity}}}}}$ is the total length of the cavity (as
mentioned above, it is $19.88$ m in our experimental setup), $\omega$
represents the optical frequency corresponding to the wavelength at the
spatial position of the SLM, and $\delta\omega$ is an offset of the central
frequency, which leads to a shift for the effective gain wavelength range of
the mode-locked fiber laser. The offset, which can be implemented
experimentally by slightly displacing the position of the hologram, provides
the precise tunability of the central wavelength for the output SMs.

Figure 1(c) displays an example of the distribution of the phase (blue) and
amplitude (red) produced by Eq. (3) with $\delta\omega=0$, where
$\tilde{n}=10\left(1+i\right)$ is set. Such one-dimensional curves can be
expanded for the two-dimensional hologram by adding the optical grating (see
Section 2 in Supplementary Material for details in hologram). Figure 1(d)
shows the distribution produced by Eq. (3) with $\delta\omega$ corresponding
to the wavelength shift $\delta\lambda=2.5$ nm.

In the present regime, the real and imaginary parts of
$\tilde{\beta}_{2,\mathrm{holo}}$ can be independently determined by the
amplitude and phase structure of the hologram. The ability to choose these
parameters offers a robust scheme for manipulating SMs in the mode-locked
system. In turn, the output produced by the system, in the form of the
photonic SMs, offers the platform for encoding data in the optical domain,
including options for multiple encoding and optical switching [19, 20]. Below,
we report experimental results for manipulations of TS in the bound soliton
pairs, quaternary encoding, and switching between several SMs, as well as the
wavelength tunability of SMs and a hysteresis effect for them.

##  3 Results

###  3.1 Manipulations of the TS (temporal separation) for SMs (soliton
molecules) in the mode-locked laser

Figure 2: The generation of the SMs (soliton molecules) in the experiment. (a)
and (b): The spectral intensity in the single-soliton ($\tilde{n}=0$) and SM
($\tilde{n}=10\left(1+i\right)$) regimes, respectively. The inserts in the top
left corner are the applied amplitude profile of hologram in SLM,
respectively.  (c) and (d): Spectral intensities measured by the optical
spectrum analyzer and their Fourier transform, which is also the
autocorrelation function, obtained for Im$\left(\tilde{n}\right)$ varying from
$0$ to $30$. (e): The evolution of spectral intensities measured by TSDFT,
where the recording length is $100$ round trips; the corresponding
autocorrelation functions are shown at the bottom. Here, the insets at the
bottom of the panels show oscillations of differences in the energy and
separations (yellow lines), and green horizontal lines represent limit values.
The basic separation is $\pm$ 100 fs. (f): The TS (temporal separation between
the bound solitons) extracted from the spectrum recorded by TSDFT. Here, we
have added the triple standard deviation to each data point. For instance, it
is $66$ fs for $\tilde{n}=10(1+i)$. Circles surrounding the points are added
for better visualization. (g): Results produced by simulations of the
theoretical model based on the CGLE and pulse-shaping equations (see Section A
in Appendix). The inset in the top left corner is the spectral intensity
produced by the simulations in the single-soliton regime, where the bandwidth
is $4.37$ nm.

We first adjust the mode-locked fiber laser to work in the regime of producing
individual solitons. To this end, first, the hologram parameter is taken as
$\tilde{n}=0$, which represents a simple transmitting mirror. Next, we
optimize parameters in the cavity, viz., the pump power and polarization, to
provide for the output in the form of a single soliton. The output spectrum of
the single-soliton is shown in Fig. 2(a), where the bandwidth is $4.40$ nm.
Two strong spectral sidebands observed in the spectrum are the usual feature
of the conventional single-soliton generation regime [46, 47, 17]. Using a
homemade autocorrelation-FROG system, we reconstruct the solitons’ temporal
and spectral profiles and their envelope phase (see Section 3 in Supplementary
Material). The so reconstructed pulse width is $670$ fs, and its chirp [41] is
very low.

Next, we switch the hologram, changing the complex control parameter,
$\tilde{n}$, defined above in Eqs. (2a) and (2b), so as to provide the
generation of soliton pairs. The SMs could be observed from spectral fringes,
that, in turn, are produced by the optical spectrum analyzer [9, 6]. We set
$\mathrm{Re}(\tilde{n})=10$ and increase $\mathrm{Im}(\tilde{n})$, starting
from $0$. The output spectra are unstable in the interval from $0$ to $6$,
which means many possible different spectra (an example is presented in Video
1 of Supporting Material). A stable SM appears when
$\mathrm{Im}(\tilde{n})=6$. Later, the output features SMs, whose TS increases
linearly up to $\mathrm{Im}(\tilde{n})\approx 29$. Fig. 2(b) shows a power
spectrum representing the SM for $\tilde{n}=10\left(1+i\right)$, where the red
and blue lines represent, respectively, the experimental result and fitting to
it (for the fitting method, see Section 4 in Supplementary Material). Figure
2(c) shows the spectral distribution for different values of
$\mathrm{Im}(\tilde{n})$, _viz_., $0$ and $6\leq\mathrm{Im}(\tilde{n})\leq 30$
with interval $\Delta\left(\mathrm{Im}(\tilde{n})\right)=2$. The color-coding
in this figure represents the relative spectral intensity, which is normalized
to take values from $0$ to $1$.

The spectral distributions provide the basic information for the SMs in the
temporal domain. In particular, the TS of the bound solitons is inversely
proportional to the spectral fringe period, and their relative phase can be
inferred from the envelope of the spectral interferogram [6]. To directly find
the TS, it is instructive to perform the Fourier transform of the spectral
intensity, which also represents the autocorrelation function [9, 6]. Figure
2(d) exhibits the Fourier transform of the corresponding autocorrelation
function. The results show that the TS grows linearly with an increase of
$\mathrm{Im}(\tilde{n})$ from $6$ to $28$. With the further increase of
$\mathrm{Im}(\tilde{n})$, the SM becomes unstable, switching to another SM
state with a large TS, as observed in Fig. 2(d) for
$\mathrm{Im}(\tilde{n})=30$. This phenomenology is similar to that
demonstrating the Hopf-type bifurcation in two-soliton bound states in other
systems [48, 5, 18]. From the spectral interferogram of SMs, one can find that
the dispersion loss strongly affects the separations, while it is practically
not related to the relative phase between two soliton pairs. For the relative
phase, many previous works show that it may be affected by the pump power or
PC orientation in the mode-locked fiber system [6]. The full dynamical
behavior of SMs supported by different holograms can be seen Video 1 in
Supporting Material.

To explore the stability of the SMs, we measured the single-shot spectrum by
means of the TSDFT technique, which may also be used to perform real-time
multiple SMs encoding and switching (the traditional optical spectrum analyzer
is unsuitable in this case, due to the slow scanning speed). For this purpose,
the output optical solitons are injected into a dispersion fiber ($5000$ m
long) [10, 49], in which the spectral profile is mapped into the time domain
[10, 50] (Also see the Section 1 of the Supplementary Material). The recorded
dispersive Fourier-transform spectra are shown in Fig. 2(e), where we present
the conventional single soliton (at $\mathrm{Im}(\tilde{n})=0$) and three SMs
(at $\mathrm{Im}(\tilde{n})=10,26,30$) with different STs. In these
measurements, the recording length (shown on the horizontal axis) is up to
$100$ round trips. Four subplots displayed in the figure are the corresponding
autocorrelation functions. Figure 2(f) shows the average TS extracted from the
dispersive Fourier-transform spectra, where the blue line is the fitting
result, produced by the linear model in the following form: TS(ps)
$=0.112\times$ $\mathrm{Im}{(\tilde{n})}+2.342$,
$6\leq\mathrm{Im}(\tilde{n})\leq 28$. Close to $\mathrm{Im}(\tilde{n})=29$, a
large jump in the TS appears, and then it does not change conspicuously with
the further increase of $\mathrm{Im}(\tilde{n})$. Here, only the range of
$\mathrm{Im}(\tilde{n})$ between $6$ and $28$ is considered, where the TS can
be tailored at will from $3.014$ to $5.478$ ps, and the SMs can be switched
gently without a hysteresis. Such features are essential and useful to perform
quaternary encoding and optical switching in the SMs, as shown in the next
section.

One can analyze the stability of the SMs from the single-shot spectrum
recorded by TSDFT. There are several methods to character the stability of the
SM states. In particular, the stability can be explored in the plane of the
separation and relative phase [51, 3, 6], through measuring the radio-
frequency (RF) spectra [52, 53], or directly observing variations of the
energy and separation of the bound pulses [3, 54, 55, 56, 57]. By considering
our demonstrations, we show the variations in pulse energy and separations
from the single-shot spectrum, which could directly reflect the stability of
SMs. For single-shot spectrum, we define the energy difference as $dQ_{m}/dm$,
where
${Q_{m}}=\frac{1}{{{Q_{ave}}}}\int{{{\left|{A(m,\omega)}\right|}^{2}}d\omega}$
is the spectrum energy for $m$-th round trip from TSDFT, and
${Q_{ave}}=\sum\limits_{m=1}^{M}{\frac{1}{M}\int{{{\left|{A(m,\omega)}\right|}^{2}}d\omega}}$
is the average value among all the recording $M$ round trips. Therefore, the
difference could be simply given as $Q_{m}-Q_{m-1}$. Actually, the definition
is the stable criterion for SMs in energy aspect [3]. We calculated the energy
difference for our recording single-shot spectrum of $\mathrm{Im}(\tilde{n})$
=0, 10, 26, and 30, which are shown at the bottom of Fig. 2 (e). Here, two
green lines set the basic energy level, which is $\pm$10% for $Q_{ave}$; the
yellow lines give the oscillations of the energy difference. The similar
calculation is performed to the separations that simply defines
$\tau_{m}-\tau_{m-1}$ . We add the separation difference at the below of
autocorrelation functions. Here, two green basic lines set $\pm$100 fs, which
correspond to the same level of the temporal resolutions from the TSDFT system
(see the Section 1 of the Supplementary Material ). The yellow line represents
the extracting difference, which is slightly below the set level $\pm$100 fs.
These results illustrate the SMs ($\mathrm{Im}(\tilde{n})=0,10,26,30$) are
stable at the current recording range (100 round trips). It is enough for the
demonstrations because the multiple SMs based encoding regime only accounts
for 20 round trips. By the way, in the Section 1 of the Supplementary Material
, we add one radio-frequency (RF) spectra for SM state of $\tilde{n}=10(1+i)$.
It shows the average signal-noise level could be up to around 50 dB.

Figure 3: The SM-based quaternary encoding and switching. (a) The measured
single-shot spectrum for string $\left\\{0|3|1|2|0|3\right\\}$, produced by
the TSDFT technique with $20$ round trips. (b): The corresponding TSs of six
SM states. (c)-(d): A $\psi$-shaped pattern in the binary and quaternary
encoding systems, respectively. (e): The measured single-shot spectrum for the
$\psi$-shaped pattern with 36 encoding holograms in the string of
$\left\\{0000|0000|0320|0020|1333|0020|0320|0000|0000\right\\}$, where each
encoding sequence is produced by TSDFT with $10$ round trips. (f): The
corresponding TSs extracted from the encoding sequence with 36 holograms.

We have performed simulations of the dissipative nonlinear system shown in
Fig. 1(a). For this purpose, we split the cavity into two main parts: in the
fiber section, it is modeled by CGLE taken as Eq. 1, while in the free-space
section the model is provided by the pulse-shaping equation. Also, we need to
address equations for the saturable absorber and gain from the erbium-doped
fiber. Details are given in Section A of Appendix. Parameters were first
adjusted to those corresponding to the single-soliton regime in Fig. 2(a). The
result obtained by the simulations for the single-soliton spectrum is shown in
the inset to Fig. 2(g), where two sidebands obviously belong to the spectral
domain. The respective bandwidth is $4.37$ nm, which is in good agreement with
the experimental observations ($4.40$ nm in Fig. 2(a)). Next, we set
$\tilde{n}=10+6i$, which leads to the creation of bound soliton pairs. A
similar simulation, for $\tilde{n}=10\left(1+i\right)$, is presented in Fig.
S1(b) of Section A in Appendix. Fig. 2(g) shows the autocorrelation function
of the output produced by the simulations after 2000 round trips (the
corresponding spectra are displayed in Fig. S1(c) of Section A in Appendix).
From these results, it is seen that the TS of the bound states increases near
linearly with the growth of $\mathrm{Im}(\tilde{n})$ in the range of
$\mathrm{Im}({\tilde{n}})$ from $6$ to $28$, similar to what is observed in
the experiments, cf. Fig. 2(d). There are two conspicuous differences between
simulations (Fig. 2(g)) and experiments (Fig. 2(d)). One is the variation for
the separations: the experiments show a linear increase of TS, while the
simulations produce a nonlinear change. The other is the values of TS: the
predicted and experimentally observed separations are not exactly equal for
each Im$(\tilde{n})$. These differences between the experiments and
simulations are due to minor differences in the GVD and dispersion loss. They
mainly originate from imperfections of the $4f$ pulse shaper, as concerns the
coupling efficiency, encoding formats, and imaging precision (Section 2 in
Supplementary Material addresses consequences of the imperfect encoding for
the hologram). The imperfections could lead to the appearance of higher-order
dispersion and dispersion losses in the experimental setup. If the higher-
order term is added to the theoretical model, the difference between the
experimentally observed and numerically predicted values of the TS becomes
smaller. For example, recent theoretical results have demonstrated that the
third-order dispersion strongly affects the separations of SMs [18].  Further,
the imperfect $4f$ pulse shaper can also limit the tuning range of separations
of SMs in the present system. Previous studies have also found that the linear
tuning range was limited – for example, in the presence of the Hopf
bifurcation [5, 18]. If, as a result, the separation undergoes a jump with the
variation of the system’s parameters, this implies instability or hysteresis
of the SM. In additional simulations, which include the TOD term (not shown
here in detail), we found that coexisting bound states appear, especially for
the larger value of Im$(\tilde{n}$), viz., Im$(\tilde{n})>35$. This fact may
quantitatively explain the jump that happens at larger Im$(\tilde{n})$ .
Generally, the theoretical model proposed here is not comprehensive for the
comparison to the experiment. Indeed, not all experimental parameters could be
exactly determined, and some assumptions are adopted to run the simulations.
Nevertheless, the main features of the experimental results and their
dependence on parameters are correctly predicted by the present model.

###  3.2 SM-based quaternary encoding and optical switching

Multiple encoding scenarios based on the use of SMs had been proposed long
ago, in the context of the multiple-soliton communications [19], storage [21],
soliton trapping [27], and optical switching [20]. Multiple encoding helps to
increase the information capacity in these applications. However, there were
few experiments to demonstrate multiple encoding. The main issue is the
difficulty in the preparation of an “on demand" SM source. Recently, the
creation of two deterministic SMs was demonstrated with the help of the pump-
modulating technique [20]. A challenging problem is to efficiently manipulate
more SMs and demonstrate possible regimes of multiple encoding and switching.
The SM regimes revealed in this work offer new possibilities to tailor the
bound states of optical solitons, which can be directly applied to the design
of SM-based quaternary encoding and optical-switching schemes [19, 20].

Multiple SMs with different values of TS offer a possibility for multiple
encoding using different states $\ket{m}$, as shown in Fig. 1(b). Indeed, it
is seen from Figs. 2(b) and (d) that the system can produce many SM varieties.
For example, the TSs in the range of $3.014$ to $5.478$ ps correspond to
$\mathrm{Im}\left(\tilde{n}\right)$ from $6$ to $28$. To avoid the crosstalk
induced by the interaction between adjacent SMs [58] and the TS jitter during
the hologram switching, we select four SMs, with
$\mathrm{Im}\left(\tilde{n}\right)=10,14,18$, and $22$, as states
$\\{\ket{0}$, $\ket{1}$, $\ket{2}$, $\ket{3}\\}$, respectively, for the
realization of quaternary encoding. Experimentally, these four states can be
controlled by appropriately using four holograms in the SLM; thus, the optical
switching can be performed by replacing the corresponding hologram. In the
experiment, we use four SMs to create a quaternary encoding scheme that forms
four optical bits. Here, we do not use the single-soliton regime
($\mathrm{Im}(\tilde{n})=0$) and SMs with a very large separation
($\mathrm{Im}(\tilde{n})\geq 30$). The main reason for the exclusion of such
states is the presence of strong hysteresis between them, which would disrupt
the implementation of the multiple-encoding and optical switching. Similar
hysteretic phenomena readily appear in the work with polarization- or power-
induced SMs [5]. However, our experimental results reveal no evident
hysteresis in the above-mentioned range of $\mathrm{Im}\left(\tilde{n}\right)$
between $6$ and $28$.

Thus, we start by using the four above-mentioned SMs, $\ket{m}$ with
$m=0,1,2,3$, and construct a simple test string
$\left\\{0|3|1|2|0|3\right\\}$. The respective set of six holograms are
consecutively displayed on the SLM, switching them automatically in sequence
with the time interval of $\delta t=$ $511$ ms, limited by the low operating
speed of the liquid-crystal SLM. Simultaneously, TSDFT records the
corresponding single-shot spectrum. In particular, the oscilloscope records
and stores four single-shot spectra for each hologram, in which the sampling
number is limited by the storage depth of the oscilloscope. Fig. 3(a) shows an
example of the spectrum, with $20$ round trips recorded for each hologram.
Further, Fig. 3(b) shows the corresponding values of TS in the bound state of
two solitons, as extracted from the single-shot spectrum. We conclude that the
single-shot spectrum is nearly stable in the framework of each encoding
sequence, and the corresponding TSs feature an acceptable jitter.

To perform quaternary encoding by four states $\ket{m}$ with different values
$\tau$ of the TS, a robust criterion should be adopted to distinguish values
of $\tau$. Here, we sort them into four states: $\ket{0}$
($\tau\in[3.0,3.64]$), $\ket{1}$ ($\tau\in[3.65,4.06]$), $\ket{2}$
($\tau\in[4.07,4.48]$), and $\ket{3}$ ($\tau\in[4.49,5.00]$), respectively.
Labels attached to the right vertical axis in Fig. 3(b) shows the respective
numbers of the so-coded states. According to the TS values in Fig. 3(b), the
output state is indeed $\left\\{0|3|1|2|0|3\right\\}$, being tantamount to the
encoding sequence which the procedure aimed to create.

We then consider a more complex encoding, viz., the $\psi$-shaped binary graph
is shown in Fig. 3(c), which includes $72$ digits in the binary format
($00000000|00000000|00111000|00001000|01111111|00001000$
$|00111000|00000000|00000000$) However, using the SM encoding, it may be
reduced to a set of $36$ digits in the quaternary-encoding basis, as shown in
Fig. 3(d), namely,
$\left\\{0000|0000|0320|0020|1333|0020|0320|0000|0000\right\\}$. We generate
the corresponding set of $36$ holograms and switch them automatically in a
sequence similar to how this was done above for the simple string,
$\left\\{0|3|1|2|0|3\right\\}$. Fig. 3(e) shows the single-shot spectrum
recorded with the help of TSDFT, using $10$ round trips for each state
sequence. The criterion to distinguish the states $\ket{m}$ with $m=0,1,2,3$
is the same as adopted in Fig. 3(b). Analyzing the autocorrelation function,
we could identify the corresponding values of the separation of SMs, which is
presented in Fig. 3(f), with the output pattern,
$\left\\{0000|0000|0320|0020|1333|0020|0320|0000|0000\right\\}$, being exactly
the same as the encoding sequence defining the $\psi$-shaped target. The
dynamic switching results for the $\psi$-shaped patterns in the quaternary-
encoding regime is presented in Supporting Material for Video 2.

Figure 4: The wavelength tunability (for $\tilde{n}=10+30i$) and hysteresis
(for $\tilde{n}=10+19i$) of the SMs. (a): The spectral fringe for SMs with
different central wavelengths, while their TS (temporal separation) is kept
nearly constant, with the avergae value $\tau=12.544$ ps, and average
variation $0.288$ ps. The corresponding autocorrelation function is shown in
the right bottom corner. (b): The central wavelength $\lambda_{0}$, extracted
from the experimental data, versus the relative position $\delta x$ of the
hologram. Here, we have added the measurement error $\pm 0.1134$ nm to each
data point, which is the wavelength resolution at $\sim$ 1555 nm of the
employed optical spectrum analyzer. Panels (c) and (d) display a wavelength-
controlled hysteresis for the SMs, in terms of the TS. (c): The
autocorrelation function of the SMs, where the top and bottom plots correspond
to varying relative $\delta x$ from left (shorter wavelength) to right (longer
wavelength), or from right (longer wavelength) to left (shorter wavelength),
respectively. (d) The TS data extracted from (c). In these results, $\delta
x=1$ means that the actual spatial shift of the hologram in SLM is 10 $\mu$m.

The demonstrated SMs can realize multiple encoding based on the set of four
states in the quarternary system and switch between them with high fidelity.
Compared to the binary encoding system with $2^{n}$ options, the quaternary
one offers $4^{n}$, which is beneficial to perform dense coding in optical
communications, imaging, and computations [59, 60]. In our system, if one
suitably reduces the criterion interval distinguishing different bound soliton
pairs, the system could realize more SMs ($d$), which means it’s possible to
perform higher-dimensional ($d^{n}$) encoding and switching. In this case, a
more accurate technique should be employed to characterize the separations of
SMs, such as the balanced optical cross-correlation method [61], to show their
time jitter on the fs scale.

###  3.3 The tunability of the central wavelength of the SMs and the
hysteresis

The ability to tune the wavelength is needed in many applications, including
optical sensing, metrology, and spectroscopy [34, 62]. The usual tuning
methods employ the nonlinear wave mixing [63], or the use of high-speed
electric optical modulators. In a mode-locked fiber laser system, one usually
changes the soliton’s wavelength by adjusting the intra-cavity tunable filter
[64, 65]. For SMs, there was no report of high-precious tuning their
wavelengths. The present setup may offer a simple platform to adjust the
wavelength of SMs. According to Eq. (3), $\delta\omega=\left(2\pi
c/\lambda_{0}^{2}\right)\delta\lambda$ represents the shift with respect to
the central wavelength $\lambda_{0}$, which leads to a shift of the effective
gain area in the cavity, and can be used to realize the precise tunability. In
the experiment, by shifting the hologram in the $x$ direction by $\delta x$,
which affects the real and imaginary parts of GVD in SLM through the space-
frequency relation of the $4f$ shaping system [36], the effective gain area
displaced by $\delta\omega$ will drag the central wavelength of the
intracavity pulse (see Fig. 1(d)).

Figure 4 displays the experimental results concerning the shift of the central
wavelength of the SMs. We first lock the fiber laser to a stable SM operation
corresponding to $\tilde{n}=10+30i$, for which the central wavelength is
$1554.6$ nm. Then, after moving the position of the hologram in SLM by $\delta
x$, the corresponding normalized spectrum is recorded by the optical spectral
analyzer, as shown in Fig. 4(a). We further apply the Fourier transform to the
spectrum, to produce the autocorrelation functions, which are shown in the
inset to Fig. 4(a). These results demonstrate that the SM’s central wavelength
changes linearly, while the period of the spectral fringes remains nearly
constant, corresponding to the TS which takes values $12.544$ $\pm$0.288 ps.
Figure 4(b) shows the relation between $\delta x$ and the wavelength shift
extracted from Fig. 4(a), where red circles represent data of experimental
measurements, and the blue solid line is the fitting produced by the linear
model, namely, $\lambda_{0}({\mu m})=0.0455\times\delta x+1.5546$. Thus, our
results exhibit the wavelength shift by $4$ nm (from $1555$ nm to $1559$ nm),
with the linear dependence on $\delta x$.  However, the linear relationship
cannot be valid in the entire range of the variation of $\delta x$. On the one
hand, the gain profile in the amplifying fiber is not as linear as the
function of the wavelength [66]. On the other hand, a larger variation of the
wavelength may lead to a large jump, such as in the case of the hysteresis
observed in Fig. 4(c). Nevertheless, it is reasonable to use the linear model
to fit the tunability in the small wavelength-tuning range. Therefore, the
present regime may offer an effective method to accurately tune the wavelength
of SMs.

In the experiment, we have also found that TS of the bound solitons features a
large jump following further increase or decrease of $\delta x$. However, the
original TS is not exactly recovered when $\delta x$ returns to the original
value, i.e., hysteresis occurs. As mentioned above, it usually takes place
between power- or polarization-dependent states in mode-locked lasers [35,
67]. Our results reveal a different phenomenon in terms of the TS in the SM
states, viz., the wavelength-dependent hysteresis. Figure 4(c) shows measured
autocorrelations of SMs, with $\tilde{n}=10+19i$, as functions of $\delta x$.
The top plot shows the situation corresponding to varying $\delta x$ from
$-10$ to $+10$, while TS $\tau$ jumps from $4.7$ to $3$ ps around $\delta
x=2$. The inverse jump happens around $\delta x=-2$, when $\delta x$ varies
from $+10$ to $-10$, as shown in the bottom plot of Fig. 4(c). Further, Fig.
4(d) shows the TS dynamics corresponding to Fig. 4(c), with blue circles and
magenta stars representing, respectively, the data collected in the experiment
with $\delta x$ varying from left (shorter wavelength) to right (longer
wavelength) and vice versa. Summarizing these observations, we conclude that
two different bound states with large and small values of TS coexist in the
middle green area, depending on the direction of the variation of $\delta x$.
Such bistability areas are key elements used in optical logic devices [68,
69].

Geneally, hysteresis is an interesting phenomenon in mode-locked fiber laser
systems. For example, the hysteresis can be induced by variations of the power
or polarization [67, 35], and by interactions with dispersive waves [70],
while our system demonstrates the wavelength-controlled hysteresis. We
conclude that the SM-based hysteresis, found here in terms of the TS and
controlled by the shift of the central wavelength, may be more robust than the
conventional optical hysteresis controlled by the pulse’s power, which is
prone to instabilities breaking the power balance [6, 7, 71, 72].

The wavelength tunability of the SMs and their hysteresis are promising for
the development of the optical encoding and logic setups. For instance, the SM
may be selected as $\ket{m}_{A}$ with wavelength $\lambda_{A}$, and as
$\ket{m}_{B}$ with another wavelength $\lambda_{B}$, both having the same TS
(Fig. 1(b)). This setting will support a larger variety of SMs in the
nonlinear dissipative system, thus expanding the available multiplicity of the
SM-based encoding and switching. Therefore, it may find applications to
optical massive data processing [73]. Further, the possibility to control the
bistability by precise adjustment of position $\delta x$ of the hologram in
the free space may be appropriate for the design of optical free-space
computational processing [74].

##  4 Discussion

The objective of this work is to propose an effective scheme for manipulations
of the SMs (soliton molecules) built as bound states of two optical solitons
in a nonlinear dissipative system, i.e., the mode-locked fiber laser. The
experimental results, in agreement with systematic simulations of the laser’s
model, based on the nonlinear equations of the CGLE types, are summarized as
functions of the real and imaginary parts of the effective GVD of the laser
cavity. The experimental findings clearly demonstrate that the management of
the dispersion loss is an efficacious tool for tailoring bound states in the
present fiber-laser system. Namely, it allows us to preselect the TS (temporal
separation) between the bound soliton in the range of ([$3.014$,$5.478$] ps,
with the resolution of $0.112$ ps). The same tool makes it possible to
accurately tune the central wavelength of the bound state, and control the
hysteresis in terms of the wavelength. The numerical simulations reported in
this work also show that the dispersion loss is an effective control
parameter, which can be used to accurately select the value of the TS. It is
worthy to note that this parameter is also an important one on other models
based on CGLEs [2, 3, 18, 38]. An explanation is that a slight variation of
the dispersion loss affects the potential of the interaction of the two
solitons, which leads to a shift of positions of the interacting solitons, and
may trigger jumps between bound states with different values of the TS. This
phenomenology is somewhat similar to the recent prediction of transitions
between different bound states, which includes a possibility of hysteresis
created by the third-order dispersion [18].

The results reveal not only the new degree of freedom in the experiments with
nonlinear dissipative systems in photonics, but also allow encoding data in
the form of multiple stable SMs with different TSs, and SMs with fixed TS but
different colors (central wavelengths). Such multiple states may be toggled in
a controllable fashion, which is important for applications of SMs to
communications [19], optical switching [20], and storage [26, 21]. Also, these
findings offer potential for the use of the optical SMs in high-capacity
schemes for optical data processing [73, 74], machine learning based on
ultrafast photonics [75, 33, 34, 76], and operations with high dimensional
temporal quantum information [77, 78].

The results demonstrate that the hologram-driven mode-locking is a robust and
accurate technique for the work with the optical SMs. The switching speed of
the present setup, $\sim$ 2 Hz, is not high, being limited by the necessity to
refresh the liquid-crystal SLM (spatial light modulator). The scheme may be
improved by using faster SLM or high-speed amplitude (or phase) modulators,
such as ones based on digital micromirror devices, that secure the speed in
the range of kHz [79]. It is also possible to employ the acousto-optic
modulator that typically operate in the $\sim 100$ MHz range [37]. In the
latter case, it is possible to monitor the dynamics of switching between
different SMs.

##  5 Appendix

###  5.1 Simulation of the soliton molecules

The possibility of the existence of soliton bound states in the mode-locked
fiber system can be predicted by the model based on the complex Ginzburg-
Landau equation (CGLE) [2, 3, 38], which is written as Eq. (1) in the main
text. Simulations of the CGLE equation produce results for SMs (soliton
molecules) similar to the experimental findings. To provide still better
proximity of the numerical results to the experimental ones, we performed
simulations of a more detailed model, based on CGLEs in the fiber sections of
the system, the pulse-shaping function in the free-space section, and
saturation transfer function of the carbon-nanotube-based absorber.

The setup of the mode-locked fiber laser is shown in Fig. 1 of the main text.
It includes two fiber segments (the single-mode one and the doped gain fiber),
which provide GVD and gain. Accordingly, the model includes CGLEs and the
additional ingredient, defined as follows:

${\;\frac{{\partial A}}{{\partial
z}}=-i{\frac{{{{\tilde{\beta}}_{2,\mathrm{fiber}}}}}{{2}}}{{{\frac{\partial^{2}}{{\partial
T^{2}}}}}}A+\frac{g}{{2\Omega_{g}^{2}}}\frac{{{\partial^{2}}}}{{\partial{T^{2}}}}A+\left({\frac{g}{2}+i\gamma|A{|^{2}}}\right)A.}$
(A1)

In the section representing the pulse shaper, the underlying equations can be
averaged and added to Eq. A1 by addressing high order dispersions [32].
However, we find that the model is made more accurate if it explicitly
includes the pulse-shaping equation. Therefore, the free-space element (the
$4f$ pulse shaper [36]) that connects to the complex dispersion coefficient,
is modeled by:

${A(\omega,z+{L_{4\mathrm{f}}})=A(\omega,z)\cdot\sqrt{\delta}\exp\left({i\sum\limits_{k}{\frac{{{{\tilde{\beta}}_{k,\mathrm{holo}}}{L_{\mathrm{cavity}}}{\omega^{k}}}}{{k!}}}}\right).}$
(A2)

Here, ${{L_{\mathrm{cavity}}}}$ is the overall length of the cavity,
${{L_{\mathrm{4f}}}}$ is the overall length of the $4f$ pulse shaper, $\delta$
is the internal linear transmission rate that mainly originates from the SLM
(spatial light modulator), optical grating, and free-space-to-fiber coupling.
The coefficient $\tilde{\beta}_{k,\mathrm{holo}}$ represents, as said above,
the additional $k$-order complex dispersion. Recall that
$\tilde{\beta}_{2,\mathrm{fiber}}$ in Eq. (A1) represents the complex
dispersion of the fiber, ignoring the higher-order terms.

To realize the stable mode-locking regime, the setup includes the homemade
saturable absorber based on carbon nanotubes, cf. Refs. [80, 17]. It is
characterized by the transfer function,

$T_{\mathrm{ab}}=\sqrt{1-(\alpha_{\mathrm{ns}}+\alpha_{0}/(1+A^{2}/I_{\mathrm{sat}}\cdot
A_{\mathrm{eff}}))}.$ (A3)

Here $\alpha_{\mathrm{ns}}$ is the unsaturated absorption, $\alpha_{0}$ is the
linear limit of the saturable absorption, $I_{\mathrm{sat}}$ is the saturation
intensity, and the effective fiber’s cross-section area is
$A_{\mathrm{eff}}=63.6\times 10^{-12}\leavevmode\nobreak\ $m2.

Using Eqs. (A1)-(A3), we simulated the evolution of the soliton in the cavity.
In particular, CGLE given by Eq. (A1) was used to simulate the propagation
through the fiber; Eq. (A2), with the appropriate GVD coefficient, described
the action of the $4f$ pulse shaper in the frequency domain; and the transfer
function defined as per Eq. (A3) determined the passage of the pulse through
the saturation absorber, $A\rightarrow A\cdot T_{\mathrm{ab}}$. Coefficient
$g$ in Eq. (A1) is the gain parameter for the erbium-doped fiber, which is
determined by the small-signal gain $g_{0}$ and saturation energy
$E_{\mathrm{sat}}$:

$g=g_{0}\exp(-E(z)/E_{\mathrm{sat}}),$ (A4)

where $E(z)=\int_{-\infty}^{+\infty}|A(z,t)|^{2}dt$ is the total energy of the
pulse traveling in the cavity.

The element of the setup described by Eq. (A2) performed pulse shaping in the
time domain by means of the hologram encoding the frequency-domain features.
The optical grating transforms the input from the time domain into its
counterpart in the frequency domain. Then, the convex lens focused different
frequencies at different spatial locations, thus introducing the spatial
dispersion in the Fourier plane. With the help of the conjugate system, the
pulse was then transformed from the spatial domain back into the temporal one.
The central frequency of the $4f$ pulse shaper, along with the surface of the
SLM, follows this simple relationship based on the diffraction equations [37,
36]:

$\frac{{\partial\lambda}}{{\partial
x}}=\frac{{d_{\mathrm{grating}}\cos({\theta_{\mathrm{diff}}})}}{{mf}}$ (A5)

Here, $m$ is the diffraction order (usually $m=1$), $\theta_{\mathrm{diff}}$
is the input angle of the optical grating (here, it is $47.5^{\circ}$),
$f=125\leavevmode\nobreak\ $mm is the focal length of the input lens, and
$d_{\mathrm{grating}}$ is the grating constant, which is $940$ per mm in our
experiment. Following these definitions, the corresponding spatial dispersion
of the SLM is expected to be $5.75$ nm$/$mm, while, in our experiment, the
measured spatial dispersion is $5.2$ nm$/$mm. Specifically, we set a double-
slit hologram in the SLM to calibrate the actual spatial-dispersion
relationship $\partial\lambda/\partial x$. The spatial separation between two
slits corresponds to the frequency (wavelength) distance between two peaks
from the output spectrum. In the simulations, we also needed to consider the
filtering effect of the SLM’s surface size, where we set the spectrum
bandwidth to be $45$ nm. Due to the consideration of this term, we ignored the
gain bandwidth of $\Omega_{g}$ in Eq. (A1).

Figure S1: The absorption curve for the carbon nanotube saturable absorber (a)
and evolution of the output pulse (b) produced by the fiber cavity. In panel
(a), the fitting parameters are $\alpha_{0}=0.095$,
$\alpha_{\mathrm{ns}}=0.4908$, $I_{\mathrm{sat}}=21.6412$ W/m2. In panel (b),
the SM’s output intensity is displayed, for $\tilde{n}=10(1+i)$; the
horizontal axis represents the number of round trips from $1$ to $2000$; the
vertical coordinate shows the temporal coordinate from $-30$ to $+30$ ps.
Here, the separation between the bound solitons takes a stable value,
$\tau=5.67$ ps. The bottom portion of the panel (b) shows the difference in
the pulse energy between the current and last round trips. The top inset shows
details for the evolution between $1$ and $200$ round trips, which corresponds
to the plot shown in the bottom portion. Panel (c) shows all normalized
spectra intensities of the final outputs, where the horizontal axis represents
the sequence of values of $\mathrm{Im}(\tilde{n})$, and the vertical axis is
the wavelength, the central one being $1553$ nm.

The parameters used in the simulations were chosen as close as possible to
their experimental counterparts. The length of the erbium-doped fiber is $3.3$
m, with $\mathrm{Re}(\tilde{\beta}_{2,\mathrm{gain-fiber}})=20$ ps${}^{2}/$km
and $\gamma=0.0032$, the length of the single-mode fiber is $15.48$ m, with
$\mathrm{Re}\left(\tilde{\beta}_{2,\mathrm{single-mode-fiber}}\right)=-21.6$
ps${}^{2}/$km and $\gamma=0.0013$, and the net dispersion of the entire fiber
is $-0.2672\leavevmode\nobreak\ $ps2. Here, as said above, we ignore the
higher-order complex dispersion in the fiber section. The length of the free-
space path is $1.1$ m, hence the total length of the cavity is $19.88$ m. In
our setup, the single-mode fiber is placed in front of the $4f$ pulse shaper,
therefore the simulations were performed in two main parts of the setup, viz.,
the fiber and free space. For the erbium-doped fiber, the small-signal gain is
$2.75$, and $E_{\mathrm{sat}}=63$ pJ for the single-soliton state (the value
is fitted by matching the bandwidth); for the SMs,
$E_{\mathrm{sat}}=80\leavevmode\nobreak\ $pJ is set (this value is fitted by
matching the TS in the SMs); and the internal linear transmission $\delta$ is
set to be $0.10$. For the homemade carbon-nanotube saturable absorber, we
measured its absorption curve, which is shown in Fig. S1(a). By fitting the
data with the help of the absorber equation, we obtained the following
parameters: $\alpha_{0}=0.095$, $\alpha_{\mathrm{ns}}=0.4908$, and
$I_{\mathrm{sat}}=21.6412$ MW/cm2.

For the $4f$ pulse shaper, which is modeled by Eq. (A2), we set
$\mathrm{Re}(\tilde{n})=10$, and vary $\mathrm{Im}(\tilde{n})$ from $6$ to
$30$, which determines $\tilde{\beta}_{k}$, as defined above. Because there
maybe a small discrepancy between the actual and theoretical parameters of the
holograms, we directly use the hologram’s phase to add the complex dispersion,
which may include higher-order terms (See Section 2 in Supplementary Material
for the discussion of the small discrepancy).

The production rate of the output coupler of the fiber cavity is set to be
$10\%$, and the overall internal loss of the $4f$ pulse shaper is about
$90\%$, due to the use the phase-amplitude encoding technique [44, 45]. With
these parameters, we have performed the simulations of the system for $2000$
round trips. In the simulations, we first generated a conventional soliton
$A_{\mathrm{soliton}}(T)$ from the random noise $\mathrm{Rand(N)}$ that which
is followed by the Gaussian distribution
[$=\mathrm{Rand(N)}\times\mathrm{\exp}(-(T/2)^{2})$], keeping the bandwidth of
the simulations as close as possible to the experimental value. Then, we used
the individual pulses to build two-soliton pairs
[${A_{\mathrm{soliton}}}(T-\tau_{0}/2)+{A_{\mathrm{soliton}}}(T+\tau_{0}/2)$,
$\tau_{0}$ being their relative separation] as the input for the next stage of
the simulations, with $\tilde{n}=10+6i$; in this case, the output produces a
stable bound state. This output is then used as the input for the next
iteration. Fig. S1(b) shows the output fields for
$\tilde{n}=10\left(1+i\right)$, between $1$ and $2000$ round trips. For
testing the stability, we calculated the difference in the pulse energy
between the given round trip and the last one, which is displayed in the
bottom portion of Fig. S1(b). The top inset in the same figure shows the
results for the round trips from $1$ to $200$. Monitoring the difference in
the pulse energy, we have concluded that the bound state remains stable after
an initial stage of strong oscillation. We extract the last output field after
$2000$ round trips and calculate their spectral intensity. All normalized
spectral intensities of the final outputs are separately shown in Fig. S1(c).
We conclude that the spectral fringes decrease with the increase of the
imaginary part of the GVD coefficient, $\mathrm{Im}(\tilde{n})$, which implies
that TS of the SMs grows accordingly. The Fourier transform of these final
output spectra produces the autocorrelation function which is displayed in
Fig. 2(g) of the main text.

The simulations were performed to verify the physical mechanism. The results
produced by the current theoretical model do not agree with experiments
exactly, as some simplifications are used in it.  For example, the difference
in separations between the experiments and simulations are due to a minor
difference between them in the GVD and dispersion loss. The difference in the
observed relative phase of the bound solitons is chiefly affected by the pump
power and polarization [6, 15]. In addition, in the simulation, we mainly
focused on the effect of the GVD and frequency-dependent loss on the TS in the
bound states, and we did not consider in detail the effects of the
polarization. Actually, the SLM can operate only in horizontal polarization.
Thus, in the stable dynamical regime, only the horizontal polarization is
effectively modulated, while the vertical polarization component is strongly
attenuated by losses. Nevertheless, the main features observed in the
experiments are correctly predicted by the present model, in which the values
of the TS in the SMs increase with the growth of the imaginary part of the
GVD, $\sim\mathrm{Im}(n)$.

###  5.2 The pulse shaper system and hologram

For the generation of a target arbitrary electric component of the optical
field, $E_{\mathrm{out}}(x,y)=A(x,y)\exp(i\phi)$, with a certain amplitude and
phase, the hologram needs to include the information concerning the amplitude
and phase. The hologram could be synthesized by the amplitude-phase encoding,
defined as [81, 44, 45]:

$\Phi(x,y)_{\mathrm{holo}}=\mathrm{sinc}(\pi-M\pi)\cdot\mathrm{Arg}(E)\cdot\Lambda(x,y),$
(A6)

where $M$ $(M=|A(x,y)|,0<M<1$) is a normalized constraint for the positive
output amplitude function, and $\mathrm{Arg}(E)=\exp(i(\phi))$ is the phase of
the output complex field. Further, $\Lambda(x,y)$ is an optical blazed grating
phase, that may be $\exp\left[{2\pi
i\left({{G_{x}}x+{G_{y}}y}\right)}\right]$, where $G_{x,y}$ are spatial
frequencies along the $x$ and $y$ directions. In our experiment, we used only
the optical grating oriented along $y$. The grating is necessary to avoid the
interference with effects of the zero-order diffraction in the blazed grating
phase.

To realize the modulation of the amplitude and phase in our study, i.e.,
$H(x)=\phi_{R}(\mathrm{Re}(\tilde{\beta}_{2,\mathrm{holo}}),x)\cdot\phi_{I}(\mathrm{Im}(\tilde{\beta}_{2,\mathrm{holo}}),x)\cdot\phi_{G}(d,x)$,
the amplitude (determined by $M$ in Eq. A6) should be taken as
$\phi_{I}(\mathrm{Im}(\tilde{\beta}_{2,\mathrm{holo}}),x)\sim\
\exp\left({-{k_{I}\mathrm{Im}(\tilde{n}){L_{\mathrm{cavity}}}\omega^{2}}/2}\right)$,
the phase information ($\mathrm{Arg}(E)$ in Eq. A6) is given by
$\phi_{R}(\mathrm{Re}(\tilde{\beta}_{2,\mathrm{holo}}),x)\sim\exp\left[{i{k_{R}\mathrm{Re}(\tilde{n})}{L_{\mathrm{cavity}}}{\omega^{2}}/2}\right]$,
and the optical blazing ($\Lambda(x,y)$) is $\phi_{G}(d,x)=\exp(-i2\pi Gy)$.
Here, $k_{R,I}$ are constant coefficients of the real and imaginary parts,
respectively. In actual encoding, we replace frequency $\omega$ by spatial
coordinate $x$, using the dispersion relationship of the $4f$ pulse shaper,
see Eq. (A5). Note that the current encoding format leads to a small
difference in the beam’s profile in comparison to the target profile [44], due
to higher-order effects. Details of the $4f$ pulse shaper, the hologram
structure and high-order effects are presented in Section 1 and 2 in
Supplementary Material.

Funding This work was supported by the International Postdoctoral Exchange
Fellowship Program 2020 by the Office of China Postdoctoral Council (No. 34
Document of OCPC,2020), the National Natural Science Foundation of China under
Grant Agreements 61705193, Fundamental Research Funds for the Central
Universities 2019QNA5003, Natural Science Foundation of Zhejiang Province
under Grant Nos. LGG20F050002. And in part, by the Israel Science Foundation
through grant No. 1286/17. Also, Shilong Liu and Ebrahim Karimi acknowledge
the support of Canada Research Chairs (CRC) and Canada First Research
Excellence Fund (CFREF) Program.

Acknowledgments We thank Prof. Xueming Liu (Zhejiang University) for providing
the part of the experimental platform, Prof. Baosen Shi and Dr. Zhiyuan Zhou
(University of Science and Technology of China) for providing the infrared
spatial light modulator and BBO crystal, Lin Huang (Zhejiang University) for
discussion of the results, and Ruimin Jie also Mengqin Gao (Zhejiang
University) for helping to measure parts of data in experiments. Also,
appreciation for Ende Zuo (University of Ottawa) for improvement of figure
esthetics.

Disclosures The auhors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

Supplemental document See Supplement 1 for supporting content.

##  6 References

## References

  * [1] F. M. Mitschke and L. F. Mollenauer, “Experimental observation of interaction forces between solitons in optical fibers,” Opt. Lett. 12, 355–357 (1987).
  * [2] B. A. Malomed, “Bound solitons in the nonlinear schrödinger/ginzburg-landau equation,” in _Large Scale Structures in Nonlinear Physics,_ (Springer, 1991), pp. 288–294.
  * [3] N. Akhmediev, A. Ankiewicz, and J. Soto-Crespo, “Multisoliton solutions of the complex ginzburg-landau equation,” Phys. Rev. Lett. 79, 4047 (1997).
  * [4] D. Tang, W. Man, H. Y. Tam, and P. Drummond, “Observation of bound states of solitons in a passively mode-locked fiber laser,” Phys. Rev. A 64, 033814 (2001).
  * [5] M. Grapinet and P. Grelu, “Vibrating soliton pairs in a mode-locked laser cavity,” Opt. Lett. 31, 2115–2117 (2006).
  * [6] G. Herink, F. Kurtz, B. Jalali, D. R. Solli, and C. Ropers, “Real-time spectral interferometry probes the internal dynamics of femtosecond soliton molecules,” Science 356, 50–54 (2017).
  * [7] 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).
  * [8] Z. 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, 1–11 (2019).
  * [9] N. K. Berger, B. Levit, V. Smulakovsky, and B. Fischer, “Complete characterization of optical pulses by real-time spectral interferometry,” Appl. Opt. 44, 7862–7866 (2005).
  * [10] A. Mahjoubfar, D. V. Churkin, S. Barland, N. Broderick, S. K. Turitsyn, and B. Jalali, “Time stretch and its applications,” Nat. Photonics 11, 341 (2017).
  * [11] J. Peng, M. Sorokina, S. Sugavanam, N. Tarasov, D. V. Churkin, S. K. Turitsyn, and H. Zeng, “Real-time observation of dissipative soliton formation in nonlinear polarization rotation mode-locked fibre lasers,” Communications Physics 1, 1–8 (2018).
  * [12] H. Liang, X. Zhao, B. Liu, J. Yu, Y. Liu, R. He, J. He, H. Li, and Z. Wang, “Real-time dynamics of soliton collision in a bound-state soliton fiber laser,” Nanophotonics 9, 1921–1929 (2020).
  * [13] M. Chernysheva, S. Sugavanam, and S. Turitsyn, “Real-time observation of the optical sagnac effect in ultrafast bidirectional fibre lasers,” APL Photonics 5, 016104 (2020).
  * [14] B. A. Malomed, “Bound states of envelope solitons,” Physical Review E 47, 2874 (1993).
  * [15] L. Gui, P. Wang, Y. Ding, K. Zhao, C. Bao, X. Xiao, and C. Yang, “Soliton molecules and multisoliton states in ultrafast fibre lasers: Intrinsic complexes in dissipative systems,” Applied Sciences 8, 201 (2018).
  * [16] W. He, M. Pang, D.-H. Yeh, J. Huang, C. Menyuk, and P. S. J. Russell, “Formation of optical supramolecular structures in a fibre laser by tailoring long-range soliton interactions,” Nat. Commun. 10, 1–9 (2019).
  * [17] X. Liu, X. Yao, and Y. Cui, “Real-time observation of the buildup of soliton molecules,” Phys. Rev. Lett. 121, 023905 (2018).
  * [18] H. Sakaguchi, D. V. Skryabin, and B. A. Malomed, “Stationary and oscillatory bound states of dissipative solitons created by third-order dispersion,” Opt. Lett. 43, 2688–2691 (2018).
  * [19] M. Stratmann, T. Pagel, and F. Mitschke, “Experimental observation of temporal soliton molecules,” Phys. Rev. Lett. 95, 143902 (2005).
  * [20] F. Kurtz, C. Ropers, and G. Herink, “Resonant excitation and all-optical switching of femtosecond soliton molecules,” Nat. Photonics 14, 9–13 (2020).
  * [21] M. Pang, W. He, X. Jiang, and P. S. J. Russell, “All-optical bit storage in a fibre laser by optomechanically bound states of solitons,” Nat. Photonics 10, 454–458 (2016).
  * [22] J. Girardot, A. Coillet, M. Nafa, F. Billard, E. Hertz, and P. Grelu, “On-demand generation of soliton molecules through evolutionary algorithm optimization,” Optics letters 47, 134–137 (2022).
  * [23] P. Grelu, F. Belhache, F. Gutty, and J.-M. Soto-Crespo, “Phase-locked soliton pairs in a stretched-pulse fiber laser,” Opt. Lett. 27, 966–968 (2002).
  * [24] T. Zhu, Z. Wang, D. Wang, F. Yang, and L. Li, “Observation of controllable tightly and loosely bound solitons with an all-fiber saturable absorber,” Photonics Res. 7, 61–68 (2019).
  * [25] Y. S. Kivshar and B. A. Malomed, “Dynamics of solitons in nearly integrable systems,” Rev. Mod. Phys. 61, 763 (1989).
  * [26] F. Leo, S. Coen, P. Kockaert, S.-P. Gorza, P. Emplit, and M. Haelterman, “Temporal cavity solitons in one-dimensional kerr media as bits in an all-optical buffer,” Nat. Photonics 4, 471–476 (2010).
  * [27] J. K. Jang, M. Erkintalo, S. Coen, and S. G. Murdoch, “Temporal tweezing of light through the trapping and manipulation of temporal cavity solitons,” Nat. Commun. 6, 1–7 (2015).
  * [28] W. He, M. Pang, D.-H. Yeh, J. Huang, and P. S. J. Russell, “Synthesis and dissociation of soliton molecules in parallel optical-soliton reactors,” Light Sci. Appl. 10, 1–15 (2021).
  * [29] B. A. Malomed, _Soliton management in periodic systems_ (Springer Science & Business Media, 2006).
  * [30] S. K. Turitsyn, B. G. Bale, and M. P. Fedoruk, “Dispersion-managed solitons in fibre systems and lasers,” Phys. Rep. 521, 135–203 (2012).
  * [31] R. I. Woodward, “Dispersion engineering of mode-locked fibre lasers,” J. Opt. 20, 033002 (2018).
  * [32] A. F. Runge, D. D. Hudson, K. K. Tam, C. M. de Sterke, and A. Blanco-Redondo, “The pure-quartic soliton laser,” Nat. Photonics 14, 492–497 (2020).
  * [33] R. Iegorov, T. Teamir, G. Makey, and F. Ilday, “Direct control of mode-locking states of a fiber laser,” Optica 3, 1312–1315 (2016).
  * [34] X. Wei, J. C. Jing, Y. Shen, and L. V. Wang, “Harnessing a multi-dimensional fibre laser using genetic wavefront shaping,” Light Sci. Appl. 9, 1–10 (2020).
  * [35] A. Komarov, H. Leblond, and F. Sanchez, “Multistability and hysteresis phenomena in passively mode-locked fiber lasers,” Phys. Rev. A 71, 053809 (2005).
  * [36] A. Monmayrant, S. J. Weber, and B. Chatel, “A newcomer’s guide to ultrashort pulse shaping and characterization,” J. Phys. B 43, 103001 (2010).
  * [37] A. M. Weiner, “Ultrafast optical pulse shaping: A tutorial review,” Opt. Commun. 284, 3669–3692 (2011).
  * [38] Y. Song, X. Shi, C. Wu, D. Tang, and H. Zhang, “Recent progress of study on optical solitons in fiber lasers,” Applied Physics Reviews 6, 021313 (2019).
  * [39] D. V. Skryabin and A. V. Gorbach, “Colloquium: Looking at a soliton through the prism of optical supercontinuum,” Rev. Mod. Phys. 82, 1287 (2010).
  * [40] N. Akhmediev, A. Ankiewicz, and J. Soto-Crespo, “Stable soliton pairs in optical transmission lines and fiber lasers,” JOSA B 15, 515–523 (1998).
  * [41] G. P. Agrawal, “Nonlinear fiber optics,” in _Nonlinear Science at the Dawn of the 21st Century,_ (Springer, 2000), pp. 195–211.
  * [42] Y. S. Kivshar and G. P. Agrawal, _Optical solitons: from fibers to photonic crystals_ (Academic press, 2003).
  * [43] P. Grelu and N. Akhmediev, “Dissipative solitons for mode-locked lasers,” Nat. Photonics 6, 84–92 (2012).
  * [44] E. Bolduc, N. Bent, E. Santamato, E. Karimi, and R. W. Boyd, “Exact solution to simultaneous intensity and phase encryption with a single phase-only hologram,” Opt. Lett. 38, 3546–3549 (2013).
  * [45] S.-L. Liu, Q. Zhou, S.-K. Liu, Y. Li, Y.-H. Li, Z.-Y. Zhou, G.-C. Guo, and B.-S. Shi, “Classical analogy of a cat state using vortex light,” Commun. Phys 2, 1–9 (2019).
  * [46] B. A. Malomed, “Propagation of a soliton in a nonlinear waveguide with dissipation and pumping,” Opt. Commun. 61, 192–194 (1987).
  * [47] S. Kelly, “Characteristic sideband instability of periodically amplified average soliton,” Electron. Lett. 28, 806–807 (1992).
  * [48] E. Ott, _Chaos in dynamical systems_ (Cambridge university press, 2002).
  * [49] Y. Cui, Y. Zhang, Y. Song, L. Huang, L. Tong, J. Qiu, and X. Liu, “Xpm-induced vector asymmetrical soliton with spectral period doubling in mode-locked fiber laser,” Laser & Photonics Reviews 15, 2000216 (2021).
  * [50] V. Torres-Company, J. Lancis, and P. Andres, “Space-time analogies in optics,” Progress in Optics 56, 1–80 (2011).
  * [51] V. V. Afanasjev, B. A. Malomed, and P. Chu, “Stability of bound states of pulses in the ginzburg-landau equations,” Physical Review E 56, 6020 (1997).
  * [52] X. Liu, “Resolving the temporal dynamics of mode-locked laser with single-shot time-microscope,” arXiv preprint arXiv:2111.00957 (2021).
  * [53] T. Xian, L. Zhan, W. Wang, and W. Zhang, “Subharmonic entrainment breather solitons in ultrafast lasers,” Physical Review Letters 125, 163901 (2020).
  * [54] X. Liu, X. Han, and X. Yao, “Discrete bisoliton fiber laser,” Sci. Rep. 6, 1–8 (2016).
  * [55] J. Peng, S. Boscolo, Z. Zhao, and H. Zeng, “Breathing dissipative solitons in mode-locked fiber lasers,” Sci. Adv. 5, eaax1110 (2019).
  * [56] M. Liu, Z.-W. Wei, H. Li, T.-J. Li, A.-P. Luo, W.-C. Xu, and Z.-C. Luo, “Visualizing the “invisible” soliton pulsation in an ultrafast laser,” Laser Photonics Rev. 14, 1900317 (2020).
  * [57] Y. Wang, C. Wang, F. Zhang, J. Guo, C. Ma, W. Huang, Y. Song, Y. Ge, J. Liu, and H. Zhang, “Recent advances in real-time spectrum measurement of soliton dynamics by dispersive fourier transformation,” Reports on Progress in Physics (2020).
  * [58] D. Tang, B. Zhao, L. Zhao, and H. Y. Tam, “Soliton interaction in a fiber ring laser,” Physical Review E 72, 016616 (2005).
  * [59] J. T. Barreiro, T.-C. Wei, and P. G. Kwiat, “Beating the channel capacity limit for linear photonic superdense coding,” Nat. Phys. 4, 282–286 (2008).
  * [60] B. Corcoran, M. Tan, X. Xu, A. Boes, J. Wu, T. G. Nguyen, S. T. Chu, B. E. Little, R. Morandotti, A. Mitchell, and D. J. Mosskomarov2012multiple, “Ultra-dense optical data transmission over standard fibre with a single chip source,” Nat. Commun. 11, 1–7 (2020).
  * [61] Y. Song, F. Zhou, H. Tian, and M. Hu, “Attosecond timing jitter within a temporal soliton molecule,” Optica 7, 1531–1534 (2020).
  * [62] U. Keller, “Recent developments in compact ultrafast lasers,” Nature 424, 831–838 (2003).
  * [63] R. W. Boyd, _Nonlinear optics_ (Academic press, 2020).
  * [64] Z. Sun, D. Popa, T. Hasan, F. Torrisi, F. Wang, E. J. Kelleher, J. C. Travers, V. Nicolosi, and A. C. Ferrari, “A stable, wideband tunable, near transform-limited, graphene-mode-locked, ultrafast laser,” Nano Research 3, 653–660 (2010).
  * [65] Y. Xiang, Y. Luo, B. Liu, Z. Yan, Q. Sun, and D. Liu, “Observation of wavelength tuning and bound states in fiber lasers,” Sci. Rep. 8, 1–9 (2018).
  * [66] P. M. Becker, A. A. Olsson, and J. R. Simpson, _Erbium-doped fiber amplifiers: fundamentals and technology_ (Elsevier, 1999).
  * [67] S. Fedorov, N. Veretenov, and N. Rosanov, “Irreversible hysteresis of internal structure of tangle dissipative optical solitons,” Phys. Rev. Lett. 122, 023903 (2019).
  * [68] S. Smith, “Optical bistability: towards the optical computer,” Nature 307, 315–316 (1984).
  * [69] S. Smith, “Optical bistability, photonic logic, and optical computation,” Appl. Opt. 25, 1550–1564 (1986).
  * [70] X. Yi, Q.-F. Yang, X. Zhang, K. Y. Yang, X. Li, and K. Vahala, “Single-mode dispersive waves and soliton microcomb dynamics,” Nat. Commun. 8, 1–9 (2017).
  * [71] A. Komarov, F. Amrani, A. Dmitriev, K. Komarov, D. Meshcheriakov, and F. Sanchez, “Multiple-pulse operation and bound states of solitons in passive mode-locked fiber lasers,” International Journal of Optics 2012 (2012).
  * [72] X. Fan, S. Wang, Y. Wang, D. Shen, D. Tang, and L. Zhao, “Pump hysteresis and bistability of dissipative solitons in all-normal-dispersion fiber lasers,” Appl. Opt. 54, 3774–3779 (2015).
  * [73] J. Feldmann, N. Youngblood, M. Karpov, H. Gehring, X. Li, M. Stappers, M. Le Gallo, X. Fu, A. Lukashchuk, A. S. Raja _et al._ , “Parallel convolutional processing using an integrated photonic tensor core,” Nature 589, 52–58 (2021).
  * [74] T. Zhou, X. Lin, J. Wu, Y. Chen, H. Xie, Y. Li, J. Fan, H. Wu, L. Fang, and Q. Dai, “Large-scale neuromorphic optoelectronic computing with a reconfigurable diffractive processing unit,” Nat. Photonics 15, 367–373 (2021).
  * [75] G. Genty, L. Salmela, J. M. Dudley, D. Brunner, A. Kokhanovskiy, S. Kobtsev, and S. K. Turitsyn, “Machine learning and applications in ultrafast photonics,” Nat. Photonics 15, 91–101 (2021).
  * [76] G. Pu, L. Yi, L. Zhang, C. Luo, Z. Li, and W. Hu, “Intelligent control of mode-locked femtosecond pulses by time-stretch-assisted real-time spectral analysis,” Light Sci. Appl. 9, 1–8 (2020).
  * [77] P. Drummond, R. Shelby, S. Friberg, and Y. Yamamoto, “Quantum solitons in optical fibres,” Nature 365, 307–313 (1993).
  * [78] M. Kues, C. Reimer, P. Roztocki, L. R. Cortés, S. Sciara, B. Wetzel, Y. Zhang, A. Cino, S. T. Chu, B. E. Little _et al._ , “On-chip generation of high-dimensional entangled quantum states and their coherent control,” Nature 546, 622–626 (2017).
  * [79] S. Scholes, R. Kara, J. Pinnell, V. Rodríguez-Fajardo, and A. Forbes, “Structured light with digital micromirror devices: a guide to best practice,” Opt. Eng. 59, 041202 (2019).
  * [80] X. Liu, D. Han, Z. Sun, C. Zeng, H. Lu, D. Mao, Y. Cui, and F. Wang, “Versatile multi-wavelength ultrafast fiber laser mode-locked by carbon nanotubes,” Sci. Rep. 3, 1–6 (2013).
  * [81] J. A. Davis, D. M. Cottrell, J. Campos, M. J. Yzuel, and I. Moreno, “Encoding amplitude information onto phase-only filters,” Appl. Opt. 38, 5004–5013 (1999).