Chirped Embedded Solitons in a Nonlinear Schrödinger Equation with a Saturable Nonlinearity and Fourth-Order Dispersion

Áurea Espinosa-Cerón, Jorge Fujioka

Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de México, México.

DOI: 10.4236/am.2026.172006                

Abstract

In this article we study the solitary wave solutions of a generalized nonlinear Schrödinger equation which contains fourth-order dispersion and a saturable nonlinearity. We obtain both: variational solutions and direct numerical solutions. The variational method leads to an averaged Lagrangian, and Euler-Lagrange equations, which contain the dilogarithm (also known as Spence’s function), which is an interesting result from a mathematical point of view, since this special function rarely appears in the description of optical solitons. The variational solutions show that the equation studied has chirped embedded solitons, and these solitons are stable solutions. The direct numerical solutions confirm that the equation under study has chirped standard and embedded solitons, but these pulses transform into chirp-free solitons as the pulses advance along the z direction. The direct numerical solutions also show that the equation studied permits the propagation of breathers.

Keywords

Solitons, Optical Solitons, Embedded Solitons, Chirped Solitons, NLS Equation, Saturable Nonlinearity, Higher-Order Dispersion

1. Introduction

The nonlinear Schrödinger (NLS) equation is one of the most interesting nonlinear partial differential equations, as it has many interesting characteristics. Among those characteristics, the following deserve to be mentioned:

• it is an integrable equation, as its initial value problem can be analytically solved by inverse scattering [1],

• it has bright, dark and grey solitons, and also multisoliton solutions,

• it has an infinite number of conserved quantities,

• it can be obtained by the least action principle, and consequently it possesses a Lagrangian density,

• it can be obtained by means of the method of multiple scales, thus implying that the NLS equation is a universal equation,

• it can be written in Hirota’s bilinear form,

• it has the Painlevé property,

• it has Bäcklund transformations,

• it has dozens of interesting generalizations,

• it is useful to describe the propagation of waves in plasmas [2], ocean waves [3], the behavior of Bose-Einstein condensates [4], the propagation of light in photorefractive systems [5], and even the propagation of defects in DNA chains [6].

• it is the cornerstone of fiber-optic telecommunications, as this technology originated when it was discovered that light pulses (in fact, infrared pulses) could propagate long distances along optical fibers without significant distortion or dispersion, and it was shown that these pulses obey the NLS equation [7].

As mentioned above, one of the distinctive characteristics of the NLS equation is that it possesses many interesting extensions. Among the first of these extensions were equations which included higher-order dispersive terms, or higher-order nonlinearities. Higher-order dispersive terms are necessary when we want to describe the propagation of ultrashort light pulses. Mollenauer, Stolen and Gordon used pulses of 7ps in their pioneering work of 1980 [7]. However, since that time it was clear that the description of shorter pulses (with durations of 1ps, or smaller) require the addition of third- or fourth-order dispersive terms in the NLS equation. On the other hand, higher-order nonlinearities are necessary whenever we want to send pulses with powers higher than a few tenths of watts, or we want to send light pulses along semiconductor-doped optical fibers, or semiconductor-doped glasses. In those cases the usual Kerr-type nonlinearity should be replaced by a saturable nonlinearity, or by a combination of cubic and quintic nonlinearities. In particular, a paper of Hayata and Koshiba, published in the Physical Review in 1995, showed that an extension of the NLS equation containing cubic and quintic nonlinearities, had interesting soliton solutions [8]. Later, in 1997 [9], it was shown that the simultaneous inclusion of fourth-order dispersion and a quintic nonlinearity permitted the propagation of a new type of optical solitons, baptized in 1999 as embedded solitons [10].

The simultaneous inclusion of the standard cubic nonlinearity, and an additional quintic term, proved to be extraordinarily interesting. We can find thousands of papers in internet dealing with the cubic-quintic NLS equation (Google Scholar finds more than 3500 articles dealing with this equation).

It is important to observe that, in the optical context, the inclusion of a quintic nonlinearity in the NLS equation is an approximate way of taking into account a saturable nonlinearity. In other words, the cubic-quintic NLS equation:

$i{u}_z+\epsilon u_{tt}+\alpha {\left|u\right|}^{2}u+\beta {\left|u\right|}^{4}u=0$ (1)

is an approximation of an equation of the form:

$i{u}_z+\epsilon u_{tt}+\frac{{\gamma }_1{\left|u\right|}^2}{1+{\gamma }_2{\left|u\right|}^2}u=0.$ (2)

Therefore, the cubic-quintic nonlinearity:

$\alpha {\left|u\right|}^{2}u+\beta {\left|u\right|}^{4}u$ (3)

can be considered as an approximation to a nonlinearity of the form:

$\frac{{\gamma }_1{\left|u\right|}^2}{1+{\gamma }_2{\left|u\right|}^2}u$ (4)

Concerning the Equations (1) and (2), two comments might be opportune. On the one hand, it is worth observing that a generalization of Equation (1), which considers nonlinearities of the form $\alpha {\left|u\right|}^{\sigma }u+\beta {\left|u\right|}^{2\sigma }u$ , instead of the cubic and quintic nonlinearities shown in Equation (1), was thoroughly studied in Ref. [11]. In this reference we may see when this equation has soliton solutions, the exact form of these solitons, and a study of the stability of these waves. On the other hand, we should remember that the existence of soliton solutions of Equation (2) was proved in Ref. [12] as early as 1991, using a variational approach. And the approximate solitons obtained by this variational approach turned out to be chirped solitons. It should be mentioned that in recent years the interest in chirped solitons has increased significantly [13]-[18].

If we now remember that the embedded solitons (ESs) were discovered in a generalized NLS equation which combines the nonlinearity (3), with second and fourth-order dispersive terms, we may wonder if ESs might be also be found in an equation containing the nonlinearity (4), and second and fourth-order dispersion. To put it another way, we wonder if the equation:

