Comparison of Models of Fast Saturable Absorption in Passively Modelocked Lasers

Fast saturable absorbers (FSAs) play a critical role in stabilizing many passively modelocked lasers. The most commonly used averaged model to study these lasers is the Haus modelocking equation (HME) that includes a third-order nonlinear FSA. However, it predicts a narrow region of stability that is inconsistent with experiments. To better replicate the laser physics, averaged laser models that include FSAs with higher-than-third-order nonlinearities have been introduced. Here, we compare three common FSA models to each other and to the HME using the recently-developed boundary tracking algorithms. The three FSA models are the cubic-quintic model, the sinusoidal model, and the algebraic model. We find that all three models predict the existence of a stable high-energy solution that is not present in the HME and have a much larger stable operating region. We also find that all three models predict qualitatively similar stability diagrams. We conclude that averaged laser models that include FSAs with higher-than-third-order nonlinearity should be used when studying the stability of passively modelocked lasers.


Introduction
Over the past few decades, ultra-short optical pulses that are produced by passively modelocked lasers have been used in a broad range of fields [1,2]. Saturable absorption plays a critical role in stabilizing in these lasers and determining their pulse shape [3,4]. There has been considerable controversy over how best to model the saturable absorber [5][6][7][8][9][10][11][12], and resolution of this issue has been hampered until recently by lack of a theoretical tool that is capable of accurately determining the stability boundaries in the parameter space for the different models of saturable absorption.
The most commonly used averaged model of passively modelocked lasers with a fast saturable absorber (FSA) and a slow saturable gain is the Haus modelocking equation (HME) [13]. This model has been used to study solid-state lasers [4,14,15], fiber lasers [9,[16][17][18], as well as semiconductor lasers [19]. The HME may be written as [20] where u(t, z) is the complex field envelope, t is the retarded time, z is the propagation distance, φ is the phase rotation per unit length, l is the linear loss coefficient, g(|u|) is the saturated gain, β is the group velocity dispersion coefficient, γ is the Kerr coefficient, ω g is the gain bandwidth, and f sa (|u|) is the fast saturable absorption. It is common in studies of the HME to set dφ /dz = 0 [6,13,14,16,21], in which case the stationary pulse solution-which we also refer to as the equilibrium solution-rotates at a constant rate as a function of z. In a stability analysis, it is more convenient to start from a stationary solution of Eq. (1), in which case φ = φ 0 must be found along with u 0 (t) [22], which is what we will do here. In the HME, it is assumed that the gain response of the medium is much longer than the roundtrip time T R , in which case the saturable gain becomes g(|u|) = g 0 1 + P av (|u|)/P sat , where g 0 is the unsaturated gain, P av (|u|) is the average power, and P sat is the saturation power. We may write P av (|u|) = T R /2 −T R /2 |u(t, z)| 2 dt/T R . In this article, we focus on the anomalous chromatic dispersion regime in which β < 0. In the HME, the FSA function f sa has the simple form The HME has been successful in exploring many of the qualitative features of passively modelocked lasers. However, the predictions of the HME for the instability thresholds are unrealistically pessimistic [21,23,24]. One particular reason is that Eq. (3)-which behaves as an ideal saturable absorber-provides unlimited gain, which grows cubically as the input pulse amplitude grows. This mechanism in theory leads to an infinite growth of the pulse energy and hence a saddle-node instability when the cubic coefficient δ becomes sufficiently large [21].
In order to better replicate the laser physics, different models of the FSA have been developed. However, there has been controversy concerning which FSA model is "best" to use in the sense that it best reproduces the broad stability region that is observed in experiments [7,24,25]. In this article, we address this issue using the recently developed boundary-tracking algorithm [22,26] to examine three common FSA models, which all have a quintic nonlinearity at second-lowest order in a Taylor expansion of the intensity. We show that while these models are quantitatively different, they produce similar stability diagrams in the parameter space. Consequently, none of these models is to be preferred on the basis of their qualitative features, and only careful comparison with experiments with an appropriate selection of the model parameters can distinguish among them. However, all three models predict the existence of a stable high-energy solution that is not present in the HME, and, as a consequence, all three models have a much larger stable operating region in the parameter space.
The reminder of this article is organized as follows: We introduce and compare an extended modelocking equation with different FSA models in Sec. 2. In Sec. 3, we compare the dynamical structure of the stability regions of the extended model equation with the three different FSA models. We show two cases in which singular solutions exist but remain stable. In Sec. 4, we compare the computational pulses with experimental results, and we conclude that modelocking models with FSA models that includes higher-than-third order terms should always be preferred to the HME when studying the stability of passively modelocked lasers with averaged models.

Models of Fast Saturable Absorption
Two physics-based models of fast saturable absorption (FSA) are commonly used. The first and oldest of these models is the algebraic model. This model has been applied particularly to solid state lasers [27], but also to analyze the noise and stability of passively modelocked lasers in general [5,7,28]. In this model, it is assumed that the absorbing medium is a two-level system in which the response time of the medium is fast compared to the pulse duration, so that the population of the upper state is proportional to |u(t)| 2 . In this case, we find [5,7,13,20,28 where ∂ u/∂ z| ab is the contribution to the loss from the absorbing mechanism, f 0 is a constant, and P ab is the saturation power of the absorber. A common way of simplifying the modelocking model is to incorporate the unsaturated loss s 0 into the linear loss l [13,20], i.e., l/2 = α + f 0 , where α donates the total loss that is not due to the material absorber, such as losses from the end mirrors and couplers. We then obtain When the saturation is relatively weak, i.e., |u(t)| 2 P ab , we find By keeping only the leading-order term, we can obtain a cubic model, i.e., Eq. (3), identifying δ = f 0 /P ab . Similarly, by truncating this expansion at the order |u| 4 , we obtain a cubic-quintic model where σ = f 0 /P 2 ab . The second physics-based model assumes that the FSA is due to a combination of nonlinear polarization rotation and polarization selective elements that attenuate low intensities more than high intensities. In this model, we have [9][10][11][12] where the constants f 0 , f 1 , µ, and ν depend on the settings of the polarization selective elements and the amount of nonlinear polarization rotation in one pass through the laser. If we may assume µ|u| 2 1, then The combination − f 0 + f 1 cos ν may be absorbed into the total linear loss. We then find that the nonlinear terms in Eq. (9) becomes similar to the Taylor expansion in Eq. (6). We note that this sinusoidal model also applies to lasers that use a nonlinear optical loop mirror (NOLM) as the FSA [29].
In this article, we study quantitatively how the different FSA models in Eqs. (5), (7), and (8) affect the stability of the modelocking equation of Eq. (1). To facilitate quantitative comparison, we match up the Taylor series of these three models through order |u| 4 by assigning in which ν = π/4 and µ = 2σ /δ . We can then write the nonlinear gain functions of the three FSA models as where σ > 0, and the sub-indices "cq", "ag", and "sn" represents the cubic-quintic model, the algebraic model, and the sinusoidal model, respectively. We refer to δ as the cubic coefficient and σ as the quintic coefficient. Compared with Eqs. (3), (5), and (8), we remove the dependence between δ and σ in Eq. (10) to broaden the generality of comparison. Here, we refer to Eq. (1) with different FSA models as the generalized Haus modelocking equation (GHME).
To compare the three FSA models, we show in Fig. 1 the nonlinear gain f sa (|u|) of the three models, in which the input |u| is assumed to be the instantaneous pulse amplitude. We show the nonlinear gain with four sets of (σ , δ ) values. For all three models, we find that the nonlinear gain is almost identical when |u| is small, and becomes increasingly different as |u| grows. The difference grows also when δ and σ increase for a given |u|. We also find that the nonlinear gain is the greatest and grows monotonically as |u| increases when using the algebraic model. Meanwhile, the nonlinear gain saturates and decreases when |u| becomes sufficiently large with either the cubic-quintic model or the sinusoidal model. Hence, the algebraic model will lead to greater nonlinear gain and stationary pulses with higher energies than is the case with the other two models. The nonlinear gain is smallest with the sinusoidal model.

Stability of the Generalized Haus Modelocking Equations
In our analysis, we use the normalization and the parameters for a soliton laser in [6,21]. The pulse u is normalized with respect to the peak power of the electrical field, the propagation variable z is normalized to the dispersion length, and the retarded time t is normalized to the pulse width. Both u and φ become become unitless with this normalization. The set of parameters is listed in Table 1. (In Table. 1 of [22], the value of ω g was given as √ 10/2 by mistake.) Table 1. Normalized values of parameters.

The cubic-quintic FSA model
In prior work [22, 26], we developed boundary tracking algorithms, and we found the stability regions of the GHME with the cubic-quintic model of Eq. (10a). The basic approach of the boundary tracking algorithms is that as parameters vary, we first find a stationary solution [u 0 (t), φ 0 ] by solving a root-finding problem. Then we determine the stability of this stationary solution by calculating the eigenvalues of the linearized evolution equations. We show two examples of the distribution of the eigenvalues on the complex plane in Appendix A. The stationary solution is unstable if any of the eigenvalues have a positive real part. We then repeat these computations while tracking the stability boundary in the parameter space. In this section, we review the results and show the stability structure in Fig. 2. We have discovered a rich dynamical structure. The δ -axis of Fig. 2 corresponds to the HME case in which the quintic coefficient is zero (σ = 0). A known analytical solution is where A h , β h , and t h are constants that can be derived from the system parameters. The stationary solution u h (t) is stable when 0.01 < δ < 0.0348 [21].

Stability of the Generalized Haus Modelocking Equations
In our analysis, we use the normalization and the parameters for a soliton laser in [6,21]. The pulse u is normalized with respect to the peak power of the electrical field, the propagation variable z is normalized to the dispersion length, and the retarded time t is normalized to the pulse width. Both u and φ become become unitless with this normalization. The set of parameters is listed in Table 1. (In Table. 1 of [22], the value of ω g was given as √ 10/2 by mistake.) Table 1. Normalized values of parameters.

The cubic-quintic FSA model
In prior work [22, 26], we developed boundary tracking algorithms, and we found the stability regions of the GHME with the cubic-quintic model of Eq. (10a). The basic approach of the boundary tracking algorithms is that as parameters vary, we first find a stationary solution [u 0 (t), φ 0 ] by solving a root-finding problem. Then we determine the stability of this stationary solution by calculating the eigenvalues of the linearized evolution equations. We show two examples of the distribution of the eigenvalues on the complex plane in Appendix A. The stationary solution is unstable if any of the eigenvalues have a positive real part. We then repeat these computations while tracking the stability boundary in the parameter space. In this section, we review the results and show the stability structure in Fig. 2. We have discovered a rich dynamical structure. The δ -axis of Fig. 2 corresponds to the HME case in which the quintic coefficient is zero (σ = 0). A known analytical solution is where A h , β h , and t h are constants that can be derived from the system parameters. The stationary solution u h (t) is stable when 0.01 < δ < 0.0348 [21].
[1] P Figure 2. The stability regions of the GHME with a cubic-quintic saturable absorber f sa,cq (|u|). The stability boundaries are marked by three curves, C 1 , C 2 , and C 3 . This figure reproduces Fig. 16  When σ = 0, the solution in Eq. (11) does not hold any more, and we have found two stable numerical pulse solutions of the GHME: a low-amplitude solution (LAS) and a high-amplitude solution (HAS). The LAS is stable in the blue-hatched region that is marked with [2 l ] as shown in Fig. 2, and it becomes unstable in region [1] (below the curve C 1 ), where the continuous modes become unstable via a Hopf bifurcation or an essential instability [30]. The amplitude eigenmode becomes unstable when we cross C 3 from region [2 l ], which corresponds to a saddlenode instability. Meanwhile, there also exists a second stationary solution [22], which we refer to as the high-amplitude solution (HAS). The HAS is stable in the red-hatched region which is marked with [2 h ], and its stability region is lower-bounded by C 2 , below which the amplitude eigenmode becomes unstable via a saddle-node bifurcation. The LAS and the HAS coexist and remain stable in region [3] which is bounded by C 2 and C 3 . A region like that of region [3] in which the LAS and the HAS coexist has recently been experimentally confirmed [31]. The LAS and the HAS merge into one single solution in region [2 h/l ] in which it remains stable. The HAS does not exist for the HME, which is the reason the behavior of the GHME is qualitatively different from the HME. We also find that the HAS remains stable until δ increases up to ∼9.51, which is almost a factor of 280 greater than the HME's stability limit [22]. For a given δ , the HAS also remains stable for any value of σ as long as σ > 0, although the stable equilibrium pulse becomes increasingly peaked and narrow as σ approaches zero [32]. We summarize which solution becomes unstable on the curves C 1 , C 2 , and C 3 and the instability mechanism in Table 2.
Curve Solution Instability Mechanism Table 2. Instability mechanisms of the GHME shown in Fig. 2, where LAS represents the low-amplitude solution, and HAS represents the high-amplitude solution.
In Fig. 3, we show an example of the pulse profiles of both the LAS and the HAS when σ = 0.004 and δ = 0.036, which is at the point P in Fig. 2. We write the stationary pulse profile When σ = 0, the solution in Eq. (11) does not hold any more, and we have found two stable numerical pulse solutions of the GHME: a low-amplitude solution (LAS) and a high-amplitude solution (HAS). The LAS is stable in the blue-hatched region that is marked with [2 l ] as shown in Fig. 2, and it becomes unstable in region [1] (below the curve C 1 ), where the continuous modes become unstable via a Hopf bifurcation or an essential instability [30]. The amplitude eigenmode becomes unstable when we cross C 3 from region [2 l ], which corresponds to a saddlenode instability. Meanwhile, there also exists a second stationary solution [22], which we refer to as the high-amplitude solution (HAS). The HAS is stable in the red-hatched region which is marked with [2 h ], and its stability region is lower-bounded by C 2 , below which the amplitude eigenmode becomes unstable via a saddle-node bifurcation. The LAS and the HAS coexist and remain stable in region [3] which is bounded by C 2 and C 3 . A region like that of region [3] in which the LAS and the HAS coexist has recently been experimentally confirmed [31]. The LAS and the HAS merge into one single solution in region [2 h/l ] in which it remains stable. The HAS does not exist for the HME, which is the reason the behavior of the GHME is qualitatively different from the HME. We also find that the HAS remains stable until δ increases up to ∼9.51, which is almost a factor of 280 greater than the HME's stability limit [22]. For a given δ , the HAS also remains stable for any value of σ as long as σ > 0, although the stable equilibrium pulse becomes increasingly peaked and narrow as σ approaches zero [32]. We summarize which solution becomes unstable on the curves C 1 , C 2 , and C 3 and the instability mechanism in Table 2.
Curve Solution Instability Mechanism Table 2. Instability mechanisms of the GHME shown in Fig. 2, where LAS represents the low-amplitude solution, and HAS represents the high-amplitude solution.
In Fig. 3, we show an example of the pulse profiles of both the LAS and the HAS when σ = 0.004 and δ = 0.036, which is at the point P in Fig. 2. We write the stationary pulse profile as where θ 0 (t) is the phase across the pulse. The pulse envelopes of both the LAS and the HAS have a nearly hyperbolic secant profile, in which the amplitude decays exponentially as t → ±∞.
Compared to the LAS, the amplitude of the HAS is visibly higher, and the pulse width is smaller. In addition, the HAS has a stronger chirp than does the LAS.

Comparison of stability with different FSA models
In Fig. 4, we compare regions of stability for the three different FSA models. We find that the dynamical structures are qualitatively similar for all three models. For regions near the origin, all three models feature the three characteristic curves C 1 , C 2 , and C 3 described in Table 2.  We first compare the location of C 1 , which marks the lower boundary where the LAS becomes unstable due to the continuous modes. When the nonlinear gain in the GHME increases, the energy of the stationary pulse solution grows, and the linear gain g(|u|) becomes increasingly saturated. Since the stability condition for continuous waves is [g(|u|) − l] < 0 [21, 22], as where θ 0 (t) is the phase across the pulse. The pulse envelopes of both the LAS and the HAS have a nearly hyperbolic secant profile, in which the amplitude decays exponentially as t → ±∞.
Compared to the LAS, the amplitude of the HAS is visibly higher, and the pulse width is smaller. In addition, the HAS has a stronger chirp than does the LAS.

Comparison of stability with different FSA models
In Fig. 4, we compare regions of stability for the three different FSA models. We find that the dynamical structures are qualitatively similar for all three models. For regions near the origin, all three models feature the three characteristic curves C 1 , C 2 , and C 3 described in Table 2. as where θ 0 (t) is the phase across the pulse. The pulse envelopes of both the LAS and the HAS have a nearly hyperbolic secant profile, in which the amplitude decays exponentially as t → ±∞.
Compared to the LAS, the amplitude of the HAS is visibly higher, and the pulse width is smaller. In addition, the HAS has a stronger chirp than does the LAS.

Comparison of stability with different FSA models
In Fig. 4, we compare regions of stability for the three different FSA models. We find that the dynamical structures are qualitatively similar for all three models. For regions near the origin, all three models feature the three characteristic curves C 1 , C 2 , and C 3 described in Table 2.  We first compare the location of C 1 , which marks the lower boundary where the LAS becomes unstable due to the continuous modes. When the nonlinear gain in the GHME increases, the energy of the stationary pulse solution grows, and the linear gain g(|u|) becomes increasingly saturated. Since the stability condition for continuous waves is  We first compare the location of C 1 , which marks the lower boundary where the LAS becomes unstable due to the continuous modes. When the nonlinear gain in the GHME increases, the energy of the stationary pulse solution grows, and the linear gain g(|u|) becomes increasingly saturated. Since the stability condition for continuous waves is [g(|u|) − l] < 0 [21, 22], the continuous wave perturbations become increasingly stable. Hence, a large nonlinear gain in the FSA model implies a large stability region. We recall that for fixed values of (σ , δ ), the nonlinear gain is largest for the algebraic model and smallest for the sinusoidal model. This behavior is consistent with what we show in Fig. 4. With the same value of σ , the value of δ on C 1 is smallest for the algebraic model, and then is larger for the cubic-quintic model, and is largest for the sinusoidal model.
Next we analyze the stability limits of the saddle-node instability that leads to an amplitude explosion [21,22]. The saddle-node instability occurs when the nonlinear gain becomes so large that it cannot be compensated by the nonlinear loss. Thus, with the same values of σ and δ , an FSA model that provides a relatively small nonlinear gain tends to be more stable against the occurrence of the saddle-node instability. In Fig. 4, the onset of the saddle-node instability of the LAS is marked by C 3 . We observe that with the same value of σ , the stability limit of δ is largest with for sinusoidal model, next largest for the cubic-quintic model, and is smallest for the algebraic model. the continuous wave perturbations become increasingly stable. Hence, a large nonlinear gain in the FSA model implies a large stability region. We recall that for fixed values of (σ , δ ), the nonlinear gain is largest for the algebraic model and smallest for the sinusoidal model. This behavior is consistent with what we show in Fig. 4. With the same value of σ , the value of δ on C 1 is smallest for the algebraic model, and then is larger for the cubic-quintic model, and is largest for the sinusoidal model.
Next we analyze the stability limits of the saddle-node instability that leads to an amplitude explosion [21,22]. The saddle-node instability occurs when the nonlinear gain becomes so large that it cannot be compensated by the nonlinear loss. Thus, with the same values of σ and δ , an FSA model that provides a relatively small nonlinear gain tends to be more stable against the occurrence of the saddle-node instability. In Fig. 4, the onset of the saddle-node instability of the LAS is marked by C 3 . We observe that with the same value of σ , the stability limit of δ is largest with for sinusoidal model, next largest for the cubic-quintic model, and is smallest for the algebraic model. In Fig. 5 we show the pulse profiles of the LAS of all three FSA models when σ = 0.004. Here, the parameter A 0 is the peak amplitude of the pulse, τ 0 is the FWHM width of the pulse solution, and b 0 is the chirp, We see that the pulse profiles with the three models are nearly identical for the LAS. As the cubic gain coefficient δ increases, we find that the pulse amplitude A 0 increases, the pulse width τ 0 decreases, and both the rate of phase rotation φ 0 and the chirp b 0 increase. Among the three models, the pulse for the algebraic model has the highest amplitude and narrowest width with a fixed (σ , δ ), while the pulse for the sinusoidal model has the smallest amplitude and the largest pulse width. In Fig. 1, we see that the curves C 1 and C 3 are quantitatively close when 0 < σ < 0.006, which is because the nonlinear gain is nearly identical when |u| is small. The small differences in the nonlinear gain of the FSA models as shown in Fig. 1 are consistent with the small differences in the pulse profiles as shown in Fig. 5(a), in which the peak amplitude A 0 is largest and the pulse width τ 0 is smallest for the algebraic model, while A 0 and τ 0 are smallest and largest respectively for the sinusoidal model. In Fig. 6, we show the variation of pulse profiles for the HAS with the three models when σ = 0.004 and δ varies. As δ decreases, the peak amplitude A 0 and then the peak quintic loss decreases, and a saddle-node bifurcation eventually occurs when the quintic loss becomes insufficient to offset the nonlinear gain. Here, the changes of the nonlinear loss and gain are dominated by the change of amplitude of the stationary pulse, which is different from the case In Fig. 5 we show the pulse profiles of the LAS of all three FSA models when σ = 0.004. Here, the parameter A 0 is the peak amplitude of the pulse, τ 0 is the FWHM width of the pulse solution, and b 0 is the chirp, We see that the pulse profiles with the three models are nearly identical for the LAS. As the cubic gain coefficient δ increases, we find that the pulse amplitude A 0 increases, the pulse width τ 0 decreases, and both the rate of phase rotation φ 0 and the chirp b 0 increase. Among the three models, the pulse for the algebraic model has the highest amplitude and narrowest width with a fixed (σ , δ ), while the pulse for the sinusoidal model has the smallest amplitude and the largest pulse width. In Fig. 1, we see that the curves C 1 and C 3 are quantitatively close when 0 < σ < 0.006, which is because the nonlinear gain is nearly identical when |u| is small. The small differences in the nonlinear gain of the FSA models as shown in Fig. 1 are consistent with the small differences in the pulse profiles as shown in Fig. 5(a), in which the peak amplitude A 0 is largest and the pulse width τ 0 is smallest for the algebraic model, while A 0 and τ 0 are smallest and largest respectively for the sinusoidal model. In Fig. 6, we show the variation of pulse profiles for the HAS with the three models when σ = 0.004 and δ varies. As δ decreases, the peak amplitude A 0 and then the peak quintic loss decreases, and a saddle-node bifurcation eventually occurs when the quintic loss becomes insufficient to offset the nonlinear gain. Here, the changes of the nonlinear loss and gain are dominated by the change of amplitude of the stationary pulse, which is different from the case of the LAS where the nonlinear terms are mainly impacted by the variation of the nonlinear of the LAS where the nonlinear terms are mainly impacted by the variation of the nonlinear coefficients σ and δ . Thus, for the HAS to remain stable, the pulse solution must have a sufficiently large peak amplitude. By comparing Figs. 6(a) and 6(c) with Fig. 5(a), we find that for each of the three FSA models, the peak amplitude of the HAS is larger than that of the LAS, and hence the peak quintic loss from σ |u| 4 u is higher than that of the LAS, which explains why the HAS remains stable at large values of δ when the LAS ceases to exist. We find in Figs. 6(a) and 6(c) that the cubic-quintic model generates pulses with higher amplitude and narrower width than does the sinusoidal model, while the algebraic model leads to pulse profiles that are significantly larger in amplitude and narrower in pulse width than for the other two FSA models. This behavior occurs because the nonlinear gain of the algebraic model increases monotonically with the input, while the nonlinear gain of the other two FSA models saturates when the input becomes sufficiently large, which is shown in Fig. 1. Hence, as shown in Fig. 4, with a fixed value of σ , the value of δ on C 2 is smallest for the algebraic model, next-smallest for cubic-quintic model, and is largest for the sinusoidal model. We find that the rate of phase rotation φ 0 increases, but the chirp b 0 decreases, as δ increases, which is different from the profiles of the LAS, as can be seen comparing Fig. 5(b) to Figs. 6(b) and 6(d). Both the rate of phase rotation and the chirp with the algebraic model are significantly larger than is the case for the other two models, as shown in Fig. 6(d). For both the LAS and the HAS with any of the FSA models, we find that the rate of phase rotation φ 0 grows together with the peak amplitude A 0 , which implies that these stationary solutions are qualitatively similar to the soliton solution of the nonlinear Schödinger equation [33].

Stability of the HAS as σ → 0
When obtaining the stability structure in Figs. 2 and 4, the boundary tracking algorithm is not able to proceed when C 2 approaches the δ -axis where σ = 0. We found that the HAS becomes increasingly narrower and more energetic as σ → 0, which implies that we are approaching a singular solution. By applying asymptotic perturbation theory to the GHME with the cubic- coefficients σ and δ . Thus, for the HAS to remain stable, the pulse solution must have a sufficiently large peak amplitude. By comparing Figs. 6(a) and 6(c) with Fig. 5(a), we find that for each of the three FSA models, the peak amplitude of the HAS is larger than that of the LAS, and hence the peak quintic loss from σ |u| 4 u is higher than that of the LAS, which explains why the HAS remains stable at large values of δ when the LAS ceases to exist. We find in Figs. 6(a) and 6(c) that the cubic-quintic model generates pulses with higher amplitude and narrower width than does the sinusoidal model, while the algebraic model leads to pulse profiles that are significantly larger in amplitude and narrower in pulse width than for the other two FSA models. This behavior occurs because the nonlinear gain of the algebraic model increases monotonically with the input, while the nonlinear gain of the other two FSA models saturates when the input becomes sufficiently large, which is shown in Fig. 1. Hence, as shown in Fig. 4, with a fixed value of σ , the value of δ on C 2 is smallest for the algebraic model, next-smallest for cubic-quintic model, and is largest for the sinusoidal model.
We find that the rate of phase rotation φ 0 increases, but the chirp b 0 decreases, as δ increases, which is different from the profiles of the LAS, as can be seen comparing Fig. 5(b) to Figs. 6(b) and 6(d). Both the rate of phase rotation and the chirp with the algebraic model are significantly larger than is the case for the other two models, as shown in Fig. 6(d). For both the LAS and the HAS with any of the FSA models, we find that the rate of phase rotation φ 0 grows together with the peak amplitude A 0 , which implies that these stationary solutions are qualitatively similar to the soliton solution of the nonlinear Schödinger equation [33]. When obtaining the stability structure in Figs. 2 and 4, the boundary tracking algorithm is not able to proceed when C 2 approaches the δ -axis where σ = 0. We found that the HAS becomes increasingly narrower and more energetic as σ → 0, which implies that we are approaching a singular solution. By applying asymptotic perturbation theory to the GHME with the cubicquintic FSA model [32], we proved that the HAS continues to exist and remains stable regard-less of the magnitude of σ . Similarly, the HAS remains stable for both the sinusoidal model and the algebraic model as σ → 0 due to the similar dynamical structure that is visible in Fig. 4.

Stability of the HAS as δ becomes sufficiently large
With the cubic-quintic model, we found earlier that as δ increases, with σ fixed, an edge bifurcation occurs in which two discrete modes bifurcate out of the continuous spectrum, which then become unstable via a Hopf bifurcation at δ ≈ 9.51 [22]. Using the same approach, described in detail in Sec. 4.B.3 of [22], we find that the sinusoidal model exhibits the same qualitative behavior and becomes unstable when δ ≈ 9.26, as shown in Fig. 7. When the instability occurs, solution of the evolution equations shows that a shelf develops, which is a phenomenon that has been experimentally observed at high pump powers [34]. Because the gain turns over as |u| increases, it becomes energetically favorable for the pulse duration to increase, rather than for its amplitude to increase. quintic FSA model [32], we proved that the HAS continues to exist and remains stable regardless of the magnitude of σ . Similarly, the HAS remains stable for both the sinusoidal model and the algebraic model as σ → 0 due to the similar dynamical structure that is visible in Fig. 4.

Stability of the HAS as δ becomes sufficiently large
With the cubic-quintic model, we found earlier that as δ increases, with σ fixed, an edge bifurcation occurs in which two discrete modes bifurcate out of the continuous spectrum, which then become unstable via a Hopf bifurcation at δ ≈ 9.51 [22]. Using the same approach, described in detail in Sec. 4.B.3 of [22], we find that the sinusoidal model exhibits the same qualitative behavior and becomes unstable when δ ≈ 9.26, as shown in Fig. 7. When the instability occurs, solution of the evolution equations shows that a shelf develops, which is a phenomenon that has been experimentally observed at high pump powers [34]. Because the gain turns over as |u| increases, it becomes energetically favorable for the pulse duration to increase, rather than for its amplitude to increase. With the algebraic model, the gain never turns over as |u| increases, and a similar instability does not occur. In Fig. 6, we show that the algebraic FSA model predicts stationary pulse solutions that are significantly shorter and more energetic than is the case for the other two FSA models. When δ > 0.055, we find that as δ increases with σ fixed, the solution becomes increasingly narrow and energetic. Similar to what happens when δ is fixed and σ → 0, we find that the HAS approaches a singular solution, so that it becomes increasingly difficult to apply the boundary tracking algorithm. However, we can apply singular perturbation theory to show that a stationary solution exists at any δ and is always stable. The detailed calculation is presented in Appendix B.

Matching Experimental Results
In Sec. 3 we showed that two stable stationary pulse solutions of the GHME coexist with an arbitrarily small value of σ , as long as σ > 0, in contrast to the HME, where there is only one stable solution in a very limited region of the cubic nonlinearity δ . Since any real system is likely to have a quintic component in its saturable absorber [35], the solutions of the GHME with higher nonlinearities should provide a better approximation to the output pulses that have been observed in experiments than does the HME. In this section, we will compare the stationary solutions that are predicted by the HME and the GHME.
In Table 1, we show the parameters that we use in the comparative study, which are estimated based on the experiments. Set 1 of the parameters corresponds to a fiber laser with nonlinear Figure 7. The stability boundary of the GHME with the cubic-quintic model and the sinusoidal model at large values of δ as σ varies. The unstable region of both models lies above each curve. The pulse solution of the GHME with the algebraic model is always stable when δ increases, as we prove in Appendix B.
With the algebraic model, the gain never turns over as |u| increases, and a similar instability does not occur. In Fig. 6, we show that the algebraic FSA model predicts stationary pulse solutions that are significantly shorter and more energetic than is the case for the other two FSA models. When δ > 0.055, we find that as δ increases with σ fixed, the solution becomes increasingly narrow and energetic. Similar to what happens when δ is fixed and σ → 0, we find that the HAS approaches a singular solution, so that it becomes increasingly difficult to apply the boundary tracking algorithm. However, we can apply singular perturbation theory to show that a stationary solution exists at any δ and is always stable. The detailed calculation is presented in Appendix B.

Matching Experimental Results
In Sec. 3 we showed that two stable stationary pulse solutions of the GHME coexist with an arbitrarily small value of σ , as long as σ > 0, in contrast to the HME, where there is only one stable solution in a very limited region of the cubic nonlinearity δ . Since any real system is likely to have a quintic component in its saturable absorber [35], the solutions of the GHME with higher nonlinearities should provide a better approximation to the output pulses that have been observed in experiments than does the HME. In this section, we will compare the stationary solutions that are predicted by the HME and the GHME.
In Table 1, we show the parameters that we use in the comparative study, which are estimated based on the experiments. Set 1 of the parameters corresponds to a fiber laser with nonlinear polarization rotation [36,37], and set 2 of the parameters corresponds to a Cr:LiSAF laser that is based on Kerr-lens modelocking [38].  We show a comparison of the computational stationary pulses and the corresponding experimental results in Fig. 8. The fiber laser in [36] generates a comb output with chirped pulses that have a duration of 210 fs and a peak power of 435 W. We show in Fig. 8(a) that using the GHME we are able to obtain a computational pulse with a full-width-half-maximum (FWHM) duration of 271 fs and peak power of 421 W. We can achieve the closest match to the experimental pulse by setting δ = 0.705 kW −1 , where the stationary pulse has an FWHM duration of 32.8 fs and a peak power of 287 W. When δ increases further, the HME solution ceases to exist. The GHME solution is visibly a better match to the experimental pulse than is the HME solution. The modeling accuracy of the GHME could be further improved by a more accurate measurement of the parameters. In particular, the value of dispersion given in [36] neglects the contribution from some cavity components.  We show a comparison of the computational stationary pulses and the corresponding experimental results in Fig. 8. The fiber laser in [36] generates a comb output with chirped pulses that have a duration of 210 fs and a peak power of 435 W. We show in Fig. 8(a) that using the GHME we are able to obtain a computational pulse with a full-width-half-maximum (FWHM) duration of 271 fs and peak power of 421 W. We can achieve the closest match to the experimental pulse by setting δ = 0.705 kW −1 , where the stationary pulse has an FWHM duration of 32.8 fs and a peak power of 287 W. When δ increases further, the HME solution ceases to exist. The GHME solution is visibly a better match to the experimental pulse than is the HME solution. The modeling accuracy of the GHME could be further improved by a more accurate measurement of the parameters. In particular, the value of dispersion given in [36] neglects the contribution from some cavity components. In Fig. 8(b), we show the comparison of the computational pulse with the experimental result corresponding to the solid state laser of [38]. In the experiment, a gain-matched output coupler is used to overcome the gain filtering effect. Using the transmission profile of the output coupler given in [38], we are able to obtain accurately the pulse profile inside the laser cavity, where the pulse energy is 14.8 nJ and the FWHM width is 30 fs. We estimate the saturation power of the saturable absorption P ab is 363 kW, and the saturable loss f 0 is 3%. The system parameters are estimated following the approach in [27]. An algebraic model was used to model the FSA. By matching the cubic and the quintic coefficients, as seen from Fig. 8(b), we are able to obtain to In Fig. 8(b), we show the comparison of the computational pulse with the experimental result corresponding to the solid state laser of [38]. In the experiment, a gain-matched output coupler is used to overcome the gain filtering effect. Using the transmission profile of the output coupler given in [38], we are able to obtain accurately the pulse profile inside the laser cavity, where the pulse energy is 14.8 nJ and the FWHM width is 30 fs. We estimate the saturation power of the saturable absorption P ab is 363 kW, and the saturable loss f 0 is 3%. The system parameters are estimated following the approach in [27]. An algebraic model was used to model the FSA. By matching the cubic and the quintic coefficients, as seen from Fig. 8(b), we are able to obtain to a computational pulse of 14.9 nJ with an FWHM width 29.2 fs, where the match is excellent. By comparison, no stable solution exists for the HME when we set σ = 0.

Summary
We have compared three common models of FSA in the GHME to each other and to the HME using boundary tracking algorithms. These three FSA models are the cubic-quintic model, the sinusoidal model, and the algebraic model. For all three models of the FSA, the stability region is greatly increased relative to the HME, in which the FSA only has a cubic nonlinearity. The behaviors of these models are qualitatively similar, but quantitatively different, and the difference is more significant as the input power increases. At low pulse energies, any of these models can be used with an appropriate choice of parameters. At high pulse energies, the model must be carefully chosen to fit experimental systems, and it becomes questionable whether any averaged model can quantitatively reproduce the stability region of a real-world laser with lumped components. We will examine this question in future work.
The HME has at most only one stable stationary solution. By contrast, the GHME with any of the three FSA models can have two stable stationary solutions-an LAS that becomes equal to the stationary solution of the HME when σ → 0 and an HAS that is qualitatively different. The LAS energy is dominated by a balance between the gain filtering, linear loss, and the cubic nonlinearity of the FSA while the HAS energy is dominated by a balance between the cubic and quintic nonlinearities of the FAS.
The existence of the HAS is the key to understanding the greatly enhanced stability region that exists in the GHME and in experiments [24]. A quintic nonlinearity exists in the FSA of any real laser. Hence, we expect that the HAS will be present in any soliton laser. We also computed the predictions for the stationary solutions of the GHME and the HME to experiments. We found that the GHME provide much better agreement. We conclude that models of fast saturable absorption that include high-order nonlinearities should always be preferred to the HME when using averaged models.
The existence and the stability of the HAS sheds further ligh on the possibility of finding high-energy and narrow pulses in soliton lasers. We show in [32] that the HAS remains stable as σ → 0, while the pulse energy increases and the pulse duration decreases. Hence, the GHME has a much larger region of stability than does the HME even when σ is small. This result suggests that it may be possible to obtain stable modelocked pulses in experiments with high energy and large bandwidth if the quintic nonlinearity in the FSA can be decreased. In practice, this parameter is difficult to control, but our results provide motivation to make the attempt.
where u = u 0 + ∆u,ū = u * 0 + ∆ū, and [∆u, ∆ū] T is a perturbation to the stationary pulse [u 0 (t), u * 0 (t)] T , and the superscript T indicates the transpose. We next formulate the eigenvalue problem where ∆u = [∆u, ∆ū] T , and L is matrix form of the right hand side of Eq. (16). Here, we suppose that ∆u and λ are an eigenvector and its corresponding eigenvalue. If any eigenvalues have positive real parts, the pulse solution is unstable. A more detailed discussion based on the GHME with the cubic-quintic FSA model can be found in Sec. 3 of [41], as well as Sec. 3 of [22]-in which the stationary solutions are referred to as equilibrium solutions. In Fig. 9, we illustrate the eigenvalues' distribution on the complex plane, which we refer to as the (linearized) spectrum, for both the NLSE and the GHME. In both cases, the spectrum includes two branches of the continuous spectrum that are symmetric about the real axis and four real discrete eigenvalues that correspond physically to perturbations of the stationary solution' which has the stationary soliton solution u 0 (t) = sech(t), φ 0 = 1/2.
Linearization of the NLSE leads to where u = u 0 + ∆u,ū = u * 0 + ∆ū, and [∆u, ∆ū] T is a perturbation to the stationary pulse [u 0 (t), u * 0 (t)] T , and the superscript T indicates the transpose. We next formulate the eigenvalue problem where ∆u = [∆u, ∆ū] T , and L is matrix form of the right hand side of Eq. (16). Here, we suppose that ∆u and λ are an eigenvector and its corresponding eigenvalue. If any eigenvalues have positive real parts, the pulse solution is unstable. A more detailed discussion based on the GHME with the cubic-quintic FSA model can be found in Sec. 3 of [41], as well as Sec. 3 of [22]-in which the stationary solutions are referred to as equilibrium solutions. In Fig. 9, we illustrate the eigenvalues' distribution on the complex plane, which we refer to as the (linearized) spectrum, for both the NLSE and the GHME. In both cases, the spectrum includes two branches of the continuous spectrum that are symmetric about the real axis and four real discrete eigenvalues that correspond physically to perturbations of the stationary solution' As shown in Fig. 9(a), for the NLSE, all four discrete eigenvalues λ t , λ φ , λ f , and λ a equal zero, and the two branches of the continuous spectrum are aligned along the imaginary axis. Hence, the NLSE soliton is neutrally stable, which means that perturbations of the stationary soliton will not decay as z increases.
By comparison, when a pulse solution of the GHME is linearly stable, the real parts of all eigenvalues are negative except for λ t and λ φ , which equal 0, as shown in Fig. 9(b). These eigenvalues must equal zero because the GHME, as well as its linearization, are invariant with respect to time and phase shifts. In the case of Fig. 9(b), the spectrum indicates that any perturbation to the pulse solution-expect for those that are proportional to the eigenfunctions of the phase shift or the time shift-will decay as the pulse propagates. The pulses's central phase and As shown in Fig. 9(a), for the NLSE, all four discrete eigenvalues λ t , λ φ , λ f , and λ a equal zero, and the two branches of the continuous spectrum are aligned along the imaginary axis. Hence, the NLSE soliton is neutrally stable, which means that perturbations of the stationary soliton will not decay as z increases.
By comparison, when a pulse solution of the GHME is linearly stable, the real parts of all eigenvalues are negative except for λ t and λ φ , which equal 0, as shown in Fig. 9(b). These eigenvalues must equal zero because the GHME, as well as its linearization, are invariant with respect to time and phase shifts. In the case of Fig. 9(b), the spectrum indicates that any perturbation to the pulse solution-expect for those that are proportional to the eigenfunctions of the phase shift or the time shift-will decay as the pulse propagates. The pulses's central phase and central time undergoes a random walk in the presence of noise. In modern comb systems, one uses electronic feedback to force a phase or a time shift to decay exponentially toward a reference phase and time [24,[43][44][45]. From a theoretical standpoint, the electronic feedback control moves the two eigenvalues λ t and λ φ to the left in Fig. 9. This effect has been experimentally observed [46]. σ (δ + σ A 2 0 ) 3/2 log We now solve Eq. (26) computationally and evaluate λ a . We show the results in Fig. 10. As δ increases, we observe that A 0 increases while λ a becomes more negative. Hence, the stationary Figure 10. As the nonlinear gain coefficient δ increases, for the GHME in Eq. (18), the variations of (a) the peak amplitude of the HAS A 0 and (b) the amplitude value λ a .
We now solve Eq. (26) computationally and evaluate λ a . We show the results in Fig. 10. As δ increases, we observe that A 0 increases while λ a becomes more negative. Hence, the stationary solution u 0 is stable against perturbations that are proportional to the amplitude eigenfunction. Considering the asymptotic behavior of both the pulse profile and the nonlinear gain of f sa (|u|), we expect λ a to become more negative as δ further increases. Therefore, all eigenvalues are negative except λ t = λ φ = 0 as δ increases. This result indicates that, for the GHME with the algebraic FSA model, the HAS remains stable for large values of δ , and its stable region is only bounded below by the Hopf bifurcation limit, which is indicated by C 2 in Fig. 4.