1 Introduction

Ultra-short optical pulses with extremely high energy are of great importance for a variety of applications. They can be used in fields like medicine, ultrasensitive laser spectroscopy, micromachining, ultrahigh bit rate communication systems, metrology, and others [1, 2]. Thus, producing shorter and more powerful pulses is vital for further progress in science and technology.

The best apparatus for generating high-energy ultrashort optical pulses appears to be a passively mode-locked laser system [3, 4]. A common approach to describe such system is an average model in which the effects of discrete laser components in the cavity are averaged over one round trip [5,6,7,8]. The resulting master equation is the cubic–quintic complex Ginzburg Landau equation (CGLE) or a generalized version of it.

The cubic–quintic CGLE allows the prediction of highly complicated pulse dynamics, such as “exploding solitons” [9], self-pulsations [10], and the formation of multi-soliton complexes [11]. Investigations based on the CGLE have been crucial in the development of the concept of dissipative solitons [12,13,14,15,16,17,18,19,20,21,22,23]. Such concept is applicable for solid-state lasers [13,14,15,16] and fiber lasers as well as for the description of optical pulses in fiber-optic communication systems in which real losses are compensated with fiber optic amplifiers [17, 18].

It has been found numerically that the energy of a dissipative soliton solution of the CGLE increases indefinitely when the equation parameters converge to a given region of the parameter space [24]. Such set of parameters was called a dissipative soliton resonance (DSR) [25]. Found in the normal dispersion regime, it initially required a positive quintic reactive nonlinearity to appear [24, 25], but then DSR was also revealed in the chromatic dispersion-free (D = 0) regime, along with negative quintic reactive nonlinearity [26]. With the increase of pump power, the energy of a DSR pulse increases mainly due to the increase of the pulse width, while keeping the amplitude at a constant level [25, 26]. As a consequence, in order to obtain high-energy ultrashort pulses, a linear pulse compression technique has to be used outside the laser cavity.

In previous works, the authors have investigated a different kind of high-energy ultrashort pulses, which correspond to the very high amplitude (VHA) soliton solutions of the CGLE [27, 28]. Such VHA solutions occur due to a singularity first predicted both numerically and using the soliton perturbation theory in Refs. [29, 30], namely, as the nonlinear gain saturation effect tends to vanish. The increase in energy of these pulses is mainly due to the increase of the pulse amplitude, whereas the pulse width becomes narrower. The region of existence of VHA pulses was found in Ref. [28], considering both the normal and the anomalous dispersion regimes, but neglecting the higher-order reactive nonlinearity effect. In the same work, VHA pulses with high energy were found mainly in the normal dispersion region. This is in agreement with the large majority of experimental results, reporting the observation of high energy pulses in passively mode-locked lasers [12, 31, 32].

In this work we further investigate the VHA pulse solutions of the CGLE in the normal dispersion regime, where the pulse energy is higher, both numerically and using the method of moments (MMs) [25, 33, 34]. A special attention is paid to the impact of the higher-order reactive nonlinearity effect. Actually, some simulations suggest that the region of existence of VHA pulses can be extended into the normal and anomalous dispersion regimes, choosing an appropriate signal for the higher-order reactive nonlinearity parameter. So, we might ask: for which values of D can the VHA pulses still form? Which is the impact of the higher-order reactive nonlinearity on the pulse amplitude, spectrum, and energy?

This paper is organized as follows: in Sect. 2 we present the governing equation, the cubic–quintic complex Ginzburg–Landau equation, as well as the approximate analytical solutions of this equation obtained using the method of moments. In Sect. 3 some predictions of the methods of moments are presented and compared with the results obtained by numerical integration of the CGLE for the normal dispersion regime, both in the absence and in presence of the higher-order reactive nonlinearity. A region of existence of the VHA pulses is obtained in the plane defined by nonlinear gain and the dispersion. Section 4 summarizes the main results of this work.

2 Theory

The propagation of ultrashort pulses in the presence of spectral filtering, linear and nonlinear gain, can be described by the following CGLE:

