Passively mode-locked fiber lasers can produce ultrafast laser outputs in the picosecond or even femtosecond range, thus they have high research value and widespread application fields, such as optical communication, biomedicine, precision machining and so on. With the performance improvement of ultrafast mode-locked fiber lasers, their composition about the passive and active fiber and structures become more complex. To guide the experimental design of the mode-locked fiber lasers, their simulation models for simulating and analyzing the ultrashort pules propagation are more complex and their simulation round trip reach to 2000 times^[1] with long time-consuming and low computation efficiency^[2-3]. As we all know, computational efficiency decreases as the accuracy of the computation increases.

The nonlinear Schrödinger equation (NLSE) is the most basic equation for simulating the propagation of nonlinear pulses in the passive fibers^[4]. It could consider and describe the effects of second-order dispersion and self-phase modulation (SPM). There is also gain and gain dispersion in active fibers, which are usually described using the Ginzburg–Landau equation (GLE). The modeling of ultrashort pulse fiber laser mainly relies on the calculation of these two equations.

In the previous work, algorithms for solving these nonlinear equations have been developed a lot. Split-step Fourier method (SSFM) has been widely used for numerical solutions in NLSE and GLE^[5-6]. It is faster than most finite difference methods, mainly due to the use of the Fast Fourier Transform (FFT) algorithm^[7]. Then on this basis, a symmetric split-step Fourier method and a fourth-order Runge-Kutta method for the nonlinear part (SSFM-RK4) are derived^[8]. The fourth-order Runge-Kutta in the Interaction Picture (RK4IP) method was originally used for Bose-Einstein condensates, but proved to be effective and precise. Because it eliminates the error that comes with calculating the dispersive and nonlinear parts separately, the two parts occur simultaneously^[9-10]. It was found in some recent studies that the numerical solutions provided by SSFM and symmetric SSFM are more stable than the numerical solutions of the RK4IP method. Furthermore, the RK4IP method introduces an unexpected energy change, which cannot reproduce the fine structure of the time distribution in the proposed soliton dynamics^[11]. SSFM has also developed variants for different situations, e.g. the embedded split-step Fourier method (embedded-SSFM) is proposed to reproduce the soliton fission dynamics^[12].

Various kinds of methods related with the SSFM including the nonlinear phase-rotation method, the logarithmic step-size method, the walk-off method and the constant step-size method have been proposed to improve the precision or computational efficiency, of which the local error method^[13] is used widely. Usually, these methods are proposed for some specific systems. For example, an adaptive split-step Fourier method (ASSFM) is proposed according to the parameters of photonic crystal fibers (PCF), which due to the minimum area mismatch (MAM) further improve the computational efficiency^[14]. The conservation quantity error (CQE) method is proposed as an enhanced local error method, which can be utilized to combine with higher-order integral algorithms^[15]. There are also adaptive methods proposed to limit the step size according to specific judgment criteria^[16].

The most commonly-used algorithms for the simulation of mode-locked lasers are the analytical methods and the simulation model based SSFM^[17-18] which brings the long computational time and low computational efficiency. An adaptive method needs to be put forward for improving the computational efficiency of the SSFM method and decreasing the computation time in ultrafast mode-locked fiber laser.

In this paper, an efficient and competitive SSFM-GELE method is proposed to improve the computational efficiency of the mode-locked fiber laser model. The estimation of the global error as the criterion to adaptively adjust the step size of the fiber is obtained from the local energy increment.

For the establishment of the theoretical model of passive mode-locked fiber lasers, there are mainly two ways. One is to normalize various system parameters into one equation to solve its self-consistent solution of the equation, which is the pulse output of the mode-locked fiber laser^[19]. Its disadvantage is that the intracavity evolution characteristics of the pulse cannot be obtained.

For the other way, the actual operation of each device and fiber in the laser cavity is individually considered, and the pulse is gradually transmitted in these devices and fibers^[20], so that the intracavity evolution characteristics can be obtained. Furthermore, the pulse tracing method is used to simulate the operation of the pulse in the cavity. Different propagation equations are used for different devices in the cavity. Specifically, NLSE is used in Single Mode Fiber (DCF) and GLE is used in active fiber. For saturable absorber, coupler filter, etc., the transfer function of the devices is used directly. In the model, the initial pulse can be set to white noise. After the pulse runs for one circle in the cavity, the output result is used as the initial pulse input for the next circle, so that the mode-locked pulse output can be realized in this model. It is worth emphasizing that this method provides insight into the formation of optical pulses and their dynamics, especially the influence of different devices in the cavity on pulse shaping, which is usually difficult to observe by experiment. Here, this is chosen to carry out numerical simulation on ultrashort pulse fiber laser in this paper.

