On-demand harnessing of photonic soliton molecules

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.


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 timestretch 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,8,[14][15][16]. Much interest has been drawn to the TS degree of freedom due to the fact that it affects numerous phenomena, including nonlinear evolution [6,7,17], short-and long-range interactions [3,16], SM complexes [8], and the Hopf bifurcation which excites robust intrinsic oscillations in SMs [5,18]. 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,15,23], and TS corresponding to tightly or loosely bound "molecules" [4,6,16,23,24]. Such SMs were predicted by theoretical models based on the nonlinear Schrödinger/Ginzburg-Landau equations [2,3,25]. 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" [6,7,15,23]. 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 [21,26], and soliton trapping [27]. Recently, the technique of active pump modulation was adopted to efficiently manipulate SMs in modelocked 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 groupvelocity 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,18,35] 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 4 f 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. 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 G 1,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 × 1080 and the pixel pitch of 6.4 µm.  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 4 f pulse shaper to morph the amplitude and phase of the spectrum. The total length of the cavity is ≈ 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 ps 2 , i.e., the mode-locked laser operates in the anomalous GVD region.

METHODS FOR ON-DEMAND GENERATION OF SMS
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 (z, t) of the output optical pulse is [32,39] ∂A ∂z Here, z and t are the propagation distance and local time; γ is the Kerr parameter; g and α represent the linear gain and loss, respectively; and H.O.E stands for terms representing higherorder 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 [3,18]. In this model,β k = Re(β k ) + i Im(β k ) includes GVD, Re(β k ), and the dispersion losses, Im(β k ), of orders k = 2, 3, 4, .... In this paper,β k are referred to as k-th order complex dispersion coefficients. In particular,β 2 represents the complex second-order dispersion. Its GVD part, Re(β 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, Im(β 2 ), being a necessary ingredient for the creation of two-and multi-soliton bound states in dissipative nonlinear systems [2,3,18,35,40], 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 Ω g , viz., Im(β 2 ) = g/2Ω 2 g [38]. In our system, β k can be precisely adjusted by means of the 4 f 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, Re β 2 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 [32,41,42], while the loss part, Im β 2 , 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β 2,fiber as the coefficient in the fiber section, andβ 2,holo as one imprinted by the hologram in the 4 f 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β 2,holo . In this context, Fig. 1(b) shows that TS between the bound solitons gradually increases, following the increase of Im(β 2,holo ) in a certain range. In particular, only the conventional single soliton, symbolized by |0 , is possible at Im(β 2,holo ) = 0. With the increase of Im(β 2,holo ), the system can create SMs ondemand with different values of TS, |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 |4 A,B , where A and B represent the central wavelength of two SMs. Experimentally, both the GVD and dispersion-loss constituents ofβ 2,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) = φ R (Re(β 2,holo ), x) · φ I (Im(β 2,holo ), x) · φ G (d, x), which includes three factors used for the implementation of the precise manipulations. Here, φ R (Re(β 2,holo ), x) and φ I (Im(β 2,holo ), x) represent real and imaginary terms (corresponding to the phase and amplitude encoding, respectively), while φ 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ñ to encode the hologram by means of its real and imaginary parts: Im β 2,holo = k I · Im (ñ) .

(2b)
Here, k R = −0.2170 × 10 −3 ps 2 /m and k I = 0.9309 × 10 −3 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). Re (ñ) represents the positive or negative number associated with GVD, while Im (ñ) 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 where L cavity is the total length of the cavity (as mentioned above, it is 19.88 m in our experimental setup), ω represents the optical frequency corresponding to the wavelength at the spatial position of the SLM, and δω 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 δω = 0, whereñ = 10 (1 + i) 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 δω corresponding to the wavelength shift δλ = 2.5 nm.
In the present regime, the real and imaginary parts ofβ 2,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.