$$ i{ }u_{Z} + \frac{D}{2}{ }u_{TT} + \left| u \right|^{2} { }u + { }\nu { }\left| u \right|^{4} { }u = i{ }\delta { }u + i{ }\beta { }u_{TT} + i\epsilon { }\left| u \right|^{2} { }u + i\mu { }\left| u \right|^{4} { }u,{ } $$
(1)

where Z is the propagation distance or the cavity round-trip number, T is the retarded time in a frame of reference moving with the group velocity, and u is the normalized complex envelope of the optical field. On the left-hand side, D represents the cavity group velocity dispersion, with D > 0 in the anomalous regime and D < 0 in the normal regime, and ν describes the active part of the reactive nonlinearity. On the right-hand side, δ represents the difference between linear gain and loss, and β accounts for spectral filtering. The terms with ε and μ (μ < 0) represent the nonlinear gain and the saturation of the nonlinear gain.

Equation (1) can be applied to the modelling of passively mode-locked lasers [15, 18,19,20,21,22]. Actually, the model is purely phenomenological, since the coefficients here could be very complicated functions of the real laser parameters and specific design factors. However, their exact representations are not really needed for understanding the laser operation at an intuitive level. The laser is assumed to be in stationary regime when all transient effects after starting the laser have vanished. Dissipative terms describe the averaged gain and loss processes taking place when the pulse moves in the cavity. Higher-order dissipative terms describe the cumulative nonlinear transmission characteristics of the cavity, including the passive mode-locking mechanism.

In order to obtain approximate stationary solutions of Eq. (1), the method of moments (MM), [25, 33, 34] will be considered. This method provides a reduction of the complete evolution problem with an infinite number of degrees of freedom to the evolution of a finite set of pulse characteristics [33]. For a localized solution with a single maximum, these characteristics include the peak-amplitude, pulse width, center-of-mass position, and phase parameters. A full description of this method can be found in Refs. [25, 33].

In general, five moments are considered, namely, the energy \( Q\), the momentum \(P\), and higher-order generalized moments, I1, I2 and I3:

$$ \begin{gathered} Q = \mathop \int \limits_{ - \infty }^{ + \infty } \left| u \right|^{2} dT, \hfill \\ P = \frac{1}{2}\mathop \int \limits_{ - \infty }^{ + \infty } \left( {uu_{T}^{*} - u^{*} u_{T} } \right) dT, \hfill \\ I_{1} = \mathop \int \limits_{ - \infty }^{ + \infty } T\left| u \right|^{2} dT, \hfill \\ I_{2} = \mathop \int \limits_{ - \infty }^{ + \infty } \left( {T - T_{0} } \right)^{2} \left| u \right|^{2} dT, \hfill \\ I_{3} = \mathop \int \limits_{ - \infty }^{ + \infty } \left( {T - T_{0} } \right)\left( {u^{*} u_{T} - uu_{T}^{*} } \right) dT. \hfill \\ \end{gathered} $$
(2)

For Eq. (1) these five moments satisfy the following evolution equations [25, 33]:

$$ \begin{gathered} \frac{dQ}{{dZ}} = i\mathop \int \limits_{ - \infty }^{ + \infty } \left( {uR^{*} - u^{*} R} \right) dT, \hfill \\ \frac{dP}{{dZ}} = - i\mathop \int \limits_{ - \infty }^{ + \infty } \left( {u_{T} R^{*} + u_{T}^{*} R} \right) dT, \hfill \\ \frac{{dI_{1} }}{dZ} = i D P + i\mathop \int \limits_{ - \infty }^{ + \infty } T\left( {uR^{*} - u^{*} R} \right) dT, \hfill \\ \frac{{dI_{1} }}{dZ} = i D P + i\mathop \int \limits_{ - \infty }^{ + \infty } T\left( {uR^{*} - u^{*} R} \right) dT, \hfill \\ \frac{{dI_{2} }}{dZ} = - i D I_{3} + i\mathop \int \limits_{ - \infty }^{ + \infty } \left( {T - T_{0} } \right)^{2} \left( {uR^{*} - u^{*} R} \right) dT, \hfill \\ \frac{{dI_{3} }}{dZ} = 2 P \frac{{dT_{0} }}{dT} + i\mathop \int \limits_{ - \infty }^{ + \infty } \left( {2 D \left| {u_{T} } \right|^{2} - \left| u \right|^{4} } \right) dT + 2 i\mathop \int \limits_{ - \infty }^{ + \infty } \left( {T - T_{0} } \right)\left( {u_{T} R^{*} + u_{T}^{*} R} \right) dT + i\mathop \int \limits_{ - \infty }^{ + \infty } \left( {uR^{*} + u^{*} R} \right) dT, \hfill \\ \end{gathered} $$
(3)