The equations describing the propagation of laser pulses in fibers are NLSE and GLE. The optical pulses transmitted in the SMF will be affected by the group velocity dispersion (GVD) and the nonlinear effect. To obtain the equation describing the evolution of the optical pulses in the SMF, we can start from the Maxwell's equations, through a series of approximations, and finally obtain the NLSE. The NLSE can be written as

$ \frac{\partial {{\boldsymbol{A}}}}{\partial \textit{z}}=-i\frac{{\beta }_{2}}{2}\frac{{\partial }^{2}{{\boldsymbol{A}}}}{\partial {T}^{2}}+i\gamma {\left|{{\boldsymbol{A}}}\right|}^{2}{{\boldsymbol{A}}}\quad. $ (1)

In formula (1), A is the amplitude of the pulse envelope, z is the distance of the pulse transmission, and $ \;{\beta }_{2} $ is the GVD coefficient. The quantity $T= t-{\textit{z}}/v_g\equiv t-\beta_1{\textit{z}}$, is the retarded time traveling at the envelope group velocity, where $ {v}_{g} $ is the group velocity and β₁ is the reciprocal of v_g. $ \gamma $ is nonlinear coefficient. This equation is the simplest nonlinear equation form to study pulse evolution in optical fibers. According to different experimental conditions, this equation also could be modified.

The active fiber generally refers to the fiber doped with rare earth ions, and its gain effect has an important influence on the evolution of the optical pulse. In addition to the nonlinear effect and the influence of GVD in the active fiber, there are also gain and gain dispersion, as cannot be described by NLSE. Usually, the GLE is used to describe the pulse propagation in the active fiber. The GLE can be written as

$ \begin{split} &\frac{\partial {{\boldsymbol{A}}}}{\partial {\textit{z}}}+\frac{\alpha }{2}{{\boldsymbol{A}}}+i\frac{{\beta }_{2}}{2}\frac{{\partial }^{2}{{\boldsymbol{A}}}}{\partial {T}^{2}}-\frac{{\beta }_{3}}{6}\frac{{\partial }^{3}{{\boldsymbol{A}}}}{\partial {T}^{3}}=\\ &i\gamma {\left|{{\boldsymbol{A}}}\right|}^{2}{{\boldsymbol{A}}}+\frac{g}{2}{{\boldsymbol{A}}}+\frac{g}{2{\mathrm{\Omega }}_{g}}\frac{{\partial }^{2}{{\boldsymbol{A}}}}{\partial {T}^{2}}\quad. \end{split} $ (2)

In formula (2), $ \alpha $ is the loss in the fiber, g is the gain coefficient of the active fiber, β₃ is the third-order dispersion, and $ {\mathrm{\Omega }}_{g} $ is the gain bandwidth.

In the mode-locked laser model, solving NLSE and GLE is the prerequisite condition. The analytical solutions of equations (1) and (2) are difficult to obtain but its numerical solution can be easy to solve. At present, the split-step Fourier method is commonly used to deal with the pulse propagation in nonlinear dispersive media. This method has the advantages of high speed and high precision.

For numerical integration, it is convenient to represent (1) and (2) in the form

$ \frac{\partial {{\boldsymbol{A}}}}{\partial {\textit{z}}}=\left(\widehat{D}+\widehat{N}\right){{\boldsymbol{A}}}\quad. $ (3)

where $ \widehat{D} $ is a linear operator that takes into account the dispersion, gain and loss in the medium, and $ \widehat{N} $ is a nonlinear operator. For the GLE (2), they are given by

$ \widehat{D}=-i\frac{{\beta }_{2}}{2}\frac{{\partial }^{2}}{\partial {T}^{2}}+\frac{g-\alpha }{2}+\frac{g}{2{\mathrm{\Omega }}_{g}}\frac{{\partial }^{2}}{\partial {T}^{2}}\quad, $ (4)

$ \widehat{N}=i\gamma {\left|\boldsymbol{A}\right|}^{2}\quad. $ (5)

Dispersion and nonlinearity work together along the length of the fiber. The SSFM assumes that dispersion and nonlinear effects work independently during the propagation of the pulse within a small distance h, so as to obtain an approximate solution. More specifically, the propagation from z to z + h is done in two steps. In the first step, the nonlinear effect works alone, that is, $ \widehat{D} $ = 0 in Eq. (3). In the second step, the dispersion effect works alone, that is, $ \widehat{N} $ = 0 in Eq. (3). Mathematically it can be expressed as:

$ {{\boldsymbol{A}}}\left({\textit{z}}+h,T\right)\approx \exp\left(h\widehat{D}\right)\mathrm{exp}\left(h\widehat{N}\right){{\boldsymbol{A}}}\left({\textit{z}},T\right) \quad.$ (6)

Two-step process is repeated over the entire fiber to approximate nonlinear pulse propagation. Since the dispersion and nonlinear operators do not commute in general, the solution (6) is only an approximation to the exact solution.

In order to improve the accuracy of SSFM, it can be optimized by using the symmetrical SSFM. Suppose there are two operators $ \widehat{a} $ and $ \widehat{b} $ that exhibit non-commutation, then according to Baker-Hausdorff formula, it can be expressed as：

$ \begin{split} &\exp\left(\widehat{a}\right)\mathrm{exp}\left(\widehat{b}\right)=\\ &\exp\left[\widehat{a}+\widehat{b}+\frac{1}{2}\left[\widehat{a},\widehat{b}\right]+\frac{1}{12}\left[\widehat{a}-\widehat{b},\left[\widehat{a},\widehat{b}\right]+\cdots \right]\right]\quad, \end{split} $ (7)

where, $ \left[\widehat{a},\widehat{b}\right]=\widehat{a}\widehat{b}-\widehat{b}\widehat{a} $. Put $ \widehat{a}=h\widehat{D} $ , $ \widehat{b}=h\widehat{N} $ into the equation (7). Currently, the accuracy of the spilt-step Fourier method is O ( $ {h}^{2} $). This method improves the calculation accuracy to some extent, and the solution is as follows:

$ {{\boldsymbol{A}}}\left({\textit{z}}+h,T\right)=\exp\left(\frac{h}{2}\widehat{D}\right)\mathrm{exp}\left(h\widehat{N}\right)\exp\left(\frac{h}{2}\widehat{D}\right){{\boldsymbol{A}}}\left({\textit{z}},T\right) \quad.$ (8)

The biggest advantage of this method is that it no longer simply separates linear and nonlinear calculations but includes nonlinear effects in the middle part of the interval. This processing method can improve the accuracy to 3rd order of step size h. Because the operators in Eq. (8) show a symmetric form, we usually call this algorithm as the symmetric SSFM. This method plays an important role in solving many optical problems. In addition, although this symmetric SSFM has advantages over other general methods in numerical solution, it should be noted that in this method for choosing appropriate parameters. The size of the window should be appropriate so that the pulses are still within this window after the transmission is completed, otherwise the pulses will be superimposed. Besides, to ensure the accuracy of the results in the SSFM with constant step, the selection of the step size should be moderate. If the selection of the step size is too long, it will lead to insufficient accuracy. If the selection of the step size is too small, it will increase the amount of calculation.

We combine global-error-local-energy (GELE) method with SSFM into SSFM-GELE method, which is proposed to solve the problem related with computational efficiency. In our laser model, a local energy can be used as a criterion to estimate the total error, so it can be used to control step size of the laser in every round trip. This idea about the local energy comes from the stable pulse energy circulated in the laser system.

When simulating the GLE equation in the active fiber, the energy contained in the pulse is calculated. We use this calculated energy to limit it to a certain range, so that the energy of the pulse within a certain step size h is not particularly high, then the overall error of the final pulse is not too much compared with the fine solution.

In the SSFM-GELE method, the step size is adaptively chosen so that the calculated δ for each turn is limited to a specified range. The δ is defined by

$ \delta =\int {\left|{u}_{{n}_{i}}\left(t\right)\right|}^{2}{\rm{d}}t-\int {\left|{u}_{({n-1})_{i}}\left(t\right)\right|}^{2}{\rm{d}}t\quad, $ (9)

$ {n}_{i}=\frac{{L}_{{\rm{AF}}}}{{h}_{i}} \quad,$ (10)

where i is the number of round trips, the total number of round trips is a constant. ${L}_{{\rm{AF}}}$ is the length of active fiber. $ {h}_{i} $ and $ {n}_{i} $ are respectively the step size and step number of i-th round. The selection of the step size adjustment coefficient will affect the efficiency of the method.