$i{u}_z+{\epsilon }_1{u}_{tt}+{\epsilon }_2{u}_{4t}+\frac{{\gamma }_1{\left|u\right|}^2}{1+{\gamma }_2{\left|u\right|}^2}u=0$ (5)

has embedded solitons among its solutions. In this communication we will investigate this issue.

In the following section we will use Anderson’s variational method [19] to obtain approximate solutions of Equation (5). These solutions will show us if solitons are likely to exist in Equation (5). Then, in Section 3, we will obtain direct numerical solutions of Equation (5) to determine, in a more precise way, how solitary waves behave according to this equation. Some final comments will be included in Sec. 4.

2. Variational Solutions

An important characteristic of Equation (5) is that this equation can be obtained by means of the least action principle, and its Lagrangian density is the following:

$ℒ=-\mathrm{Im}\left(u^{*}{u}_z\right)-{\epsilon }_1{\left|u_t\right|}^2+{\epsilon }_2{\left|u_{tt}\right|}^2+\frac{{\gamma }_1}{{\gamma }_2}{\left|u\right|}^2-\frac{{\gamma }_1}{{\gamma }_2^2}\ln \left(1+{\gamma }_2{\left|u\right|}^2\right)$ (6)

Now we will take advantage of the existence of this Lagrangian density to obtain approximate solitary wave solutions of Equation (5). To this end we will apply Anderson’s variational method, which is frequently used in the study of optical solitons.

In order to apply Anderson’s method, we will propose a trial function of the form:

$u\left(z,t\right)=A\left(z\right)\mathrm{sech}\left[\frac{t}{w\left(z\right)}\right]{\text{e}}^{i\left[h\left(z\right)+b\left(z\right)t^2\right]}$ (7)

Substituting this function in the right-hand side of Equation (6), and calculating the integral:

$L={\displaystyle {\int }_{-\infty }^{+\infty }ℒ\text{d}t}$ (8)

we find that the Lagragian $L$ (sometimes referred to as averaged Lagrangian), has the following form:

$\begin{array}{c}L=-2h^{\prime }Aw+\frac{2{\epsilon }_2{A}^2}{w^3}+2{\gamma }_1{A}^{2}w/{\gamma }_2\\ +\frac{{\pi }^2}{6}\left(-b^{\prime }A{w}^3-4{\epsilon }_1{b}^2{A}^2{w}^3+8{\epsilon }_2{b}^2{A}^{2}w\right)+\frac{8{\epsilon }_2}{5}\frac{A^2}{w^3}\\ +\frac{14{\pi }^4}{15}{\epsilon }_2{b}^{4}A{w}^5-\frac{2}{3}\left(\frac{{\epsilon }_1}{w}+\frac{4{\epsilon }_2}{w^3}\right)A^2-\frac{{\gamma }_1}{{\gamma }_2^2}wF\left(A,{\gamma }_2\right)\end{array}$ (9)

where we have defined:

$\begin{array}{c}F\left(A,{\gamma }_2\right)=-\frac{{\pi }^2}{6}-L{i}_2\left(-1-2{\gamma }_2{A}^2-2\sqrt{{\gamma }_2{A}^2\left(1+{\gamma }_2{A}^2\right)}\right)\\ -L{i}_2\left(-1-2{\gamma }_2{A}^2+2\sqrt{{\gamma }_2{A}^2\left(1+{\gamma }_2{A}^2\right)}\right)\end{array}$ (10)

and $L{i}_n\left(z\right)$ is the polylogarithm function (also known as Jonquière’s function):

$L{i}_n\left(z\right)={\displaystyle \sum }_{k=1}^{\infty }\frac{z^k}{k^n}$ (11)

In particular, when $n=2$ , the function $L{i}_2\left(z\right)$ is referred to as dilogarithm, and also as Spence’s function. The occurrence of the dilogarithm in the averaged Lagrangian (9) is an interesting result (from a mathematical point of view), as it is extremely rare that this special function enters in the description of optical solitons. Up to the present moment, in the study of optical solitons, only two equations are known whose Lagrangians are connected with the dilogarithm: the Equation (5), and a generalized NLS equation which describes the propagation of light pulses in photorefractive crystals (the Equation (7) in Ref. [5]).

Having the Lagrangian $L$ shown in (9), we can write the Euler-Lagrange equations corresponding to the functions $A\left(z\right)$ , $w\left(z\right)$ , $h\left(z\right)$ and $b\left(z\right)$ .

From the Lagrange equation corresponding to $h\left(z\right)$ it follows that the product $A\left(z\right)w\left(z\right)$ is a constant of the motion, i.e., we have that:

$A\left(z\right)w\left(z\right)=A\left(0\right)w\left(0\right)=\text{constant}\equiv E$ (12)

and from the remaining three Lagrange equations we can obtain the following expressions for the derivatives $A^{\prime }\left(z\right)$ , $b^{\prime }\left(z\right)$ and $h^{\prime }\left(z\right)$ :

$A^{\prime }=-4{\epsilon }_{1}b{A}^2-\frac{8{\epsilon }_2}{E^2}b{A}^4+\frac{56}{5}{\pi }^2{E}^2{\epsilon }_2\frac{b^3}{A}$ (13)

$\begin{array}{c}{b}^{\prime }=-\frac{54{\epsilon }_2}{{\pi }^2{E}^6}{A}^7-\frac{6}{{\pi }^2{E}^2}\frac{{\gamma }_1}{{\gamma }_2}{A}^3-2\left({\epsilon }_1{b}^{2}A+\frac{2{\epsilon }_2}{E^2}{b}^2{A}^3\right)\\ +\frac{56}{5}{\pi }^2{\epsilon }_2{E}^2\frac{b^4}{A^2}+\frac{1}{{\pi }^2}\left(\frac{6{\epsilon }_1}{E^4}{A}^4+\frac{40{\epsilon }_2}{E^6}{A}^6\right)A\\ -\frac{3}{{\pi }^2{E}^2}\frac{{\gamma }_1}{{\gamma }_2^2}AF+\frac{3}{{\pi }^2{E}^2}\frac{{\gamma }_1}{{\gamma }_2^2}{A}^2\frac{\partial F}{\partial A}\end{array}$ (14)