with \(R = i{ }\delta { }u + i{ }\beta { }u_{TT} + { }i\epsilon { }\left| u \right|^{2} { }u + \left( { - \nu + i\mu } \right){ }\left| u \right|^{4} { }u\).

In the present work, the dynamical system considered was obtained using the MM for a trial function with a quartic term [28]:

$$ u\left( {T,Z} \right) = A exp\left( { - \frac{{T^{4} }}{{w^{4} }}} \right)exp\left( {icT^{2} } \right), $$
(4)

where A is the soliton amplitude, w is the soliton width, and c is the soliton chirp.

This kind of function is considered more appropriate for pulses with very high values of energy, other than higher-order Gaussian for example [25].

Since the solutions we are looking for are symmetric, with zero (transverse) velocity, some of the moments are identically zero, and the nonzero moments obtained for the quartic trial function considered above, are, respectively:

$$ Q = 1.524 A^{2} w,\quad I_{2 } = 0.239 Q w^{2} ,\quad I_{3 } = 4 i c I_{2} , $$
(5)

where Q represents the pulse energy, I2 and I3 are higher order moments [25, 28, 33].

Considering the trial function given by Eq. (4), and inserting it into the equations of the method of moments [25, 33], we obtained the following evolution equations for the soliton parameters [28]:

$$ \begin{gathered} Q_{Z} = Q\left( {2\delta - \frac{2.868 \beta }{{w^{2} }} - 1.912 c^{2} w^{2} \beta + \frac{1.103 Q\epsilon }{w} + \frac{{0.654 Q^{2} \mu }}{{w^{2} }}} \right), \hfill \\ w_{Z} = w\left( {2cD + \frac{0.388 \beta }{{w^{2} }} - 1.136 c^{2} w^{2} \beta - \frac{0.162 Q\epsilon }{w} - \frac{{0.138 Q^{2} \mu }}{{w^{2} }}} \right), \hfill \\ c_{Z} = \frac{1}{{w^{2} }}\left( { - 2c^{2} w^{2} D + \frac{3D}{{w^{2} }} - \frac{0.577 Q}{w} - \frac{{0.456 Q^{2} \nu }}{{w^{2} }} - 16.74 c\beta } \right) , \hfill \\ \end{gathered} $$
(6)

where \(Q_{Z} , w_{Z } {\text{and }}c_{z}\) denote the derivative of \(Q, w {\text{and}} c\) in order to Z.

In order to check the predictions obtained by the method of moments (MMs), Eq. (1) was numerically solved with a split—step Fourier symmetric method, in Fortran, for time steps ΔT between 0.0002 and 0.01 and distance steps ΔZ between 0.000005 and 0.0001, as described in [35]. Absorbing boundary conditions were used, as presented in [36]. An initial condition of the form 20 sech (5 T), has been considered. Some results were verified using a different computer and a different code with a Runge–Kutta 4th order pseudo-spectral method, in Matlab.

3 Results

Our present work is based on the predictions from the method of moments (MM), and on numerical simulations carried out with Eq. (1). We focus mainly in the normal dispersion regime, where VHA pulses with high energy have been found, both theoretically and experimentally [12, 28, 31, 32].

3.1 Absence of Higher-Order Reactive Nonlinearity