In our method, the step size is chosen by keeping the δ within a specified range (1/2 $ {\delta }_{\rm G} $, $ {\delta }_{\rm G} $), where $ {\delta }_{\rm G} $ is the goal energy increment standard of a step size, using $ {\delta }_{\rm G} $ to control the change of step size can make the total error of the target be controlled within the goal value. If δ is in the range ( $ {\delta }_{\rm G} $, $ {2\delta }_{\rm G} $), the step size is divided by a factor of $ {2}^{\tfrac{1}{3}} $ for the next step. If ${\delta < \dfrac{1}{2\delta }}_{\rm G}$, the step size is multiplied by a factor of $ {2}^{\tfrac{1}{3}} $ for the next step. The value of $ {\delta }_{\rm G} $ should be adjusted in the simulation according to the acceptable global error.

Using the global error to control the adaptive change of the step size, which could increase the computational time. Here, we choose δ the local energy increment as a substitute of the global error, which can avoid unnecessary computational time consumption. Furthermore, when solving the GLE, the local pulse energy in the active fiber must be calculated. Thus, the proposed SSFM-GELE method as a global error optimization method can be utilized in the laser model.

In this section, the performance of the adaptive step size algorithm described in Section 2 is compared with the traditional constant step size algorithm. To compare the different methods, first a fine solution $ {A}_{{\rm{ref}}} $ is computed with SSFM-constant method using a fine small step size called as SSFM-fine, which step size is smaller than the smallest step size chosen by the adaptive method. Next, the numerical solutions $ {A}_{{\rm{calc}}} $ is obtained by various kinds of methods and the global error ξ is defined by：

$ \xi =\frac{\displaystyle\sum _{k=1}^{N}\left|{\left|{A}_{{\rm{calc}}}\right|}^{2}-{\left|{A}_{{\rm{ref}}}\right|}^{2}\right|/N}{\mathrm{max}\left({\left|{A}_{{\rm{ref}}}\right|}^{2}\right)} \quad,$ (11)

where N is the number of frequency grid points.

In an all-normal dispersion laser, high-energy dissipative soliton pulses can be output. In addition to the balance of nonlinear effects and dispersion, dissipative soliton also needs to consider the balance of gain and loss. These multiple effects work together to produce a fixed local solution. It is worth noting that this fixed local solution is independent of the initial pulse condition. Dissipative solitons are very sensitive to numerical errors, and thus it is requiring an efficient adaptive algorithm.

Here an all-normal dispersion passive mode-locked erbium-doped fiber laser is modelled for outputting a dissipative soliton, where its structure is shown in Fig. 1.

The fast saturable absorption model is applied to simulate Saturable Absorber (SA), which absorption coefficient is described as:

$ \alpha \left(t\right)=\frac{{\alpha }_{0}}{1+\dfrac{{E}_{{\rm{p}}}}{{E}_{{\rm{sat}}}}}+{\alpha }_{ns}\quad, $ (12)

where, $ {\alpha }_{0} $ is the modulation depth, $ {E}_{{\rm{p}}} $ is the pulse energy, $ {E}_{{\rm{sat}}} $ is the saturation energy of SA, and $ {\alpha }_{{\rm{ns}}} $ is the unsaturated loss. Due to the important role of the saturable absorber in the evolution of dissipative solitons, the modulation depth $ {\alpha }_{0} $ is set to 30%. Recovery time is 500 fs, $ {E}_{{\rm{sat}}} $ is 9.6133 pJ, and the ${\alpha }_{ns}$ is 5%. The GVD parameter of 4 m long DCF and 1 m long EDF is −36 ps/(nm·km), nonlinear coefficient is 0.004 W⁻¹·m⁻¹. The gain fiber bandwidth was set to 50 nm, the central wavelength was 1550 nm, gain coefficient is 4, and the saturation energy was 8 nJ. Fibers are described using the NLSE and GLE mentioned above, respectively. The Gaussian-shaped filter has a bandwidth of 12 nm and a central wavelength of 1550 nm. The coupler power output ratio is 30%. Filter and coupler are directly multiplied by their transfer functions in frequency domain and time domain respectively. The grid points are 8 192. The pulse evolves from initial noise, and after many cycles, the pulse evolves into a stable dissipative soliton. In the simulation, to verify the advantage of our improved method, we choose a large initial step size. The number of cycles is set to 300 to enable stable mode-locking to be established.

Before simulating with the SSFM-GELE method, it is necessary to observe the relationship between global error and the parameter δ. Ideally, the global error should be linear with the local energy increment δ. The curve in Fig. 2 is shown that the global error is related with the local energy increment δ. It is approximated as a straight line, which proves that they are indeed linear relationship. The local energy increment δ can be as a good indicator to evaluate the global error. To obtain a low global error, the local energy increment δ should be controlled as a befitting range. The parameter $ {\delta }_{G} $ it is utilized to adjust the local energy increment δ for improving the computing accuracy of the simulation method.