$\begin{array}{c}{h}^{\prime }=\frac{81{\epsilon }_2}{10E^4}{A}^5+\frac{5{\gamma }_1}{2{\gamma }_2}A-\frac{{\pi }^2}{24}\left(12{\epsilon }_1{E}^2\frac{b^2}{A}-40{\epsilon }_2{b}^{2}A\right)\\ -\frac{7}{15}{\pi }^4{\epsilon }_2{E}^4\frac{b^4}{A^4}-\left(\frac{7{\epsilon }_1}{6E^2}{A}^2+\frac{6{\epsilon }_2}{E^4}{A}^4\right)A\\ +\frac{1}{4}\frac{{\gamma }_1}{{\gamma }_2^2}\frac{F}{A}-\frac{3}{4}\frac{{\gamma }_1}{{\gamma }_2^2}\frac{\partial F}{\partial A}\end{array}$ (15)

Solving the Equations (13)-(15) with different coefficients ${\epsilon }_n$ and ${\gamma }_n$ , and different initial conditions, we found that these equations do indeed accept soliton-like solutions. As a first example, in Figure 1 and Figure 2 we can see the amplitude $A\left(z\right)$ and chirp $b\left(z\right)$ of the solitary wave (7), for the coefficients and initial conditions shown in the figure captions. When the NLS equation is written in dimensionless form, the second-order dispersive term is usually written with a coefficient equal to 1/2 (as, for example, in [7]), or $\beta /2$ (as, for example, in [17]), because these values simplify the analytical solutions of the dimensionless NLS equation. Therefore, we used the value ${\epsilon }_1=3/2$ in the numerical examples presented in this section. The value of ${\epsilon }_2=24/49$ was chosen because this value simplifies the analytical solution of the cubic-quintic NLS equation (as shown in Ref. [9]). On the other hand, the coefficients ${\gamma }_1=1$ and ${\gamma }_2=1.3$ were chosen in order to obtain adequate damping rates in the insets shown in Figure 1 and Figure 2. However, the selection of these values is not restrictive, as the values of the coefficients $\alpha$ and $\beta$ , which appear in the cubic-quintic nonlinearity (3), can be normalized in such a way that only the signs of these coefficients are significant (see [11]). And consequently, as ${\gamma }_1=\alpha$ and ${\gamma }_2=-\beta /\alpha$ , the particular values of ${\gamma }_1$ and ${\gamma }_2$ can be changed by using a different normalization. As we can see in Figure 1 and Figure 2, $A\left(z\right)$ and $b\left(z\right)$ are practically constants, except for the extremely small oscillations which can be seen in the insets of these figures. We can see in these insets that these oscillations only occur at the beginning of the propagation of the wave along the z axis, and they decay very rapidly.

Figure 1. Amplitude $A\left(z\right)$ of the solitary wave (7), for ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=1.0735$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ . The inset shows the small oscillations that occur when the wave begins to propagate along the z axis.

Figure 2. Chirp $b\left(z\right)$ of the solitary wave (7), for ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=1.0735$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ . The inset shows the small oscillations that occur when the wave begins to propagate along the z axis.

Moreover, in Figure 3 we can see the behavior of the function $h\left(z\right)$ , for the same parameters shown in the figure captions of Figure 1 and Figure 2. Despite its apparent simplicity, this is an interesting figure. It shows that the derivative of this function, $h^{\prime }\left(z\right)$ , is a positive constant, and this constant is the wavenumber corresponding to the solitary wave (7). As the linear dispersion relation corresponding to Equation (5) is:

$k={\epsilon }_2{\omega }^4-{\epsilon }_1{\omega }^2,$ (16)

if ${\epsilon }_1$ and ${\epsilon }_2$ are both positive, the range of wavenumbers permitted by this function contains all the values of k satisfying the condition:

$k>-\frac{{\epsilon }_1^2}{4{\epsilon }_2}\equiv k_{\min }$ (17)

Consequently, the positive value of $h^{\prime }\left(z\right)$ is contained (or embedded) in this range of wavenumbers, thus implying that the soliton solution of the Equations (13)-(15) corresponding to the coefficients and initial conditions shown in the figure captions of Figures 1-3 is an embedded soliton.

Figure 3. Function $h\left(z\right)$ which appears in the solitary wave (7), for ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=1.0735$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ .

The solutions of the Equations (13)-(15) shown in Figures 1-3 clearly show that Anderson’s variational method indicates that the Equation (5) does indeed have embedded solitons. Moreover, we can use the variational method to show that these embedded solitons are stable solutions.

In Figure 4 and Figure 5 we can see the behavior of the functions $A\left(z\right)$ and $b\left(z\right)$ , respectively, when we increase in 10% the value of $A\left(0\right)$ [in comparison to the value of $A\left(0\right)$ used to obtain the solutions shown in Figures 1-3].

These graphs show that the perturbed pulse rapidly tends to a stationary solution (i.e., a soliton). And the way in which the initial perturbed pulse is attracted towards the stationary solution is better appreciated if we plot the trajectory of a point with coordinates $\left(A\left(z\right),b\left(z\right)\right)$ , as shown in Figure 6 (for $0<z<4$ ).

In a similar way, we can calculate the behavior of a perturbed soliton whose initial amplitude is 10% less that the amplitude of the soliton solution whose behavior is shown in Figures 1-3. In Figure 7 and Figure 8 we can see the behavior of the functions $A\left(z\right)$ and $b\left(z\right)$ in this case.