A. Manipulations of the TS (temporal separation) for SMs (soliton molecules) in the mode-locked laser
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ñ = 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 [17,46,47]. 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,ñ, 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 [6,9]. We set Re(ñ) = 10 and increase Im(ñ), 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 Im(ñ) = 6. Later, the output features SMs, whose TS increases linearly up to Im(ñ) ≈ 29. Fig. 2(b) shows a power spectrum representing the SM forñ = 10 (1 + i), 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 Im(ñ), viz., 0 and 6 ≤ Im(ñ) ≤ 30 with interval ∆ (Im(ñ)) = 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 [6,9]. Figure 2(d) exhibits the Fourier transform of the corresponding autocorrelation function. The results show that the TS grows linearly with an increase of Im(ñ) from 6 to 28. With the further increase of Im(ñ), the SM becomes unstable, switching to another SM state with a large TS, as observed in Fig. 2(d) for Im(ñ) = 30. This phenomenology is similar to that demonstrating the Hopf-type bifurcation in two-soliton bound states in other systems [5,18,48]. 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 singleshot 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 Fouriertransform spectra are shown in Fig. 2(e), where we present the conventional single soliton (at Im(ñ) = 0) and three SMs (at Im(ñ) = 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× Im(ñ) + 2.342, 6 ≤ Im(ñ) ≤ 28. Close to Im(ñ) = 29, a large jump in the TS appears, and then it does not change conspicuously with the further increase of Im(ñ). Here, only the range of Im(ñ) 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 [3,6,51], 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 = 1 Q ave |A(m, ω)| 2 dω is the spectrum energy for m-th round trip from TSDFT, and 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 Im(ñ) =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 ±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 τ m − τ m−1 . We add the separation difference at the below of autocorrelation functions. Here, two green basic lines set ±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 ±100 fs. These results illustrate the SMs (Im(ñ) = 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ñ = 10(1 + i). It shows the average signal-noise level could be up to around 50 dB.
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ñ = 10 + 6i, which leads to the creation of bound soliton pairs. A similar simulation, forñ = 10 (1 + i), 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 Im(ñ) in the range of Im(ñ) 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 The other is the values of TS: the predicted and experimentally observed separations are not exactly equal for each Im(ñ). 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 4 f 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 4 f 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(ñ), viz., Im(ñ) > 35. This fact may quantitatively explain the jump that happens at larger Im(ñ) . 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.

B. 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 |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 Im (ñ) 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 Im (ñ) = 10, 14, 18, and 22, as states {|0 , |1 , |2 , |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 (Im(ñ) = 0) and SMs with a very large separation (Im(ñ) ≥ 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 Im (ñ) between 6 and 28.
Thus, we start by using the four above-mentioned SMs, |m with m = 0, 1, 2, 3, and construct a simple test string {0|3|1|2|0|3}. The respective set of six holograms are consecutively displayed on the SLM, switching them automatically in sequence with the time interval of δ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.
We then consider a more complex encoding, viz., the ψ-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, {0000|0000|0320|0020|1333|0020|0320|0000|0000}. 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, {0|3|1|2|0|3}. 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 |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, {0000|0000|0320|0020|1333|0020|0320|0000|0000}, being exactly the same as the encoding sequence defining the ψ-shaped target. The dynamic switching results for the ψ-shaped patterns in the quaternary-encoding regime is presented in Supporting Material for Video 2.
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.

C. 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 highprecious tuning their wavelengths. The present setup may offer In these results, δx = 1 means that the actual spatial shift of the hologram in SLM is 10 µm.
a simple platform to adjust the wavelength of SMs. According to Eq. (3), δω = 2πc/λ 2 0 δλ represents the shift with respect to the central wavelength λ 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 δx, which affects the real and imaginary parts of GVD in SLM through the space-frequency relation of the 4 f shaping system [36], the effective gain area displaced by δω 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ñ = 10 + 30i, for which the central wavelength is 1554.6 nm. Then, after moving the position of the hologram in SLM by δ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 ±0.288 ps. Figure 4(b) shows the relation between δ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, λ 0 (µm) = 0.0455 × δx + 1.5546. Thus, our results exhibit the wavelength shift by 4 nm (from 1555 nm to 1559 nm), with the linear dependence on δx. However, the linear relationship cannot be valid in the entire range of the variation of δ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 δx. However, the original TS is not exactly recovered when δ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ñ = 10 + 19i, as functions of δx. The top plot shows the situation corresponding to varying δx from −10 to +10, while TS τ jumps from 4.7 to 3 ps around δx = 2. The inverse jump happens around δx = −2, when δ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 δ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 δx. Such bistability areas are key elements used in optical logic devices [68,69].
Geneally, hysteresis is an interesting phenomenon in modelocked fiber laser systems. For example, the hysteresis can be induced by variations of the power or polarization [35,67], 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 |m A with wavelength λ A , and as |m B with another wavelength λ 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 δx of the hologram in the free space may be appropriate for the design of optical free-space computational processing [74].

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 [21,26]. 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 [33,34,75,76], and operations with high dimensional temporal quantum information [77,78].
The results demonstrate that the hologram-driven modelocking is a robust and accurate technique for the work with the optical SMs. The switching speed of the present setup, ∼ 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 ∼ 100 MHz range [37]. In the latter case, it is possible to monitor the dynamics of switching between different SMs.

A. 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 carbonnanotube-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 singlemode one and the doped gain fiber), which provide GVD and gain. Accordingly, the model includes CGLEs and the additional ingredient, defined as follows: ∂A ∂z = −iβ 2,fiber 2 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 pulseshaping equation. Therefore, the free-space element (the 4 f pulse shaper [36]) that connects to the complex dispersion coefficient, is modeled by: (A2) Here, L cavity is the overall length of the cavity, L 4f is the overall length of the 4 f pulse shaper, δ 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β k,holo represents, as said above, the additional k-order complex dispersion. Recall thatβ 2,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. [17,80]. It is characterized by the transfer function, Here α ns is the unsaturated absorption, α 0 is the linear limit of the saturable absorption, I sat is the saturation intensity, and the effective fiber's cross-section area is A eff = 63.6 × 10 −12 m 2 . 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 4 f 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 → A · T 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 sat : where E(z) = +∞ −∞ |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 4 f pulse shaper, along with the surface of the SLM, follows this simple relationship based on the diffraction equations [36,37]: Here, m is the diffraction order (usually m = 1), θ diff is the input angle of the optical grating (here, it is 47.5 • ), f = 125 mm is the focal length of the input lens, and d 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 ∂λ/∂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 Ω g in Eq. (A1). 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 Re(β 2,gain−fiber ) = 20 ps 2 /km and γ = 0.0032, the length of the single-mode fiber is 15.48 m, with Re β 2,single−mode−fiber = −21.6 ps 2 /km and γ = 0.0013, and the net dispersion of the entire fiber is −0.2672 ps 2 . 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 singlemode fiber is placed in front of the 4 f 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 sat = 63 pJ for the single-soliton state (the value is fitted by matching the bandwidth); for the SMs, E sat = 80 pJ is set (this value is fitted by matching the TS in the SMs); and the internal linear transmission δ 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: α 0 = 0.095, α ns = 0.4908, and I sat = 21.6412 MW/cm 2 .
For the 4 f pulse shaper, which is modeled by Eq. (A2), we set Re(ñ) = 10, and vary Im(ñ) from 6 to 30, which determines β 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 4 f 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 soliton (T) from the random noise Rand(N) that which is followed by the Gaussian distribution [= Rand(N) × 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 soliton (T − τ 0 /2) + A soliton (T + τ 0 /2), τ 0 being their relative separation] as the input for the next stage of the simulations, withñ = 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ñ = 10 (1 + i), 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 dis-played 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, Im(ñ), 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, ∼ Im(n).

B. The pulse shaper system and hologram
For the generation of a target arbitrary electric component of the optical field, E out (x, y) = A(x, y) exp(iφ), 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 [44,45,81]: Φ(x, y) holo = sinc(π − Mπ) · Arg(E) · Λ(x, y), where M (M = |A(x, y)|, 0 < M < 1) is a normalized constraint for the positive output amplitude function, and Arg(E) = exp(i(φ)) is the phase of the output complex field. Further, Λ(x, y) is an optical blazed grating phase, that may be exp 2πi G x x + G y y , 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) = φ R (Re(β 2,holo ), x) · φ I (Im(β 2,holo ), x) · φ G (d, x), the amplitude (determined by M in Eq. A6) should be taken as φ I (Im(β 2,holo ), x) ∼ exp −k I Im(ñ)L cavity ω 2 /2 , the phase information (Arg(E) in Eq. A6) is given by φ R (Re(β 2,holo ), x) ∼ exp ik R Re(ñ)L cavity ω 2 /2 , and the optical blazing (Λ(x, y)) is φ G (d, x) = exp(−i2πGy). Here, k R,I are constant coefficients of the real and imaginary parts, respectively. In actual encoding, we replace frequency ω by spatial coordinate x, using the dispersion relationship of the 4 f pulse shaper, see Eq. (A5). Note that the current encoding format leads to Fig. 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 α 0 = 0.095, α ns = 0.4908, I sat = 21.6412 W/m 2 . In panel (b), the SM's output intensity is displayed, forñ = 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, τ = 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 Im(ñ), and the vertical axis is the wavelength, the central one being 1553 nm. a small difference in the beam's profile in comparison to the target profile [44], due to higher-order effects. Details of the 4 f pulse shaper, the hologram structure and high-order effects are presented in Section 1 and 2 in Supplementary Material.