Fig. 3 (color online) shows the evolution of pulse energy in a local region of the active fiber. The local energy will rise a lot, due to the relatively large initial step size used in our proposed adaptive method. However, due to the limitation of δ, the local pulse energy is quickly optimized to the level of reference solution obtained by the constant method at the 42nd round. That is, δ is a parameter used to control the step size. The result for red line is obtained by SSFM with the extremely small constant step size of 0.0001 m as the fine solution, which can be called as SSFM-fine by us. Two simulation methods about SSFM-GELE and SSFM-fine with small step have been compared in the following.

Firstly, the pulse spectra are obtained by the two simulation methods, as shown in Fig. 4. In Fig. 4, the spectral shapes are rectangular, flat in the middle and steep on both sides, which proves that obtained dissipative soliton mode-locked pulses accord with the structure design of the mode-locked fiber laser. In Fig. 4(a) spectral bandwidth is 13.4 nm, and in Fig. 4(b) spectral bandwidth is 13.1 nm. Their central wavelengths are both at 1550 nm. The difference between them is very small. It can be concluded that the adaptive method (SSFM-GELE) could present the same accuracy of the spectra solution with small step constant method (SSFM-fine).

The temporal pulses have been obtained by the two simulation methods, as shown in Fig. 5. The difference in the simulation results between the solution obtained by the SSFM-GELE method and SSFM-fine is shown in the inserted picture in Fig. 5. The pulse widths obtained by the SSFM-GELE method and SSFM-fine method are 11.83 ps and 11.72 ps, respectively. Their time-bandwidth products are 158.5 and 153.5, implying that the pulses have a very large positive chirp. Two intervals basically coincide so the adaptive method brings little error to the result. Currently, the global error between them is 8.261 4×10⁻⁴.

Next, we observed the evolution of the mode-locked pulses when using different methods. Fig. 6(a) and 6(c) (color online) show the temporal evolution of the dissipative soliton generation process calculated by the SSFM-GELE method. Fig. 6(b) and 6(d) (color online) show the temporal evolution by SSFM-fine. By comparison, the proposed adaptive method has little effect on the establishment process of the dissipative soliton. The pulse is gradually narrowed and amplified from the initial noise, and finally mode-locked in the dynamic balance of gain and loss, dispersion and nonlinearity, and a stable dissipative soliton is outputted.

Finally, to analyse the advantage of our proposed method, the computation time and computation accuracy have been calculated to judge the simulation method. The global error tolerable depends largely on the nature of the method being studied. For the type of dissipative solitons that is simulated here, global error of $ {10}^{-3} $ or smaller are required to numerically achieve stable results that are adequate for studying the physics of the process or for comparison with experimental results. In Fig. 7, it can be clearly seen that the efficiency of the SSFM-GELE method is always better than that of the SSFM with a constant step size method. As expected, our proposed adaptive methods always save unnecessary computation time. We can also see that when the calculation accuracy is higher, the SSFM-GELE method saves more time than the normal SSFM with constant step method. The global error calculated by the SSFM-GELE method can reach 2.81×10⁻⁶ with 255 s, while the SSFM with small constant step size method needs to calculate 3855 s to achieve the same order of computation accuracy. However, limited by the accuracy of the SSFM algorithm, the optimization accuracy of the SSFM-GELE method is also limited. Thus, higher accuracy is required, higher-order simulation algorithms would also be needed to combine with our SSFM-GELE method.

In this paper, a SSFM-GELE method has been proposed to simulate the numerical solution of passive mode-locked fiber laser. It is combined with the SSFM method to solve problems involving both NLSE and GLE equations. Such method is utilized to optimize the step size by limiting the energy in the local area of the active fiber, so that the global error between the output pulse of each turn and the reference solution decrease dreamily. It solves the conflict problem between long computational time and high accuracy in passively mode-locked fiber laser simulation.

The performance of the mode-locked fiber laser based on the SSFM-GELE method and SSFM-constant method have been compared. The comparation results show that the SSFM-GELE method always outperforms SSFM-constant method in time efficiency because the characteristic of our method to dynamically adjust the step size leads to high computational efficiency. The global error calculated by the SSFM-GELE method can reach 2.81×10⁻⁶ with 255 s, while the SSFM with small constant step size method needs to calculate 3855 s to achieve the same order of computation accuracy. At the same time, the SSFM-GELE method also ensures that the calculation accuracy would not be lost. Our proposed method is not only combined with the SSFM method, but also has the potential to combine with other high-order algorithms, such as RK4IP, Adams, predictor–corrector, etc for improving the accuracy.