Moreover, in Figure 9 we can see the trajectory of a point whose coordinates are $\left(A\left(z\right),b\left(z\right)\right)$ , and $A\left(z\right)$ and $b\left(z\right)$ are the solutions of the Equations (13)-(14), where the coefficients of these equations, and the initial conditions, are the same used to obtain Figure 7, Figure 8.

Figure 4. Form of the function $A\left(z\right)$ corresponding to a perturbed soliton, whose initial amplitude is 10% higher than the amplitude of the soliton used to produce Figures 1-3. The parameters used to obtain this figure are: ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=1.18085$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ .

Figure 5. Form of the function $b\left(z\right)$ corresponding to a perturbed soliton, whose initial amplitude is 10% higher than the amplitude of the soliton used to produce Figures 1-3. The parameters used to obtain this figure are: ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=1.18085$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ .

Figure 6. Trajectory of a point with coordinates $\left(A\left(z\right),b\left(z\right)\right)$ , where $A\left(z\right)$ and $b\left(z\right)$ are solutions of the Equations (13)-(14) for $0<z<4$ . The coefficients of these equations, and the initial conditions, are the same used to obtain Figure 4 and Figure 5.

Figure 7. Form of the function $A\left(z\right)$ corresponding to a perturbed soliton, whose initial amplitude is 10% smaller than the amplitude of the soliton used to produce Figures 1-3. The parameters used to obtain this figure are: ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=0.96615$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ .

Figure 8. Form of the function $b\left(z\right)$ corresponding to a perturbed soliton, whose initial amplitude is 10% smaller than the amplitude of the soliton used to produce Figures 1-3. The parameters used to obtain this figure are: ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=1.3$ , $A\left(0\right)=0.96615$ , $b\left(0\right)=0.823068$ , $h\left(0\right)=0$ and $E=0.7$ .

Figure 9. Trajectory of a point with coordinates $\left(A\left(z\right),b\left(z\right)\right)$ , where $A\left(z\right)$ and $b\left(z\right)$ are solutions of the Equations (13)-(14) for $0<z<4$ . The coefficients of these equations, and the initial conditions, are the same used to obtain Figure 7 and Figure 8.

Figure 6 and Figure 9 clearly show that the perturbed solitons are attracted towards a stationary solution (an exact soliton) with a nonzero chirp. That the soliton whose behavior is shown in Figures 1-3 has a nonzero chirp was already evident from the horizontal line shown in Figure 3. However, when Figure 3 was obtained, it might be considered that the chirp $b\left(z\right)$ just maintained passively its initial value. On the contrary, Figure 6 and Figure 9 show that the structure of the soliton requires a well-defined nonzero value of the chirp. Therefore, this soliton solution of the Equations (13)-(15) is not only an embedded soliton, but a chirped embedded soliton, and this is an interesting result, since it is the first time that solitons of this type have been found in a single extension of the NLS equation.

Now we would like to present a second example of a solitary wave solution of the variational Equations (13)-(15). In this second example we want to show that the damped oscillations of the amplitude $A\left(z\right)$ and chirp $b\left(z\right)$ of a perturbed soliton which is approaching to an exact soliton, may be significantly different from the oscillations shown in Figure 4, Figure 5, Figure 7 and Figure 8. For example, if we consider the same coefficients used to obtain Figure 4, Figure 5, and the initial conditions $A\left(0\right)=1.35432$ , $b\left(0\right)=0.324$ , $h\left(0\right)=10.48325$ and $E=0.975977$ , the solutions $A\left(z\right)$ , $b\left(z\right)$ and $h\left(z\right)$ of the Equations (13)-(15) now have the forms shown in Figures 10-12, and a point with coordinates $\left(A\left(z\right),b\left(z\right)\right)$ now describes the trajectory shown in Figure 13.

We can see that the oscillations of the functions $A\left(z\right)$ and $b\left(z\right)$ , shown in Figure 10 and Figure 11, are appreciably different from the oscillations that we see in Figure 4, Figure 5, Figure 7 and Figure 8. However, Figure 13 shows that also in this case the point with coordinates $\left(A\left(z\right),b\left(z\right)\right)$ approaches to a stationary point in the plane $\left(A,b\right)$ , which corresponds to the shape of a soliton solution of the variational Equations (13)-(15). Moreover, the graph presented in Figure 12 shows that, after some initial small-amplitude oscillations, $h^{\prime }\left(z\right)$ [the derivative of $h\left(z\right)$ ] is always positive, thus implying that the oscillating pulse whose behavior is shown in Figures 10-13, is evolving towards a chirped embedded soliton solution of the Equations (13)-(15), as in the cases shown in Figure 6 and Figure 9.

Figure 10. Form of the function $A\left(z\right)$ corresponding to a perturbed soliton, whose initial conditions are: $A\left(0\right)=1.35432$ , $b\left(0\right)=0.324$ , $h\left(0\right)=10.48325$ and $E=0.975977$ , and the same coefficients used in Figures 1-9.

Figure 11. Form of the function $b\left(z\right)$ corresponding to a perturbed soliton, whose initial conditions are: $A\left(0\right)=1.35432$ , $b\left(0\right)=0.324$ , $h\left(0\right)=10.48325$ and $E=0.975977$ , and the same coefficients used in Figures 1-9.

Figure 12. Form of the function $h\left(z\right)$ corresponding to a perturbed soliton, whose initial conditions are: $A\left(0\right)=1.35432$ , $b\left(0\right)=0.324$ , $h\left(0\right)=10.48325$ and $E=0.975977$ , and the same coefficients used in Figures 1-9.