Figure 1 shows (a) the pulse amplitude and (b) the energy against the dispersion parameter, D, obtained from computations carried out with the dynamical system obtained with the MM (Eqs. (6)), for values of D ≤ 0, i.e., in the normal dispersion regime. The results were obtained for three different values of the nonlinear gain saturation parameter, μ, namely, μ = − 0.0001, μ = − 0.00005 and μ = − 0.00001, in the absence of higher-order reactive nonlinearity (ν = 0). The other parameter values considered in the present work, are δ = − 0.1, β = 0.2, and ε = 0.25. The net linear gain parameter, δ, should be negative in order to avoid the amplification of linear waves (background instability). On the other hand, the nonlinear gain coefficient, ε (which arises, for example, from saturable absorption), must be positive to compensate for the various types of losses in the system, namely the loss corresponding to the fiber attenuation and the loss induced by spectral filtering. Finally, the coefficient corresponding to the saturation of nonlinear gain, μ, should be negative in order to guarantee a stable soliton propagation [17, 37]. Actually, the key property of the nonlinearity in gain is to give an effective gain to the soliton and a suppression to the noise. Nonlinear optical loop mirrors (NOLM’s), nonlinear amplifying loop mirrors (NALM’s), multi-quantum-well absorbers, or nonlinear birefringent fibers with polarizers may be used as a saturable absorber for this purpose.

Fig. 1
figure 1

Pulse equilibrium amplitude (a) and energy (b), versus the dispersion parameter, D, in the normal dispersion regime. Three different values of the nonlinear gain saturation are considered, namely, μ =  − 0.00001 (thin solid curve), μ =  − 0.00005 (dotted curve), and μ =  − 0.0001 (dashed-dotted curve). The other parameter values are δ = − 0.1, β = 0.2, = ε0.25 and ν = 0. These results were obtained using Eqs. (4), based on the method of moments

It is worth to mention that the stable propagation of ultrashort pulses can be achieved even in the absence of saturation of the nonlinear gain, under the influence of some higher-order effects, namely the intra-pulse Raman scattering, the self-steepening, or the third-order dispersion effects. This fact was reported for the first time by the authors in 2005 [38] and confirmed in later works [27, 39]. Such higher-order effects can have also a significant impact on the dynamics of both stationary and pulsating ultrashort dissipative solitons [37, 40,41,42].

From Fig. 1a, it can be seen that the amplitude evolution with D is analogous for the three values of μ considered. In fact, as D changes from − 5 to 0, the pulse amplitude decreases in a nonlinear fashion, until decay is reached for D ≈ 0. The rate of amplitude decrease is slow until D ≈ − 2. For − 2 < D < 0, the amplitude suffers a fast decrease and decay is observed for D ≈ 0.

The energy evolution with D is represented in Fig. 1b. It decreases significantly with D in a nonlinear fashion, from values of 107, at D = − 5, to values 102 near D ≈ 0, vanishing at D = 0. It also can be seen that the curve concavity changes near D ≈ − 2, where the rate of the decrease of amplitude changes. Both the amplitude and the energy reach higher values for lower absolute values of μ.

To confirm the predictions based on the MM, numerical computations have been carried out with Eq. (1), assuming the same values used to obtain the Fig. 1. It was found numerically that pulse formation occurs for − 2.3 ≲ D ≲ − 0.75. For values of D > − 0.75 pulse decay was observed, while for D < − 2.3 unstable front expansion occurs. The present results are in agreement with results obtained in our previous work [28], assuming the same parameter values but = ε0.35, in wich pulse formation was observed for values of D in the interval [− 1.65, − 0.3]. Actually, assuming a lower value for the nonlinear gain (ε = 0.25), it seems that the region of existence of VHA pulses is shifted to higher absolute values of the (normal) dispersion.

