## 4 Numerical Simulation of Soliton Molecules  
The solutions of soliton molecules were first discovered mainly by theoretical analysis of nonlinear Schrodinger equation and complex Ginzburg-Landau equation, which are distributed propagation models and allow a good approximation of the dominating effects in passively mode-locked fibre lasers. Different approaches such as perturbation theory [11], inverse scattering method [57], etc. enable simplified theoretical framework to address the interaction mechanisms and the general behaviours of bound solitons, which were rather controversial in terms of the stability properties [14; 15]. Nevertheless, it is noteworthy that those theoretical researches opened exciting avenues for understanding of multipulse interactions in fibre lasers and proved possible capacity scaling in optical communication, inspiring a number of successive experimental findings. For details, see the pioneer work by [PERSON], [PERSON], [PERSON] etc. [11; 12; 13; 14; 58; 59] and references therein. Later, more and more numerical analysis was made based on lumped model which enables better comparison with experiments since real parametric values such as dispersion and nonlinear Kerr coefficient are involved. Up till now, a plethora of fibre laser cavities with distinct dispersion regimes [19; 44], mode-locking techniques [4; 60], etc. have been numerically investigated, giving rise to different features of bound states of multiple pulses, for instance, different envelopes and modulation fringes of the optical spectra. The general dynamics of soliton molecules and crucial influences of some cavity parameters have been unravelled in detail, providing fundamental insights into the underlying physics of multisoliton interactions. Some instructive work can be found in refs. [4; 5; 19; 44; 60; 61]. In this review we would like to demonstrate the basic method of carrying out the numerical simulation of two-soliton molecules. To provide a direct comparison with experiments, impact from several experimentally controllable variables including initial condition, gain and higher-order dispersion are shown respectively. Here we use a modified nonlinear Schrodinger equation based on lumped description of functions of different components. We neglect birefringence and assume a fast saturable absorber in order to simplify numerical calculation without losing the most important impact from cavity parameters. The pulse propagation in gain fibre can be described by a scalar formula  
\[i\frac{\partial A}{\partial z}=\frac{\beta_{2}}{2}\frac{\partial^{2}A}{ \partial T^{2}}+i\frac{\beta_{3}}{6}\frac{\partial^{3}A}{\partial T^{3}}- \gamma|A|^{2}A+i\frac{1}{2}\cdot\frac{g_{0}}{1+E_{pulse}/E_{sat}}(1+T_{2}^{2} \frac{\partial^{2}}{\partial T^{2}})A \tag{6}\]  
Figure 14: Autocorrelation traces of **(a)** a five-soliton bound state and **(b)** a twin-pulse bound state. (Reproduced with permission from [47]. © The Optical Society of America (Washington, DC, USA), 2007). where \(A\) denotes the pulse envelope and \(T\) is the time in a frame of reference moving with the group velocity \(v_{g}\). \(\beta_{2}\) and \(\beta_{3}\) are the second- and third-order dispersion, \(\gamma\) is the nonlinear Kerr coefficient and \(g_{0}\) is the small-signal gain. \(E_{pulse}\) is the pulse energy while \(E_{sat}\) is the saturation energy. Gain dispersion is taken into account in the manner \(T_{2}=2\pi/(ck^{2}\Delta\lambda_{g})\) with \(\Delta\lambda_{g}\) being gain bandwidth, in which \(c\) is the speed of light in vacuum and \(k\) is the wave vector. \(g_{0}\) is set to zero when the pulses propagate in passive fibres. The simulated laser cavity is primarily comprised of three segments of optical fibres schematically shown in Figure 15. The erbium-doped fibre provides small-signal gain \(g_{0}\) = 3.5 m\({}^{-1}\), with saturation energy \(E_{sat}\) = 0.5 nJ and gain bandwidth of 80 nm.