Figure 13. Trajectory of a point with coordinates $\left(A\left(z\right),b\left(z\right)\right)$ , where $A\left(z\right)$ and $b\left(z\right)$ are solutions of the Equations (13)-(14) for $0<z<4$ . The coefficients of these equations are the same used to obtain Figures 1-9, and the initial conditions are the same used to obtain Figures 10-12.

3. Direct Numerical Solutions

The results presented in the previous section show that the variational approximation indicates that Equation (5) is likely to have interesting soliton solutions: not just standard solitons, but embedded ones. And not the usual type of embedded solitons (ESs), found, for example in the cubic-quintic NLS equation with second- and fourth-order dispersive terms [9], or in the complex-modified KdV equation [20], but chirped embedded solitons (CHESs).

Now, in the present section, we will obtain direct numerical solutions of Equation (5) in order to verify if this equation indeed has soliton solutions, if these solitons may be embedded, and if Equation (5) permits the propagation of chirped embedded solitons.

In the following we will present numerical solutions of Equation (5) obtained with initial conditions of the form:

$u\left(z=0,t\right)=A_0\mathrm{sech}\left(\frac{t}{w_0}\right){\text{e}}^{i\left(h_0+b_0{t}^2\right)}$ (18)

As we shall see, some of these solutions are qualitatively similar to those obtained with the variational method, but some of them contain significant differences.

Let us begin by observing the solution shown in Figure 14. This solution was obtained using the coefficients ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=3.25$ , ${\gamma }_2=-2.25$ , and introducing in (18) the values $A_0=4.6272$ , $w_0=0.2/A_0$ , $h_0=0$ and $b_0=0.823068$ . This figure confirms that Equation (5) does indeed permit the propagation of solitary waves of constant height, i.e., soliton-like solutions. However, this soliton has an important difference in comparison with the solitons obtained with the variational method. The variational solitons presented in Sec. 2 were embedded solitons. On the contrary, the soliton shown in Figure 14 is not embedded, it is a standard soliton. This information can be obtained by examining the real and imaginary parts of the function $u\left(z,t=0\right)$ . These two functions show regular oscillations, and from the wavelength of these oscillations, and the phase difference between these functions, it follows that the solution of Equation (5) whose amplitude is shown in Figure 14 has the following wavenumber:

$k_s=-17$

As in this case the range of wavenumbers permitted by linear dispersion relation shown in (16) are the wavenumbers satisfying the condition:

$k>k_{\min }=-\frac{{\epsilon }_1^2}{4{\epsilon }_2}=-1.1484$ (19)

and $k_s=-17$ does not satisfy this condition, we conclude that the soliton shown in Figure 14 is not embedded. Therefore, the solution presented in Figure 14 shows that Equation (5) permits the propagation of standard solitons. Whether Equation (5) also allows the existence of embedded solitons is a question that remains to be answered.

Figure 14. Shape of the function $\left|u\left(z,t\right)\right|$ obtained from the numerical solution of Equation (5), with the coefficients ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=3.25$ , ${\gamma }_2=-2.25$ , and the initial condition (18) with the values $A_0=4.6272$ , $w_0=0.2/A_0$ , $h_0=0$ and $b_0=0.823068$ .

Now, in Figure 15, we present another solution of Equation (5) which sheds some light on the stability of the soliton solutions of this equation. This solution corresponds to the coefficients ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=-1.5$ and the parameters $A_0=1.0735$ , $w_0=0.7/A_0$ , $h_0=0$ and $b_0=0.823068$ . In this case we can see that the initial condition was not equal to an exact soliton, but the wave rapidly modifies its form, in order to acquire the shape of a true soliton. Therefore, the initial pulse, which can be considered as a perturbed soliton, is not destroyed by the perturbation, but evolves towards an exact soliton, thus implying that these solitons are stable solutions. Moreover, it is worth observing that the pulse emitted a small amplitude radiation wave during the initial transient, when the pulse was transformed into a true soliton.

An important characteristic of the pulse shown in Figure 15 is that this pulse is a standard soliton (i.e., it is not an embedded one). This result can be obtained by studying the real and imaginary parts of $u\left(z,t=0\right)$ . From the wavelength of the oscillations seen in these functions, and the phase difference between them, it follows that the wavenumber of the soliton shown in Figure 15 is $k_s=-6.1$ . As this wavenumber is out of the range of wavenumbers permitted by the linear dispersion relation of the system (i.e. the wavenumbers satisfying the condition $k>k_{\min }=-1.1484$ ), it follows that this soliton is not an embedded one.

Another interesting characteristic of the soliton shown in Figure 15 can be discovered if we plot the shape of the function $Re\left[u\left(z=0.4,t\right)\right]$ . This function can be seen in Figure 16. This figure clearly shows that the period of the oscillations of the function $Re\left[u\left(z=0.4,t\right)\right]$ grows as time advances, and this growth implies that the soliton shown in Figure 15 contains a frequency chirp. In Figure 17 we can see how the frequency of the oscillations seen in Figure 16 decreases with time. Therefore, for small values of z ($z\lesssim 0.4$ ) the soliton seen in Figure 15 is a chirped standard soliton. However, as this pulse advances along the z direction, it transforms into a chirp-free soliton, as occurs with some of the chirped solitons found in Ref. [18]. In Figure 18, for example, we can see the shape of the function $Re\left[u\left(z=3.8,t\right)\right]$ , which shows that at $z=3.8$ the chirp has disappeared. In Figure 19 we can see the frequency of the oscillations shown in Figure 18, and this figure confirms that now the frequency no longer decreases in time.

Figure 15. Shape of the function $\left|u\left(z,t\right)\right|$ obtained from the numerical solution of Equation (5), with the coefficients ${\epsilon }_1=1.5$ , ${\epsilon }_2=24/49$ , ${\gamma }_1=1$ , ${\gamma }_2=-1.5$ , and the initial condition (18) with the values $A_0=1.0735$ , $w_0=0.7/A_0$ , $h_0=0$ and $b_0=0.823068$ .