Figure 2 shows (a) the VHA pulse amplitude profiles and (b) the corresponding power spectral densities found numerically for a dispersion value D = − 2.3, at which the pulses have the highest energies. In Fig. 2a, the peak amplitude values are 156.6, 70, 49.5 for the three values of nonlinear gain saturation parameter considered (μ =  − 0.00001, − 0.00005 and − 0.0001), whereas the MM predicts 151.5, 67.8, 47.9, respectively. In these circumstances the agreement can be considered very good, since the difference between the two set of results is less than 3.5% in each case.

Fig. 2
figure 2

a Amplitude profiles and b power spectral densities for D = − 2.3, and three different values of nonlinear gain saturation parameter, μ, namely, − 0.00001 (solid curves), − 0.00005 (dashed curves) and − 0.0001 (dashed-dotted curves). The other parameter values are δ = − 0.1, β = 0.2, = ε0.25 and ν = 0.

It can be seen fom Fig. 2a, that the pulses become narrower as μ → 0 − . When compared with the results obtained in Fig. 6c of [28] for D = − 1.5 and = ε0.35, we verify that the amplitude of the pulses is smaller, whereas the width is almost twice. This is due to a smaller value of the nonlinear gain parameter (ε = 0.25), and a higher absolute value of dispersion (D = − 2.3.) assumed in Fig. 2.

The pulse power spectral density versus frequency, f, \( \left( {f = \frac{\omega }{2\pi }} \right),\) is represented in Fig. 2b. As observed for the amplitude peak values, both the peak power and the spectral range increase as μ → 0 − . The pulse energies, given by \(Q = \mathop \int \limits_{ - \infty }^{ + \infty } \left| u \right|^{2} dT, \) corresponding to the three values of the nonlinear gain (μ = − 0.00001, − 0.00005, and − 0.0001), are Q = 18,873, 9550 and 7089, whereas the MM predicts 15,187, 6790 and 4800, respectively. In this case, the agreement between the two set of values is not so good. When compared with the results presented in Fig. 6d of [28] for the case D = − 1.5 and for = ε0.35, it can be seen that he spectral range is shorter, but the peak power is much higher for D = − 2.3 and = ε0.25. Actually, for each value of μ, the pulse energy is more than twice for D = − 2.3 and = ε0.25, compared with the case D = − 1.5 and = ε0.35, in spite of a smaller value of the nonlinear gain parameter.

The validity of some of the above numerical results was verified using a different code and a good agreement has been found.

3.2 Impact of the Higher-Order Reactive Nonlinearity

Since the present work is focused in the normal regime, in order to prevent pulses expansion, a positive value of ν has to be chosen [28]. The effect of the higher-order reactive nonlinearity term, ν, was characterized in Refs. [26, 33]. Specifically, it was found that the slowest energy increase in a small region of ε occurs for ν > 0. For ν > 0, the pulse becomes narrower and of high intensity. Whereas, for ν < 0, the fastest energy increase is observed, mainly due to an increase of the pulse width.

In order to study the impact of the higher-order reactive nonlinearity on the pulses formation, the dynamical system given by Eqs. (6) was numerically solved for ν = 0.00035, and for the same set of parameter values considered in Figs. 1 and 2. The value of ν = 0.00035 was chosen in order to guarantee that we observe pulse formation for μ = − 0.0001 at D = − 5. That should become clear in Sect. 3.3.

Figure 3 shows (a) the pulse equilibrium amplitude and (b) the energy versus dispersion parameter, D, given by Eqs. (6), derived using the method of moments, for values of D ≤ 0 in the presence of higher-order reactive nonlinearity, ν = 0.00035. Comparing the Figs. 1a and 3a, it is apparent that the rate of decrease of the amplitude with D is higher in the presence of higher-order reactive nonlinearity, for all of the three values of μ considered. Decay is observed in Fig. 3a for D ≈ − 0.2, while in the case for ν = 0 it was reached for D ≈ 0. It can also be seen that, the amplitude suffers a significant peak reduction in the case μ =  − 0.00001 (thin solid curve).

Fig. 3
figure 3