Figure 16. Shape of the function $Re\left[u\left(z=0.4,t\right)\right]$ obtained from the numerical solution of Equation (5), using the same the coefficients and the same initial condition used to obtain Figure 15.

Figure 17. Frequency of the oscillations seen in Figure 16.

Figure 18. Shape of the function $Re\left[u\left(z=3.8,t\right)\right]$ obtained from the numerical solution of Equation (5), using the same the coefficients and the same initial condition used to obtain Figure 15.

Figure 19. Frequency of the oscillations seen in Figure 18.

The encounter of the chirped soliton shown in Figure 15 reminds us that we had also found chirped solitons in Sec. 2, in that case by means of Anderson’s variational method. But there is an interesting difference. The variational chirped soliton which originated Figures 1-9, shown in Sec. 2, was an embedded soliton. On the contrary, the chirped soliton shown in Figure 15 is a standard (i.e., not embedded) one. This difference immediately reminds us of a question that had already appeared earlier in this section [in the paragraph below Equation (19)]: does Equation (5) have embedded solitons? In the following we will address this question.

In Figure 20 we can see the numerical solution of Equation (5) corresponding to the coefficients ${\epsilon }_1=1.7$ , ${\epsilon }_2=18/49$ , ${\gamma }_1=3.27$ , ${\gamma }_2=-3.46$ and the initial condition (18) with the parameters $A_0=4.6272$ , $w_0=0.99/A_0$ , $h_0=0$ and $b_0=0.823068$ This figure shows that the initial pulse was not an exact soliton, but as the pulse advances along the system, it modifies its form, until it acquires the shape of a true soliton. We can see that as the initial height of the pulse is decreasing towards a constant value, the pulse emits small amplitude radiation waves. The emission of these waves is part of the process which transforms the initial pulse into a soliton.

If we now examine the real and imaginary parts of the function $u\left(z,t=0\right)$ we will find two oscillatory functions. From the wavelength of these functions, and the phase difference between them, it follows that the function $u\left(z,t=0\right)$ has the following wavenumber:

$k_{s2}=-1.87$ (20)

As in this case ${\epsilon }_1=1.7$ and ${\epsilon }_2=18/49$ , the wavenumbers permitted by linear dispersion relation (16) are the wavenumbers which satisfy the condition:

$k>k_{\min }=-\frac{{\epsilon }_1^2}{4{\epsilon }_2}=-1.97$ (21)

As $k_{s2}=-1.87$ satisfies this condition, we conclude that in this case, the soliton shown in Figure 20 is an embedded soliton. Therefore, Equation (5) does indeed permit the propagation of embedded solitons, as the variational method predicted.

Now let’s examine another solution. In Figure 21 we can see the solution of Equation (5) which corresponds to the coefficients ${\epsilon }_1=2.1$ , ${\epsilon }_2=13/49$ , ${\gamma }_1=3.29$ , ${\gamma }_2=-3.4$ and the initial condition (18) with the parameters $A_0=4.6272$ , $w_0=0.99/A_0$ , $h_0=0$ and $b_0=0.823068$ As we can see in this figure, in this case the shape of the initial pulse is very different from an exact soliton. However, even though the initial pulse is too high in comparison to an exact soliton of a similar width, as the pulse advances along the z direction, it evolves towards an exact soliton. This process shows that the solitons of Equation (5) are stable solutions. We can also observe that as the pulse decreases its height, it emits small amplitude radiation waves, and the energy taken away by these waves probably contributes to decrease the height of the pulse.

Figure 20. Shape of the function $\left|u\left(z,t\right)\right|$ obtained from the numerical solution of Equation (5), with the coefficients ${\epsilon }_1=1.7$ , ${\epsilon }_2=18/49$ , ${\gamma }_1=3.27$ , ${\gamma }_2=-3.46$ , and the initial condition (18) with the values $A_0=4.6272$ , $w_0=0.99/A_0$ , $h_0=0$ and $b_0=0.823068$ .

Figure 21. Shape of the function $\left|u\left(z,t\right)\right|$ obtained from the numerical solution of Equation (5), with the coefficients ${\epsilon }_1=2.1$ , ${\epsilon }_2=13/49$ , ${\gamma }_1=3.29$ , ${\gamma }_2=-3.4$ , and the initial condition (18) with the values $A_0=4.6272$ , $w_0=0.99/A_0$ , $h_0=0$ and $b_0=0.823068$ .

When we calculate the wavenumber of the soliton shown in Figure 21, we obtain the value:

$k_{s3}=-3.87$ (22)

and, as this value satisfies the condition:

$k_{s3}>k_{\min }=-\frac{{\epsilon }_1^2}{4{\epsilon }_2}=-4.15$ (23)

we conclude that the soliton shown in Figure 21 is an embedded soliton.

If we now plot the shape of the function $Re\left[u\left(z=0.03,t\right)\right]$ (corresponding to the solution shown in Figure 21) we obtain the function shown in Figure 22. This figure shows that the period of the oscillations of the function $Re\left[u\left(z=0.03,t\right)\right]$ is not constant. At first (when $t<-10$ ) the period decreases as time advances, and afterwards (when $-10<t$ ) the period begins to grow again. The important point is that the frequency of these oscillations is not constant, which implies that the pulse shown in Fig. has a frequency chirp. Therefore, for small values of z (i.e., $z\lesssim 0.03$ ), the pulse shown in Figure 21 is a chirped embedded soliton. However, as the pulse advances in the z direction, it transforms into a chirp-free soliton, as it occurred with the standard soliton shown in Figure 15, and some of the chirped solitons found in Ref. [18]. Nevertheless, it is interesting that, at least within certain ranges of z, Equation (5) permits the propagation chirped standard solitons, and chirped embedded ones.

Figure 22. Shape of the function $Re\left[u\left(z=0.03,t\right)\right]$ obtained from the numerical solution of Equation (5), using the same coefficients and the same initial condition used to obtain Figure 21.

It should be noticed that Figure 21 is particularly interesting. This figure shows a pulse that evolves into a true soliton (a pulse of constant height), even though the initial condition is quite different from the shape of the final soliton. The initial pulse is much higher than the final soliton. Nevertheless, the pulse is not destroyed, on the contrary, it transforms into an exact soliton. And this process shows that the soliton presented in Figure 21 is a very stable solution. This result is particularly interesting, as the pulse shown in this figure is an embedded soliton, and this characteristic might suggest that a slight perturbation of this soliton might produce a resonance between the soliton and the small-amplitude linear waves capable of propagating in the system, and this resonance might destroy the soliton. In fact, when embedded solitons were first discovered in 1997 [9], it was thought that these solitons would always be unstable solutions due to these resonances. However, later on, in 2003, it was proved that stable embedded solitons may indeed exist, due to a precise balance between dispersion and nonlinearity [20]. Consequently, the solution of Equation (5) shown in Figure 21 shows that the combination of second- and fourth-order dispersion, and a saturable nonlinearity, permits the existence of stable embedded solitons, which is an interesting result.

And there is another issue concerning the standard soliton shown in Figure 15 and the embedded soliton shown in Figure 21, which deserves to be discussed. These two solitons have initially a frequency chirp. This is not a surprise, as the variational analysis presented in Sec. 2 indicates that the solitons of Equation (5) are chirped solitons. However, even though the variational analysis predicts that the solitons’ chirp remains constant, the direct numerical solutions of Equation (5) show that this chirp disappears as the pulses advance along the z direction. To understand this discrepancy we must remember that Anderson’s variational method usually employs chirped trial functions [as the trial function (7)], even though the true solitons that the method is trying to describe might be chirpless. In fact, Anderson’s pioneering work [19] is an example of this situation: in [19] Anderson uses a chirped trial function, while the NLS solitons (that Anderson is trying to describe), are chirpless pulses. Therefore, it shouldn’t be a surprise that the chirp which appears in our variational solutions, is not present in the true solitons of Equation (5).

It should be emphasized, however, that even though the true solitons of Equation (5) seem to be chirpless, the apparition of chirped embedded solitons (CHESs) with a constant chirp in the variational solutions of this equation is an interesting result, since CHESs have not been found in any other nonlinear partial differential equation (NLPDE). Therefore, the occurrence of CHESs in the variational analysis of Equation (5) triggers an interesting question: are there NLPDEs which have CHESs (with a constant chirp) among their true solutions? At present this question remains open.

Finally, it is worth mentioning that Equation (5) also permits the propagation of solitary waves whose heights are not constant. A solution of this type can be seen, for example, in Figure 23. This solution corresponds to the coefficients ${\epsilon }_1=1$ , ${\epsilon }_2=18/49$ , ${\gamma }_1=3.27$ , ${\gamma }_2=-3.46$ , and an initial condition of the form (18), with the parameters $A_0=4.6272$ , $w_0=0.99/A_0$ , $h_0=0$ and $b_0=0.823068$ . The solution, in this case, is a breather. It should be noticed that height of the initial pulse is much larger that the amplitude of the breather. But the initial pulse rapidly adjusts its size to the amplitude of a breather with a similar width. We should observe that in order to reduce the height of the initial pulse, a certain amount of energy must be disposed of, and the emission of the small-amplitude radiation waves which are seen in Figure 23 probably contributes to get rid of this excess energy.

Figure 23. Shape of the function $\left|u\left(z,t\right)\right|$ obtained from the numerical solution of Equation (5), with the coefficients ${\epsilon }_1=1$ , ${\epsilon }_2=18/49$ , ${\gamma }_1=3.27$ , ${\gamma }_2=-3.46$ , and the initial condition (18) with the values $A_0=4.6272$ , $w_0=0.99/A_0$ , $h_0=0$ and $b_0=0.823068$ .

4. Conclusions

In this article we investigate if a generalized NLS equation [Equation (5)] with second- and fourth-order dispersive terms, and a saturable nonlinearity, permits the propagation of embedded solitons (ESs). The possibility that Equation (5) may have ESs amongst its solutions is suggested by the fact that the saturable nonlinearity contained in Equation (5) can be approximated by the sum of a cubic and a quintic nonlinearities [shown in Equation (3)], and the generalized NLS equation containing these nonlinearities has embedded solitons.

In order to investigate if Equation (5) has ESs, we obtain both, variational and numerical solutions of this equation.

The variational approach shows that Equation (5) has chirped embedded solitons (CHESs). This is an interesting result, since this is the first time that CHESs have been found in a single extension of the NLS equation. Moreover, the variational method shows that these CHESs are stable solutions, which is an important result. In Figure 6 and Figure 9 we can see how perturbed CHESs evolve towards exact CHESs. It is also worth mentioning that the averaged Lagrangian (9), calculated in Sec. 2, and the differential Equations (14) and (15), which were obtained by the variational method, all involve the dilogarithm (also known as Spence’s function), which is an interesting result (from a mathematical point of view), as this special function rarely appears in the study of optical solitons.

On the other hand, direct numerical solutions confirmed that Equation (5) permits the propagation of chirped standard solitons, and chirped embedded ones. However, these pulses transform into chirp-free solitons (standard as well as embedded) as they advance along the z direction, as it occurs with some of the chirped solitons found in Ref. [18]. Additionally, it was also shown that Equation (5) permits the propagation of breathers.

As a final comment, we would like to emphasize that both, the variational and the direct numerical solutions of Equation (5), confirmed that Equation (5) has embedded solitons. Therefore, now we know that the embedded solitons can indeed propagate in the presence of a saturable nonlinearity.

Acknowledgements