Pulse equilibrium amplitude (a) and energy (b), versus dispersion parameter, D. The same three different values of the nonlinear gain saturation considered in Fig. 1 are assumed, namely, μ =  − 0.00001 (thin solid curve), μ =  − 0.00005 (dotted curve), and μ =  − 0.0001 (dashed-dotted curve). The other parameter values are δ = − 0.1, β = 0.2, = ε0.25 and ν = 0.00035. These results were obtained obtained using Eqs. (4), based on the method of moments

Figure 3b shows that the energy decreases with D in a nonlinear way. However, comparing Fig. 1b and 3b it follows that the energy suffers a significant reduction in the presence of higher-order reactive nonlinearity, especially for higher absolute values of D. For example, for D = − 5 the energy is ~ 104, in contrast with 107, in the absence of higher-order reactive nonlinearity (ν = 0). On the other hand, while the amplitude achieves higher values for lower absolute values of μ, for any value of D, a different behavior is shown by the energy. Actually, the energy shows a similar behavior only for values of D between − 3 and − 0.3, approximately, but the reverse is observed for values of D ≲ − 3.

Figure 4 shows the VHA pulse amplitude profiles and the corresponding power spectral density found numerically in the presence of higher-order reactive nonlinearity (ν = 0.00035), for D = − 2.3. This is the same value of D, for wich the highest pulse energy was found, in the absence of the higher-order reactive nonlinearity effect.

Fig. 4
figure 4

a Amplitude profiles and b power spectral density for D = − 2.3, in the presence of higher-order reactive nonlinearity term (ν = 0.00035). Three different values of nonlinear gain saturation parameter μ, have been considered, namely, − 0.00001 (solid curves), − 0.00005 (dashed curves) and − 0.0001 (dashed-dotted curves). The other parameter values are δ = − 0.1, β = 0.2, and = ε0.25.

Comparing Figs. 2a and 4a, a great difference on pulse amplitudes is observed. The higher-order reactive nonlinearity reduces significantly the pulse amplitude for the lower absolute value of μ (μ =  − 0.00001), as well as the pulse width for the other two values of μ, namely, − 0.00005 and − 0.0001. As a consequence, the pulse energies are drastically reduced: from Q = 7089, 9550, 18,873 to Q = 1078, 1099, 1159, respectively. Moreover, the pulse energies become similar for the three values of μ and they are not significantly enhanced when μ → 0 − , as happens for ν = 0.

The impact on pulse energy of the higher-order reactive nonlinearity is well illustrated by Figs. 2b and 4b. Particularly for the case μ =  − 0.00001, both the peak power and the spectral range are significantly reduced. Another major difference can be observed from the power spectral density: the single peak spectra become dual peak in the presence of higher-order reactive nonlinearity. The pulses spectral range shrink, but becomes larger as μ → 0 − , as expected.

The peak amplitudes obtained numerically in Fig. 4a for the cases μ =  − 0.0001, − 0.00005, and − 0.00001 are 39.1, 47.7, and 62.2, while the MMs predicts 43.8, 57.4, and 88, respectively. On the other hand, the pulse energies obtained numerically are Q = 1078, 1099, 1159, whereas the MM predictions are Q = 2340.2, 2458.7 and 2545.4, respectively. It can be seen that the agreement is not so good as in the previous cases. In particular, the energy of each pulse obtained numerically is roughly half of the one predicted with MM.

Maybe this can be due to the choice of our trial function that seems to be more appropriated for VHA pulses with very high energy.

Assuming a value \(\nu =\) 0.00035 for the higher-order reactive nonlinearity coefficient, pulse formation was observed until D = − 5, D = − 4.6 and D = − 3.5, for μ =  − 0.0001, μ =  − 0.00005, and μ =  − 0.00001, respectively. For higher absolute values of D, expansion with instability was observed. The amplitude profiles and the corresponding power spectral densities in such limit-cases are presented in Fig. 5.

Fig. 5
figure 5