We thank DGTIC-UNAM (Dirección General de Cómputo y de Tecnologías de Información y Comunicación de la Universidad Nacional Autónoma de México) for granting us access to the super-computer Miztli through Project LANCAD-UNAM-DGTIC-164, in order to carry out this work.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Zakharov, V.E. and Shabat, A.B. (1972) Exact Theory of Two-Dimensional Self-Focusing and One-Dimensional Self-Modulation of Waves in Nonlinear Media. Journal of Experimental and Theoretical Physics, 34, 62-69. https://api.semanticscholar.org/CorpusID:131766213
[2]	Lee, J.H., Pashaev, O.K., Rogers, C. and Schief, W.K. (2007) The Resonant Nonlinear Schrödinger Equation in Cold Plasma Physics. Application of Bäcklund-Darboux Transformations and Superposition Principles. Journal of Plasma Physics, 73, 257-272.[CrossRef]
[3]	Akhmediev, N., Ankiewicz, A. and Taki, M. (2009) Waves That Appear from Nowhere and Disappear without a Trace. Physics Letters A, 373, 675-678.[CrossRef]
[4]	Flores-Calderón, R., Fujioka, J. and Espinosa-Cerón, A. (2021) Soliton Dynamics of a High-Density Bose-Einstein Condensate Subject to a Time Varying Anharmonic Trap. Chaos, Solitons & Fractals, 143, Article ID: 110580.[CrossRef]
[5]	Betancur-Silvera, C.A., Espinosa-Cerón, A., Malomed, B.A. and Fujioka, J. (2024) Regular, Beating and Dilogarithmic Breathers in Biased Photorefractive Crystals. Axioms, 13, Article 338.[CrossRef]
[6]	Zdravković, S. and Satarić, M.V. (2008) Nonlinear Schrödinger Equation and DNA Dynamics. Physics Letters A, 373, 126-132.[CrossRef]
[7]	Mollenauer, L.F., Stolen, R.H. and Gordon, J.P. (1980) Experimental Observation of Picosecond Pulse Narrowing and Solitons in Optical Fibers. Physical Review Letters, 45, 1095-1098.[CrossRef]
[8]	Hayata, K. and Koshiba, M. (1995) Algebraic Solitary-Wave Solutions of a Nonlinear Schrödinger Equation. Physical Review E, 51, 1499-1502.[CrossRef] [PubMed]
[9]	Fujioka, J. and Espinosa, A. (1997) Soliton-Like Solution of an Extended NLS Equation Existing in Resonance with Linear Dispersive Waves. Journal of the Physical Society of Japan, 66, 2601-2607.[CrossRef]
[10]	Yang, J., Malomed, B.A. and Kaup, D.J. (1999) Embedded Solitons in Second-Harmonic-Generating Systems. Physical Review Letters, 83, 1958-1961.[CrossRef]
[11]	Micallef, R.W., Afanasjev, V.V., Kivshar, Y.S. and Love, J.D. (1996) Optical Solitons with Power-Law Asymptotics. Physical Review E, 54, 2936-2942.[CrossRef] [PubMed]
[12]	Herrmann, J. (1991) Propagation of Ultrashort Light Pulses in Fibers with Saturable Nonlinearity in the Normal-Dispersion Region. Journal of the Optical Society of America B, 8, 1507-1511.[CrossRef]
[13]	Messouber, A., Triki, H., Liu, Y., Biswas, A., Yıldırım, Y., Alghamdi, A.A., et al. (2023) Chirped Spatial Solitons on a Continuous-Wave Background in Weak Nonlocal Media with Polynomial Law of Nonlinearity. Physics Letters A, 467, Article ID: 128731.[CrossRef]
[14]	Rizvi, S.T.R., Seadawy, A.R., Mustafa, B., Ali, K. and Ashraf, R. (2022) Propagation of Chirped Periodic and Solitary Waves for the Coupled Nonlinear Schrödinger Equation in Two Core Optical Fibers with Parabolic Law with Weak Non-Local Nonlinearity. Optical and Quantum Electronics, 54, Article No. 545.[CrossRef]
[15]	Seadawy, A.R., Rizvi, S.T.R., Mustafa, B., Ali, K. and Althubiti, S. (2022) Chirped Periodic Waves for an Cubic-Quintic Nonlinear Schrödinger Equation with Self Steepening and Higher Order Nonlinearities. Chaos, Solitons & Fractals, 156, Article ID: 111804.[CrossRef]
[16]	Daoui, A.K., Messouber, A., Triki, H., Zhou, Q., Biswas, A., Liu, W., et al. (2021) Propagation of Chirped Periodic and Localized Waves with Higher-Order Effects through Optical Fibers. Chaos, Solitons & Fractals, 146, Article ID: 110873.[CrossRef]
[17]	Triki, H. and Kruglov, V.I. (2021) Chirped Periodic and Localized Waves in a Weakly Nonlocal Media with Cubic-Quintic Nonlinearity. Chaos, Solitons & Fractals, 153, Article ID: 111496.[CrossRef]
[18]	Tedjani, A.H., Seadawy, A.R., Rizvi, S.T.R. and Solouma, E. (2023) Construction of Chirped Propagation with Jacobi Elliptic Functions for the Nonlinear Schrödinger Equations with Quadratic Nonlinearity with Inter-Modal and Spatio-Temporal Dispersions. The European Physical Journal Plus, 138, Article No. 973.[CrossRef]
[19]	Anderson, D. (1983) Variational Approach to Nonlinear Pulse Propagation in Optical Fibers. Physical Review A, 27, 3135-3145.[CrossRef]
[20]	Rodríguez, R.F., Reyes, J.A., Espinosa-Cerón, A., Fujioka, J. and Malomed, B.A. (2003) Standard and Embedded Solitons in Nematic Optical Fibers. Physical Review E, 68, Article ID: 036606.[CrossRef] [PubMed]