a Amplitude profiles and b power spectral density for D = − 5,(μ =  − 0.0001), D = − 4.6 (μ =  − 0.00005), and D = − 3.5, (μ =  − 0.00001), in the presence of higher-order reactive nonlinearity effect. The three different profiles correspond to three different values of the nonlinear gain saturation parameter μ, namely, − 0.00001 (solid curves), − 0.00005 (dashed curves) and − 0.0001 (dashed-dotted curves), respectively. The other parameter values are δ = − 0.1, β = 0.2, = ε0.25 and ν = 0.00035.

Two qualitatively different kinds of pulses can be observed in Fig. 5a and b. On one hand, we have a wide pulse, with a relatively high energy (10,367) and a narrow band (− 10 < f < 10) single peak spectrum for μ =  − 0.0001. On the other hand, there are two narrow pulses, with small energies (4615 and 2356) and dual peak spectra, with large spectral bands, for μ =  − 0.00005, and − 0.00001, respectively.

The numerical results obtained using a different code were found to be slightly different from those presented above, especially for the narrowest VHA pulses. However, such differences do not justify the discrepancies observed between the numerical results and the predictions based on the method of moments.

3.3 Region of Existence of the VHA Pulses

In this section we study the impact of a changing of the nonlinear gain in some of the VHA pulses found in Sect. 3.2. The set of parameter values considered are: δ = − 0.1, β = 0.2, μ =  − 0.0001 and ν = 0.00035.

Figure 6 shows the results for (a) the pulse equilibrium amplitude and (b) the energy, as a function of dispersion, D, obtained using the dynamical system of Eqs. (6), based on the MM. Different values are considered for the nonlinear gain parameter ε, namely, 0.075, 0.1, 0.2, 0.4, 1.0, and 1.5. The other parameter values are δ = − 0.1, β = 0.2, μ =  − 0.0001 and ν = 0.00035.

Fig. 6
figure 6

a Pulse equilibrium amplitude and b energy against the dispersion parameter, D, for different values of the nonlinear gain parameter, namely = ε0.075, 0.1, 0.2, 0.4, 1.0, and 1.5. The other parameter values are δ = − 0.1, β = 0.2, μ =  − 0.0001 and ν = 0.00035. These results were obtained using the dynamical system of Eqs. (4), based on the MM.

From Fig. 6a, we observe that the amplitude of the pulses evolve with D in a similar way for all the values of ε. For the higher values of the nonlinear gain (ε > 0.1), as the magnitude of D decreases, the pulse amplitude slightly decreases until D ≈ 0, where an abrup decay is observed. However, for small values of ε, (ε < 0.1), the amplitude decay is observed for – 2 ≲ D ≲ − 1. Figure 6b shows that the pulse energy decreases significantly as the magnitude of D decreases, getting values of \(\sim 10^{5}\) for D = − 5 and \(\sim 10^{2}\), for D = 0. Both the amplitude and the energy of the pulse increase for higher values of the nonlinear gain parameter.

The region of existence of VHA pulses has been investigated in the plane (ε, D). Both stationary pulses (SP) and multi-pulses (MP) have been found in this plane, as shown in Fig. 7. In all the computations an initial condition u(0,T) = 20 sech (5 T) has been considered.

Fig. 7
figure 7

Region of existence in the plane (ε, D) of VHA pulses in the presence of higher-order reactive nonlinearity effect (ν = 0.00035). The other parameter values are: δ = − 0.1, β = 0.2, and μ =  − 0.0001. The triangle and the square correspond to the solitons represented in Fig. 4 (D = − 2.3) and Fig. 5 (D = − 5) for ε = 0.25, in that order. The circles correspond to the pulses represented in Fig. 8(D = − 1), for ε = 0.4, 0.7 and 1, respectively

Above the upper ε limit of the regions shown in Fig. 7, the solutions tend to expand. Below the lower limit the solutions tend to decay. If we compare the region of existence of stationary pulses with the MM predictions, given in Fig. 6a, we can see that the upper ε limit is below the expectations. On the other hand, the lower limit agrees qualitatively with the MMs predictions as the magnitude of D decreases, especially when D approaches 0. Nevertheless, when D = 0 no decay is observed for ε > 0.4. Moreover, when ε > 0.7 multi-pulses formation is observed.

Figure 8 shows (a) the amplitude profiles and (b) the power spectral density for three stationary pulses (SP), found in the parameter plane (ε, D) at D = − 1, corresponding to the three circles, for = ε0.4 ( dashed-dotted curve), 0.7 (dashed curve) and 1.0 (solid curve). We observe that, as ε increases, the pulses become narrower, whereas their energy and spectral range increase.

Fig. 8
figure 8

a Amplitude profiles and b power spectral density for three stationary pulses (SP), found in the parameter plane (ε, D) at D = − 1, corresponding to the three circles, for = ε0.4 (dashed-dotted curve), 0.7 (dashed curve) and 1.0 (solid curve). The other parameter values are δ = − 0.1, β = 0.2, μ = − 0.0001 and ν = 0.00035

The formation of multi-pulses was also observed near D = 0, for values of the nonlinear gain > ε0.7. An example of multi-pulse formation is represented in Fig. 9 for D = 0, and = ε10. These results were obtained with the pseudo-spectral code (a 4th order code).

Fig. 9
figure 9

a Multi-pulse formation for D = 0 and = ε10. b With a proper initial condition, (150 sech(700 t)), a single pulse is obtained. The other parameter values are: δ = − 0.1, β = 0.2, μ =  − 0.0001 and ν = 0.00035

With a proper initial condition, 150 sech(700t), only one pulse was formed, as can be seen from Fig. 9b.

Figure 10a and b show the final pulse amplitude profile and its spectral density, for z = 50, respectively.

Fig. 10
figure 10

a Amplitude profile and b power spectral density for the pulse represented in Fig. 9b, for D = 0, and = ε10. The other parameter values are: δ = − 0.1, β = 0.2, μ =  − 0.0001 and ν = 0.00035

The pulse amplitude and spectral profiles are symmetric. In the temporal domain, it can be seen (Fig. 10a), that the pulse is very narrow. Its spectral range \( {\Delta } f \approx 800\).

4 Conclusions

In this work, the impact of the higher-order reactive nonlinearity on very-high amplitude (VHA) stationary solutions of the CGLE was investigated. The VHA pulses were found in previous works and are due to a singularity when the saturation of the nonlinear gain tends to zero. We considered especially the normal dispersion regime, where the VHA pulses have the highest energies.

In the present work we extended some of our previous results [28], using the method of moments (MM) and numerically solving the cubic–quintic CGLE, for a smaller value of the nonlinear gain parameter, = ε0.25. In particular we observed pulse formation until D reaches − 2.3, a higher absolute value than in our previous work (D = − 1.65, for = ε0.35). Assuming the same values for the saturation of the nonlinear gain, we also found much higher values for the pulse energies (more than twice in some cases), even if their spectral range suffer a srink. The agreement between the MM predictions and the numerical results is relatively good in this case.

In the presence of the higher-order reactive nonlinearity effect, the amplitude and the energy of the VHA pulses decrease, while their spectral range shrinks, for a given value of dispersion. Moreover, for higher absolute values of the normal dispersion, the pulse energy no longer grows significantly when μ → 0 − . This result was predicted by MM. The presence of the higher-order reactive nonlinearity effect allows the pulse formation for much higher absolute values of D. In particular, for ν = 0.00035, pulse formation was observed until D = − 5, for μ =  − 0.0001. However, since the pulses become narrower and the amplitudes are reduced, their energies are not so high as in the absence of such effect. We verified also that spectral range of each pulse still increases as μ → 0 − .

A region of existence of VHA pulses was found in the plane defined by the normal dispersion and nonlinear gain, in the presence of the higher-order reactive nonlinearity. In general, the predictions of MM for higher values of the nonlinear gain parameter are not observed and instead expansion occurs, especially for higher absolute values of D. In general, an increase of the nonlinear gain parameter, ε, has a major impact on the pulse spectral range, while the pulse energies are less affected, as D → 0 − . Multi-pulses formation was observed near D = 0 for high values of the nonlinear gain.