Abstract

In this article, we study a generalized nonlinear Schrödinger equation with second- and fourth-order diffraction, a weak nonlocality, and cubic and quintic nonlinearities. We show that this equation possesses several exact analytical solutions. In particular, it has embedded and standard bright solitons, kinks, and standard and asymmetric dark solitons. The Vakhitov-Kolokolov stability criterion shows that the embedded solitons may be VK-stable or VK-unstable, the standard solitons are VK-unstable, and the dark solitons are VK-stable. It is shown that if the embedded solitons are perturbed, they emit monochromatic radiation.

Keywords

Nonlinear Schrödinger Equation, Weak Nonlocality, Bright and Dark Solitons, Fourth-Order Diffraction, Vakhitov-Kolokolov Stability Criterion

Share and Cite:

1. Introduction

There are two variants of the spatial nonlinear Schrödinger (NLS) equation which have shown to be particularly interesting. The first one is an extension of the NLS equation which takes into account a weak nonlocality, resulting from the nonlocal dependence of the refractive index on the field intensity [1]:

$i{u}_z+{\epsilon }_1{u}_{xx}+{\gamma }_1{\left|u\right|}^{2}u+\mu \frac{{\partial }^2{\left|u\right|}^2}{\partial x^2}u=0$ (1)

where $z$ is the propagation distance, $x$ is a transverse coordinate, and the last term in this equation, originates from the approximation of the integral nonlinearity which is usually used to introduce a nonlocal contribution in the NLS equation. Generalizations of Equation (1) have also been studied in Refs. [2] [3].

A second variant of the NLS equation that we would like to mention takes into account fourth-order diffraction and a parity-time (PT) symmetric potential [4]:

$i{u}_z+{\epsilon }_1{u}_{xx}+{\epsilon }_2{u}_{4x}+\left[V_1\left(x\right)+i{V}_2\left(x\right)\right]u+{\gamma }_1{\left|u\right|}^{2}u=0,$ (2)

where $u_{xx}$ and $u_{4x}$ describe second- and fourth-order diffraction, and $V_1+i{V}_2$ is the PT-symmetric potential. A generalization of Equation (2) which includes higher-order nonlinearities, has also been studied in [5].

The observation of Equations (1) and (2) suggests that it might be interesting to investigate the consequences of taking into account simultaneously the weak nonlocality that appears in Equation (1) and the fourth-order diffraction that appears in Equation (2). Recently, Li, Ge and Shen considered this idea, and in Ref. [6] they used Anderson’s variational method [7] to obtain several approximate solutions of the equation:

$i{u}_z+{\epsilon }_1{u}_{xx}+{\epsilon }_2{u}_{4x}+{\gamma }_1{\left|u\right|}^{2}u+\mu \frac{{\partial }^2{\left|u\right|}^2}{\partial x^2}u=0$ (3)

This model is interesting, but it has an important drawback: the soliton solutions of Equation (3) do not have exact analytical expressions.

In the present communication, we show that if we generalize Equation (3) by including a quintic nonlinearity, the resulting model has different types of solitons, all of which possess exact analytical expressions. Therefore, in the following sections, we will show that the equation:

$i{u}_z+{\epsilon }_1{u}_{xx}+{\epsilon }_2{u}_{4x}+{\gamma }_1{\left|u\right|}^{2}u+{\gamma }_2{\left|u\right|}^{4}u+\mu \frac{{\partial }^2{\left|u\right|}^2}{\partial x^2}u=0$ (4)

has bright and dark solitons, kinks and asymmetric dark solitons, and we will present the exact analytical expressions of these solitons. It is worth remembering that nonlocal nonlinearities and higher-order diffraction appear in several types of nonlinear media. For example, nonlocality exists in nematic liquid crystals with molecular reorientation [8] and lead glass with heat conduction [9]. And fourth-order diffraction is taken into account in the study of photonic crystals [10], semiconductor microcavities [11] and metamaterials [12]-[15].

The structure of this article is described in the following. In Section 2, we present the analytical expression of the bright solitons of Equation (4). We will see that, depending on the values of the coefficients that appear in Equation (4), the bright solitons of this equation may be embedded solitons [16]-[25], or standard ones. Then it will be shown that the Vakhitov-Kolokolov (VK) criterion [26] [27] indicates that, depending on the coefficients of Equation (4), these embedded solitons may stable or unstable solutions. In Section 3, we present the analytical expression of the dark solitons of Equation (4), and it will be shown that the VK criterion indicates that these dark solitons are stable solutions. In Section 4, we use an algebraic method (also referred to as auxiliary equation method) [28]-[30] to obtain additional exact analytical solutions of Equation (4), in terms of the solutions of a reduced Riccati equation with constant coefficients. Using this method, we will see that Equation (4) permits the propagation of kinks and asymmetric dark solitons. Then, using a slightly different auxiliary equation, an additional solution of Equation (4) will be found. Finally, in Section 5, we present a summary and the conclusions of this work.

2. Standard and Embedded Bright Solitons

Direct substitution shows that Equation (4) has an exact solution of the form:

$u\left(x,z\right)=A\mathrm{sech}\left(\frac{x}{w}\right){\text{e}}^{ipz}$ (5)

where $w$ , $A$ and $p$ are defined by the following equations:

$\begin{array}{l}24{\epsilon }_2{\left({\gamma }_1{w}^4+4\mu w^2\right)}^2+{\gamma }_2{w}^4{\left(2{\epsilon }_1{w}^2+20{\epsilon }_2\right)}^2\\ -6\mu w^2\left(2{\epsilon }_1{w}^2+20{\epsilon }_2\right)\left({\gamma }_1{w}^4+4\mu w^2\right)=0,\end{array}$ (6)

$A^2=\frac{2{\epsilon }_1{w}^2+20{\epsilon }_2}{{\gamma }_1{w}^4+4\mu w^2} ,$ (7)

$p=\frac{{\epsilon }_1}{w^2}+\frac{{\epsilon }_2}{w^4}$ (8)

It should be noticed that (5) will be a bright soliton solution of Equation (4) only if the coefficients ${\epsilon }_1$ , ${\epsilon }_2$ , ${\gamma }_1$ , ${\gamma }_2$ and $\mu$ are such that Equation (6) has, at least, one real solution for $w$ . If Equation (6) only had complex solutions (or $w=0$ ), then Equation (4) would not have bright solitons of the form (5).

We must also observe that if we had substituted the tentative solution:

$u\left(x,z\right)=A\mathrm{sech}\left(\frac{x-vz}{w}\right){\text{e}}^{i\left(pz+qx+r{x}^2\right)}$ (9)

in Equation (4), we would have found that $v$ , $q$ and $r$ are necessarily equal to zero. Therefore, the bright solitons of Equation (4) cannot move along the x axis (which is an unexpected result).

Now let’s investigate if the bright solitons of Equation (4) are embedded solitons, or standard (i.e. not embedded) ones.

Direct substitution shows that the function:

$u\left(x,z\right)={\text{e}}^{i\left(Pz-kx\right)}$ (10)

satisfies the linear part of Equation (4) if $P$ and $k$ satisfy the linear dispersion relation:

$P={\epsilon }_2{k}^4-{\epsilon }_1{k}^2$ (11)

If ${\epsilon }_1$ and ${\epsilon }_2$ are both positive, this equation describes a curve qualitatively similar to the graph shown in Figure 1, and the minimum of the function $P\left(k\right)$ is the following value:

$P_{min}=-\frac{{\epsilon }_1^2}{4{\epsilon }_2}$ (12)

As the soliton’s wavenumber $p$ [defined in Equation (8)] is positive (if ${\epsilon }_1$ and ${\epsilon }_2$ are both positive), then we will have that:

$p=\frac{{\epsilon }_1}{w^2}+\frac{{\epsilon }_2}{w^4}>P_{min}$ (13)

and this inequality implies that the bright solitons of Equation (4) are always embedded solitons (whenever ${\epsilon }_1$ and ${\epsilon }_2$ are both positive).

Figure 1. Dispersion relation [Equation (11)] corresponding to ${\epsilon }_1=0.75$ and ${\epsilon }_2=0.1$ .

However, if ${\epsilon }_1$ and ${\epsilon }_2$ are not both positive, then, depending on the specific values of the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right\}$ , the bright soliton defined by Equations (5)-(8) might be embedded or standard. For example, if we consider the coefficients:

${\epsilon }_1=0.1$ , ${\epsilon }_2=-5$ , ${\gamma }_1=-1$ , ${\gamma }_2=1$ , $\mu =1$ (14)

then Equations (6)-(8) imply that $p=-0.033507$ . As in this case (with ${\epsilon }_2<0$ ) the dispersion relation (11) is a downward-opening quartic parabola, the negative value of $p$ is contained (embedded) in the range of wavenumbers permitted by the dispersion relation, thus implying that the soliton defined by the Equations (5)-(8) is an embedded soliton. On the other hand, if we choose the coefficients:

${\epsilon }_1=0.2$ , ${\epsilon }_2=-1$ , ${\gamma }_1=-1$ , ${\gamma }_2=1$ , $\mu =1$ (15)

the value of $p$ defined by Equation (8) is now $p=0.005600$ . As this positive value is not contained in the range of negative wavenumbers permitted by the dispersion relation (11), then, in this case, the bright soliton defined by Equations (5)-(8) is a standard (i.e., not embedded) soliton.

We have thus seen that the bright solitons of Equation (4) may be embedded or standard. Now let us analyze the stability of these solitons by means of the Vakhitov-Kolokolov (VK) criterion of stability. In order to apply this criterion, we have to determine the sign of the derivative $\text{d}N/\text{d}p$ , where $N$ is the total power defined as:

$N={\displaystyle {\int }_{-\infty }^{+\infty }{\left|u\right|}^2\text{d}x}$ (16)

Substituting the bright soliton (5) into Equation (16) we find that:

$N=2A^{2}w$ (17)

and consequently, in order to calculate the sign of the derivative $\text{d}N/\text{d}p$ , we need to express $A^2$ and $w$ as functions of $p$ .

From Equation (8), we can obtain $w$ as function of $p$ , but the precise expression of this function depends on the signs of ${\epsilon }_1$ , ${\epsilon }_2$ and $p$ . If ${\epsilon }_1$ , ${\epsilon }_2$ and $p$ are all positive, then Equation (8) implies that:

$w\left(p\right)={\left[\frac{{\epsilon }_1+\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}}{2p}\right]}^{1/2}$ (18)

and substituting (18) in (7) we obtain:

$A^2\left(p\right)=\frac{4{\epsilon }_{1}p\left({\epsilon }_1+\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}\right)+80{\epsilon }_2{p}^2}{{\gamma }_1{\left({\epsilon }_1+\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}\right)}^2+8\mu p\left({\epsilon }_1+\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}\right)}$ (19)

Then, substituting $w\left(p\right)$ and $A^2\left(p\right)$ in (17), we can obtain the function $N\left(p\right)$ .

In the case when ${\epsilon }_1>0$ , ${\epsilon }_2<0$ and $p<0$ , the function $w\left(p\right)$ has the form:

$w\left(p\right)={\left[\frac{{\epsilon }_1-\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}}{2p}\right]}^{1/2}$ (20)

and from Equations (7) and (20) we obtain:

$A^2\left(p\right)=\frac{4{\epsilon }_{1}p\left({\epsilon }_1-\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}\right)+80{\epsilon }_2{p}^2}{{\gamma }_1{\left({\epsilon }_1-\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}\right)}^2+8\mu p\left({\epsilon }_1-\sqrt{{\epsilon }_1^2+4p{\epsilon }_2}\right)}$ (21)

Then, substituting (20) and (21) in (17) we can obtain $N\left(p\right)$ .

On the other hand, when ${\epsilon }_1>0$ , ${\epsilon }_2<0$ and $p>0$ , it is not evident if $w\left(p\right)$ and $A^2\left(p\right)$ are given by the Equations (18)-(19), or the Equations (20)-(21). To determine the correct set of equations, we must observe that the Equations (6)-(8) define the values of $w$ , $A^2$ and $p$ (for a given set of coefficients ${\epsilon }_1$ , ${\epsilon }_2$ , ${\gamma }_1$ , ${\gamma }_2$ and $\mu$ ), and therefore we can calculate the positions of the points $P_1=\left(p,w\right)$ and $P_2=\left(p,A^2\right)$ . Then, we should investigate which set of equations, (18)-(19) or (20)-(21), define curves $w\left(p\right)$ and $A^2\left(p\right)$ that pass through the points $P_1$ and $P_2$ , respectively.

Let us begin by studying the stability of an embedded soliton in a case when ${\epsilon }_1$ and ${\epsilon }_2$ are both positive. In particular, let us consider the coefficients:

${\epsilon }_1=1$ , ${\epsilon }_2=1$ , ${\gamma }_1=1$ , ${\gamma }_2=-1$ , $\mu =1$ (22)

For these coefficients, the soliton’s wavenumber has the value $p=0.0893409$ . In this case, the function $N\left(p\right)$ has the form shown in Figure 2. This figure shows that $\text{d}N/\text{d}p>0$ and, according to the VK criterion, this inequality implies that the embedded soliton corresponding to the coefficients shown in (22) is a stable solution. To emphasize that the stability of this embedded soliton was determined by means of the VK criterion, we will say that this embedded soliton is VK-stable.

Now let us study the stability of an embedded soliton in a case when ${\epsilon }_1>0$ , ${\epsilon }_2<0$ and $p<0$ . If we consider the coefficients shown in (14), the functions $w\left(p\right)$ and $A^2\left(p\right)$ will be given by the Equations (20) and (21). Using these equations, we find that $N\left(p\right)$ has the form shown in Figure 3. We can see that in this case $N\left(p\right)$ has a negative slope at $p=-0.033507$ , thus implying that

Figure 2. Graph of the function $N\left(p\right)$ corresponding to the coefficients shown in (22).

Figure 3. Graph of the function $N\left(p\right)$ corresponding to the coefficients shown in (14).

the embedded soliton corresponding to the coefficients (14) is now an unstable solution. In other words, this embedded soliton is VK-unstable.

The two examples analyzed above show that the VK criterion indicates that the embedded solitons of Equation (4) may be VK-stable or VK-unstable. At first sight this result (i.e., the possibility of having VK-stable and VK-unstable embedded solitons), might seem strange. However, this result reflects the fact that the behavior of embedded solitons is rather unusual: it is known that these solitons are linearly stable, but they can be nonlinearly stable or unstable [17]. Therefore, the stability of embedded solitons is not a simple dichotomic concept that may only have two excluding possibilities: stable or unstable. Precisely because the stability of embedded solitons is not a dichotomic concept, Yang, Malomed and Kaup say that these solitons are semi-stable objects [18]. And the existence of VK-stable and VK-unstable embedded solitons is a consequence of this semi-stability.

Next let us analyze the stability of the standard soliton corresponding to the coefficients shown in (15). In this case, we use Equations (6)-(11) to determine the values: $p=0.005600$ , $w=2.45196$ and $A^2=1.45449$ . Then we can verify that the Equations (20)-(21) define curves $w\left(p\right)$ and $A^2\left(p\right)$ which pass through the points $\left(p,w\right)$ and $\left(p,A^2\right)$ , respectively. Therefore, we use the Equations (24)-(25) to determine the function $N\left(p\right)=2A^2\left(p\right) w\left(p\right)$ , and Figure 4 shows the form of this function.

Figure 4. Graph of the function $N\left(p\right)$ corresponding to the coefficients shown in (15).

Figure 4 shows that the value of the derivative $\text{d}N/\text{d}p$ (at the point $p=0.005600$ ) is negative, thus implying that the standard soliton corresponding to the coefficients shown in (15) is VK-unstable. In this way, we have seen that Equation (4) has VK-stable embedded solitons, VK-unstable embedded solitons, and VK-unstable standard solitons.

It should be noticed that the VK criterion indicates if a soliton is linearly stable, but it does not guarantee that oscillatory or nonlinear instabilities may appear. Therefore, direct numerical solutions of Equation (4) would be necessary to prove, beyond doubts, the stability of the solitons of Equation (4). Moreover, more information about the stability of these solitons might be obtained by transforming Equation (4) into a four-dimensional dynamical system, by studying the behavior of solutions of the form:

$u\left(x,z\right)=F\left(x\right){\text{e}}^{ikz}$ (23)

and defining the functions $G$ , $Q$ and $R$ as follows:

$F^{\prime }=G$ (24)

$G^{\prime }=Q$ (25)

$Q^{\prime }=R$ (26)

Substituting these definitions into Equation (4) we can see that Equation (4) transforms into:

${\epsilon }_2{R}^{\prime }=kF-{\epsilon }_{1}Q-{\gamma }_1{F}^3-{\gamma }_2{F}^5-\mu F\left(2G^2+FQ\right)$ (27)

In this way the Equation (4) has been transformed into a four-dimensional dynamical system: the system (24)-(27). As the solutions of this system define curves in a four-dimensional space, it is not trivial to imagine the forms of these curves. However, by projecting these curves on 2D planes we can obtain useful information about the behavior of these solutions. For example, if we impose the conditions $F=G=0$ , the system (24)-(27) reduces to the 2-dimensional system:

$Q^{\prime }=R\equiv h_1\left(Q,R\right)$ (28)

$R^{\prime }=-\frac{{\epsilon }_1}{{\epsilon }_2}Q\equiv h_2\left(Q,R\right)$ (29)

and in Figure 5 we can see some of the solutions of this system in the phase plane $\left(Q,R=Q^{\prime }\right)$ , for the coefficients ${\epsilon }_1=0.2$ and ${\epsilon }_2=2$ .

Figure 5. Phase plane showing solutions $\left(Q\left(x\right),R\left(x\right)\right)$ of the system (28)-(29) with ${\epsilon }_1=0.2$ and ${\epsilon }_2=2$ , and using the initial conditions: $\left(Q\left(0\right),R\left(0\right)\right)=\left(-5.5,1\right)$ [exterior curve], $\left(Q\left(0\right),R\left(0\right)\right)=\left(-5,0.75\right)$ [intermediate curve], and $\left(Q\left(0\right),R\left(0\right)\right)=\left(-4.5,0.5\right)$ [interior curve].

The Equations (28)-(29) show that $\left(Q,R\right)=\left(0,0\right)$ is an equilibrium point of the system, and Figure 5 confirms this result. Moreover, the Jacobian matrix of this system is:

$\left(\begin{array}{cc}\frac{\partial h_1}{\partial Q}& \frac{\partial h_1}{\partial R}\\ \frac{\partial h_2}{\partial Q}& \frac{\partial h_2}{RQ}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -{\epsilon }_1/{\epsilon }_2& 0\end{array}\right)$ (30)

and the determinant of this matrix, is $J\left(Q,R\right)=\left({\epsilon }_1/{\epsilon }_2\right)$ , and, as we used ${\epsilon }_1=0.2$ and ${\epsilon }_2=2$ to generate Figure 5, it follows that $J>0$ , which confirms that $\left(Q,R\right)=\left(0,0\right)$ is a center, showing neutral stability and periodic behavior. As the curves seen in Figure 5 are just 2-dimensional projections of the four-dimensional trajectories that constitute the actual solutions of the dynamical system (24)-(27), Figure 5 does not prove conclusively the existence of a stable solution of this system. However, Figure 5 strongly suggests that system (24)-(27) may indeed possess stable solutions. In this way, the phase plane analysis of the dynamical system (24)-(27) may give information about the stability of the solutions of Equation (4).

As we mentioned above, to confirm the stability predictions obtained by means of the VK criterion, requires the calculation of direct numerical solutions of Equation (4), as well as a complete stability analysis of the dynamical system (24)-(27). These calculations may be the subject of future communications, and they will be reported elsewhere.

Now, to close this section, let us study one of the most distinctive characteristics of embedded solitons: their response to perturbations. Usually, when an embedded soliton is perturbed, it resonates with the small-amplitude radiation waves capable of propagating in the system, and the perturbed soliton emits a monochromatic radiation wave with a wavenumber $k=2\pi /\lambda$ , where the value of $k$ is determined by the dispersion relation (11), where $p$ is the soliton’s wavenumber. In the following, we will see if the embedded solitons of Equation (4) behave in this way.

Let us consider again the embedded soliton corresponding to the coefficients shown in (22). With these coefficients the bright soliton of Equation (4) has the form (5), with $A=0.475929$ , $w=3.48093$ and $p=0.0893409$ . Now let us calculate numerically the solution of Equation (4) corresponding to an initial condition of the form:

$u\left(x,z=0\right)=\left(1.1\right)A\mathrm{sech}\left(\frac{x}{w}\right)$ (31)

thus implying that the pulse (31) has an amplitude 10% higher that the exact soliton.

The numerical solution of Equation (4), corresponding to the initial condition (31) and the coefficients (22), shows that the pulse emits small-amplitude radiation, as shown in Figure 6. And in Figure 7 we can see the amplitude of the Fourier transform (with respect to $x$ ) of this solution, at $z=200$ .

Figure 6. Solution $\left|u\left(x,z=200\right)\right|$ corresponding to the coefficients (22) and the initial condition (31).

Figure 7. Absolute value of the Fourier transform of $u\left(x,z=200\right)$ , whose amplitude is shown in Figure 6.

The spectrum shown in Figure 7 shows two small peaks located at the spatial frequencies $1/\lambda =\pm 0.16$ . On the other hand, if we substitute the soliton’s wavenumber $p=0.0893409$ in the dispersion relation (11), we find that $k=\pm 1.04044$ , thus implying that $1/\lambda =k/2\pi =\pm 0.165$ , thus corroborating that the small peaks seen in Figure 7 correspond to the resonance of the perturbed embedded soliton with the small-amplitude radiation waves capable of propagating in the system.

Therefore, the perturbation of the embedded solitons of Equation (4) triggers the emission of monochromatic radiation, which is the distinctive behavior of embedded solitons.

3. Dark Solitons

Now let us investigate if Equation (4) has dark soliton solutions.

Direct substitution of the function:

$u\left(x,z\right)=A\tanh \left(\frac{x}{w}\right){\text{e}}^{ipz}$ (32)

into Equation (4) shows that (32) is an exact solution of this equation if $A$ , $w$ and $p$ satisfy the following equations:

$A^2=\frac{2{\epsilon }_1{w}^2-40{\epsilon }_2}{28\mu w^2-{\gamma }_1{w}^4}$ (33)

$w^2=\left(-\frac{3\mu }{{\gamma }_2}\pm \sqrt{\frac{9{\mu }^2}{{\gamma }_2^2}-\frac{24{\epsilon }_2}{{\gamma }_2}}\right) \frac{1}{A^2}$ (34)

$p={\gamma }_1{A}^2+{\gamma }_2{A}^4+\frac{12\mu A^2}{w^2}$ (35)

and from Equations (33) and (34) we can obtain an equation which permits us to calculate the value of $w$ if we know the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,{\mu }_2\right\}$ :

$\frac{2{\epsilon }_1{w}^2-40{\epsilon }_2}{28\mu -{\gamma }_1{w}^2}=-\frac{3\mu }{{\gamma }_2}\pm \sqrt{\frac{9{\mu }^2}{{\gamma }_2^2}-\frac{24{\epsilon }_2}{{\gamma }_2}}$ (36)

Depending on the values of the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,{\mu }_2\right\}$ , the system (33-35) may have two different solutions $\left(A_i,w_i,p_i\right)$ (with $i=1,2$ and real $A_i$ , $w_i$ and $p_i$ ), it may have only one solution $\left(A,w,p\right)$ , or it may have no (real) solutions at all.

It is easy to see that if all the coefficients are positive, then Equation (34) will not give us any real and positive value for $w^2$ . Therefore, Equation (4) does not have dark solitons if all the coefficients are positive.

If ${\epsilon }_2$ and ${\gamma }_2$ are negative, and the remaining coefficients are positive, then Equation (4) might have two different dark solitons. It should be noticed, however, that the negativeness of ${\epsilon }_2$ and ${\gamma }_2$ is not enough to guarantee the existence of dark solitons, two additional conditions must be satisfied. The first condition is the following inequality:

$\frac{9{\mu }^2}{{\gamma }_2^2}-\frac{24{\epsilon }_2}{{\gamma }_2}>0$ (37)

If this condition is not satisfied, then Equation (34) will give us a complex value for $w^2$ . And the second condition is that Equation (33) must give us positive values for $A^2$ . One example of coefficients that satisfy these conditions is the following:

${\epsilon }_1=0.5,{\epsilon }_2=-1,{\gamma }_1=1,{\gamma }_2=-\frac{1}{3},\mu =1$ (38)

With these coefficients the system (33)-(35) has the following two solutions:

$A_1=0.725966271,w_1=4.771711514,p_1=0.712198685$ (39)

and:

$A_2=0.572821961,w_2=4.276179871,p_2=0.507568357$ (40)

This example thus proves that Equation (4) may indeed have two different dark solitons (for the same set of coefficients).

Now we will see that it is also possible that Equation (4) has just one dark soliton. This may occur when ${\epsilon }_2$ is negative, and the remaining coefficients are positive, because in this case, it is possible that the parenthesis that appears on the right-hand-side of Equation (34) acquires a positive value with one (and only one) of the signs in front of the square root. A set of coefficients which leads to this situation is the following:

${\epsilon }_1=1,{\epsilon }_2=-3,{\gamma }_1=1,{\gamma }_2=1,\mu =1$ (41)

With these coefficients the Equations (33)-(35) give us the following values of the soliton’s parameters:

$A=1,w=\sqrt{6},p=4$ (42)

Therefore, in this case [with the coefficients shown in (41)], Equation (4) has only one dark soliton of the form (32), with the parameters shown in (42).

Now we will calculate the stability of these dark solitons by means of the VK criterion. However, in the case of dark solitons, we cannot use the standard definition of the total energy as presented in Equation (20), because this integral has an infinite value when $u\left(x,z\right)$ is a dark soliton of the form (32). To avoid this problem, we use a renormalized total energy, defined as follows:

$N_r={\displaystyle {\int }_{-\infty }^{+\infty }\left(A^2-{\left|u\right|}^2\right)\text{d}x}$ (43)

When Equation (32) is substituted in (43), the normal result is recovered:

$N_r=2A^{2}w$ (44)

Consequently, as in the case of the bright solitons studied in Section 2, to determine the sign of the derivative $\text{d}{N}_r/\text{d}p$ , we need to express $A^2$ and $w$ as functions of the soliton’s wavenumber $p$ .

Substituting Equation (33) in (35) we obtain an equation for $w$ , and solving this equation we can find the function $w\left(p\right)$ . And once with $w\left(p\right)$ , the function $A^2\left(p\right)$ can be determined with Equation (34). This calculation, however, is not as straightforward as it seems at first sight. And the difficulty is that when we substitute Equation (33) in (35) we obtain an equation of eighth degree for $w$ , and when we solve this equation, we obtain eight functions $w_1\left(p\right)$ , $\cdots$ , $w_8\left(p\right)$ . To determine which of these eight functions is the correct one, we must observe that given the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,{\mu }_2\right\}$ , we can find the value of $w$ with Equation (36), and then the values of $A^2$ and $p$ with Equations (34) and (35). In this way we can determine the coordinates of the point $\left(p,w\right)$ , and then we can check which of the curves $w_1\left(p\right)$ , $\cdots$ , $w_8\left(p\right)$ passes through this point. We will see that only one of these eight curves passes through the point $\left(p,w\right)$ . In this way, we will determine the correct function $w\left(p\right)$ . And once having $w\left(p\right)$ , we will be able to determine $A^2\left(p\right)$ with Equation (33) or Equation (34), and then $N_r\left(p\right)$ with (44).

Following the procedure described in the previous paragraph, let us determine the VK-stability of the two dark solitons corresponding to the coefficients shown in (38). As we mentioned above, the substitution of Equation (33) in Equation (35) leads to an eighth-order algebraic equation for $w\left(p\right)$ . In the case of the coefficients (38), this equation has two complex solutions [$w_1\left(p\right)$ and $w_2\left(p\right)$ ], and six real solutions $w_3\left(p\right)$ , $\cdots$ , $w_8\left(p\right)$ . These six real solutions can be seen in Figure 8.

Figure 8. Real solutions $w_3\left(p\right)$ , $\cdots$ , $w_8\left(p\right)$ of the equation obtained by substituting (33) in (35), using the coefficients shown in (38).

The Equations (39) and (40) show that the two dark solitons corresponding to the coefficients (38) have the parameters $\left(p_1,w_1\right)=\left(0.71,4.77\right)$ and $\left(p_2,w_2\right)=\left(0.50,4.27\right)$ . And examining the six curves shown in Figure 8 we find that only the curve defined by the function $w_6\left(p\right)$ passes through the points $\left(p_1,w_1\right)$ and $\left(p_2,w_2\right)$ . Using this function, we can calculate $A^2\left(p\right)$ and $N_r\left(p\right)$ [with Equations (34) and (44)], and the form of $N_r\left(p\right)$ is shown in Figure 9. This figure shows that $\text{d}{N}_r/\text{d}p>0$ at $p_1$ and $p_2$ , and consequently the two dark solitons that exist when the coefficients take the values shown in (38) are VK-stable.

Figure 9. Graph of the function $N_r\left(p\right)$ corresponding to the coefficients shown in (38).

Now let us study the stability of the dark soliton corresponding to the coefficients shown in (41). As in the previous case [when we considered the coefficients shown in (38)], the substitution of Equation (33) in (35) leads to an eighth-order equation for $w\left(p\right)$ that has eight solutions, two of which are complex [$w_1\left(p\right)$ and $w_2\left(p\right)$ ], and the remaining six are real. The real solutions are shown in Figure 10.

As we see in (42), the dark soliton corresponding to the coefficients shown in (41) has the parameters $\left(p,w\right)=\left(4,\sqrt{6}\right)\approx \left(4,2.45\right)$ , and the only one of the six functions shown in Figure 10 that passes through this point is $w_4\left(p\right)$ . Therefore, using $w_4\left(p\right)$ , we can calculate the function $N_r\left(p\right)$ , and Figure 11 shows the form of this function.

In Figure 11, we can see that the slope of the function $N_r\left(p\right)$ [corresponding to the coefficients shown in (41)], may be positive or negative, depending on the value of $p$ . As the dark soliton corresponding to these coefficients has $p=4$ [as seen in (42)], and $\text{d}{N}_r/\text{d}p>0$ at $p=4$ , it follows that this dark soliton is a VK-stable solution.

4. Kinks and Asymmetric Dark Solitons

In the last 30 years, several techniques have been devised to obtain exact analytical solutions of nonlinear partial differential equations (NLPDEs). In many of these

Figure 10. Real solutions $w_3\left(p\right)$ , $\cdots$ , $w_8\left(p\right)$ of the equation obtained by substituting (33) in (35), using the coefficients shown in (41).

Figure 11. Graph of the function $N_r\left(p\right)$ corresponding to the coefficients shown in (41).

techniques the solution of the NLPDE is expressed as a sum of powers of a function $\phi \left(\xi \right)$ , where $\xi$ is a function of the independent variables of the NLPDE which may be used to generate a similarity reduction, and $\phi \left(\xi \right)$ is a solution of a nonlinear ordinary differential equation (NODE) that possesses exact analytical solutions. This NODE is the auxiliary equation mentioned by some authors [28]-[30]. In the following, we will use a method of this type to find additional solutions of Equation (4), different from the bright and dark solitons that we have already found in the past two sections.

In the case of Equation (4), we only have two independent variables ($x$ and $z$ ), and then we can try to find travelling wave solutions of this equation that can be expressed in the form:

$u\left(x,z\right)=U\left(\xi \right){\text{e}}^{i\theta }$ (45)

where:

$\xi =x-Vz$ (46)

$\theta =Px+Qz$ (47)

with $V$ , $P$ and $Q$ real constants, and $U\left(\xi \right)$ is a function of the form:

$U\left(\xi \right)={\displaystyle \sum }_{k=0}^n{a}_k{\phi }^k\left(\xi \right)$ (48)

where $\phi \left(\xi \right)$ is the solution of the auxiliary equation. As explained in [28], there are multiple options for choosing the auxiliary equation. Following the example of Fan [31], here we will choose as auxiliary equation a Riccati equation of the form:

${\phi }^{\prime }=b+{\phi }^2$ (49)

To find the value of the parameter $n$ [the upper limit in the sum (48)], we substitute (49) in (48), and then (45) in Equation (4). In this way Equation (4) is transformed into a long equation, and then we identify in this equation which are the terms with the highest powers of $\phi$ in the nonlinear term of highest order $\left({\left|u\right|}^{4}u=U^5{\text{e}}^{i\theta }\right)$ , and in the highest-order dispersive term $\left(u_{4x}\right)$ . These two terms turn out to be ${\phi }^{5n}$ and ${\phi }^{n+4}$ , respectively. As these two terms must cancel each other [for Equation (4) to be satisfied], it is necessary that $5n=n+4$ , thus implying that $n=1$ , and consequently Equation (48) reduces to:

$U=a_0+a_1\phi$ (50)

Now, introducing the value $n=1$ in the long equation mentioned above, we obtain that Equation (4) reduces to an equation of the form:

${\displaystyle \sum }_{k=0}^5{c}_k{\phi }^k=0$ (51)

where the coefficients $c_k$ depend on the constant $b$ which appears in (49), the constants $a_0$ and $a_1$ which appear in (48), the constants $V$ , $P$ and $Q$ which appear in (46)-(47), and the coefficients of Equation (4): ${\epsilon }_1$ , ${\epsilon }_2$ , ${\gamma }_1$ , ${\gamma }_2$ and $\mu$ . That is, the coefficients $c_k$ are functions of the form:

$c_k=c_k\left(b,a_0,a_1,V,P,Q,{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right)$ (52)

For Equation (51) to be satisfied, each of the coefficients $c_k$ must be equal to zero. At first sight it might seem that this requirement implies that we have 6 equations to be satisfied, but it is not so. The coefficients $c_0$ . $c_2$ and $c_4$ are complex expressions, and therefore their real and imaginary parts must be zero, and consequently the condition $c_k=0$ (for $k=0,\cdots ,5$ ) produces 9 equations. In particular, the conditions $Im\left(c_0\right)=0$ and $Im\left(c_4\right)=0$ imply that:

$V=0$ and $P=0$ (53)

Then, when we introduce these values in the equation $Im\left(c_2\right)=0$ , we find that this equation reduces to an identity $\left(0=0\right)$ , and therefore 6 equations remain. These equations are the following:

$a_0^2{\gamma }_1+a_0^4{\gamma }_2+2a_1^2{b}^2\mu =Q$ (54)

$2b{\epsilon }_1+16b^2{\epsilon }_2+3a_0^2{\gamma }_1+5a_0^4{\gamma }_2+4a_0^{2}b\mu +2a_1^2{b}^2\mu =Q$ (55)

$3{\gamma }_1+10a_0^2{\gamma }_2+12b\mu =0$ (56)

$2{\epsilon }_1+40b{\epsilon }_2+a_1^2{\gamma }_1+10a_0^2{a}_1^2{\gamma }_2+8a_1^{2}b\mu +4a_0^2\mu =0$ (57)

$a_1^2{\gamma }_2+2\mu =0$ (58)

$24{\epsilon }_2+a_1^4{\gamma }_2+6a_1^2\mu =0$ (59)

From (58) and (59) we can obtain the value of $a_1^2$ , and a condition that the coefficients ${\epsilon }_2$ , ${\gamma }_2$ and $\mu$ must satisfy:

$a_1^2=-\frac{2\mu }{{\gamma }_2}$ (60)

$3{\epsilon }_2{\gamma }_2={\mu }^2$ (61)

These two equations imply that the solutions that we are going to find in this section only exist when $\mu$ and ${\gamma }_2$ have opposite signs, and ${\epsilon }_2{\gamma }_2>0$ . And four equations remain to be satisfied: Equations (54)-(57). From these four equations, we should determine the values of the remaining parameters that appear in Equations (46)-(47) and (49)-(50), namely: $b$ , $a_0$ and $Q$ . We have, therefore, an overdetermined system (4 equations and three unknowns), and consequently we can’t find values of $b$ , $a_0$ and Q which satisfy (54)-(57) for arbitrary values of the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,{\mu }_2\right\}$ . In fact, we already know that we cannot choose arbitrarily the five coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right\}$ if the function $u\left(x,z\right)$ , defined by the Equations (45)-(49), is a solution of Equation (4), because the condition (61) must hold. Therefore, we are free to choose two of the coefficients $\left\{ {\epsilon }_2,{\gamma }_2,\mu \right\}$ , but the third coefficient is determined by (61). Then, to eliminate the overdetermination of the system (54)-(57), we can regard one of the two remaining coefficients $\left\{{\epsilon }_1,{\gamma }_1\right\}$ as an undetermined variable. For example, if we consider that ${\gamma }_1$ is undetermined, the system (54)-(57) will have four unknowns: $b$ , $a_0$ , Q and ${\gamma }_1$ , and the values of these variables will be found by solving the system.

For example, if we choose the coefficients:

${\epsilon }_1=0.2$ (62)

${\epsilon }_2=-0.2$ (63)

$\mu =0.23$ (64)

the values of ${\gamma }_2$ and $a_1$ are determined by (61) and (60):

${\gamma }_2=-0.088167$ (65)

$a_1=2.284160$ (66)

and solving the system (54)-(57) we find the values of the remaining four unknowns:

$b=-0.022401$ (67)

$a_0=0.448427$ (68)

$Q=0.012125$ (69)

${\gamma }_1=0.072036$ (70)

We have, therefore, the values of all the parameters that appear in the equations (45)-(47) and (49)-(50), and what we need now is the form of the solution $\phi \left(\xi \right)=\phi \left(x\right)$ of Equation (49), and when $b<0$ , as in the present case, the solution of (49) is:

$\phi \left(x\right)=-\sqrt{\left|b\right|} \frac{{\alpha }_2-{\alpha }_1{\text{e}}^{-2x\sqrt{\left|b\right|}}}{{\alpha }_2+{\alpha }_1{\text{e}}^{-2x\sqrt{\left|b\right|}}}$ (71)

where ${\alpha }_1$ and ${\alpha }_2$ are two arbitrary constants.

In the particular case when ${\alpha }_1={\alpha }_2=5$ Equations (67) and (71) define the form of $\phi \left(x\right)$ , then (66), (68), (71) and (50) give us $U\left(x\right)$ , and finally substituting $U\left(x\right)$ and the value of Q given in (69) in Equation (45), we obtain an exact solution of Equation (4). In Figure 12, we can see the shape of the function ${\left|u\right|}^2$ as a function of x.

Figure 12. Squared modulus of the solution (45) of Equation (4), corresponding to ${\epsilon }_1=0.2$ , ${\epsilon }_2=-0.2$ , ${\gamma }_1=0.072036$ , ${\gamma }_2=-0.088167$ and $\mu =0.23$ .

The solution shown in Figure 12 is an asymmetric dark soliton. It is an interesting solution, as the standard NLS equation does not have this kind of solitons. However, it is worth observing that in the case of the discrete NLS (the DNLS equation), lattice solitons similar to the solution shown in Figure 12 do indeed exist [32].

Proceeding in the same way, but using different coefficients, we can obtain additional exact analytical solutions of Equation (4). Let us generate a second solution using this method. In this case, we will consider the following coefficients:

${\epsilon }_1=5$ (72)

${\epsilon }_2=-1$ (73)

$\mu =1$ (74)

and from (60)-(61) it follows that the values of ${\gamma }_2$ and $a_1$ are:

${\gamma }_2=-1/3$ (75)

$a_1=\sqrt{6}$ (76)

Then, solving the system (54)-(57), we find the values of the remaining four unknowns:

$b=-0.0156851$ (77)

$a_0=0.94023$ (78)

$Q=0.371582$ (79)

${\gamma }_1=0.711664$ (80)

Finally, to obtain a particular solution $\phi \left(x\right)$ of Equation (49), we need to choose values for the constants ${\alpha }_1$ and ${\alpha }_2$ which appear in (71). If we choose ${\alpha }_1={\alpha }_2=5$ (as in the previous example), we find a new solution of Equation (4), whose squared modulus (as a function of x), is shown in Figure 13.

Figure 13. Squared modulus of the solution (45) of Equation (4), corresponding to ${\epsilon }_1=5$ , ${\epsilon }_2=-1$ , ${\gamma }_1=0.711664$ , ${\gamma }_2=-1/3$ and $\mu =1$ .

The solution shown in Figure 13 is a kink, but it is not the usual type of kinks that exist in the cubic-quintic NLS equation (CQNLS), as those kinks connect zero and nonzero equilibria, i.e., ${\left|u\right|}^2\to 0$ as $x\to \infty$ or $x\to -\infty$ [33]. Therefore, the kink shown in Figure 13 is an interesting solution of Equation (4).

We can obtain additional solutions of Equation (4) if we replace the auxiliary Equation (49) by the slightly more general equation:

${\phi }^{\prime }=b_0+b_1\phi +b_2{\phi }^2$ (81)

If we again look for solutions of Equation (4) in the form (45), with $\xi$ and $\theta$ given by Equations (46)-(47), and we consider that the relation between $U$ and $\phi$ is still of the form (48), we find that the equation $n=1$ still holds in this case, and therefore Equation (50) is still valid.

Proceeding as we explained at the beginning of this section, we find equations similar to Equations (51)-(52), but with different coefficients:

${\displaystyle \sum }_{k=0}^5{\alpha }_k{\phi }^k=0$ (82)

${\alpha }_k={\alpha }_k\left(a_0,a_1,b_0,b_1,b_2,V,P,Q,{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right)$ (83)

In this case five of the coefficients ${\alpha }_k$ are complex, and therefore the condition ${\alpha }_k=0$ leads to 11 real equations. However, from the equations $Im\left({\alpha }_2\right)=0$ and $Im\left({\alpha }_4\right)=0$ it follows that $P=V=0$ , and these values of $P$ and $V$ imply that $Im\left({\alpha }_k\right)$ is also zero for $k=0,1,3,5$ . In this way, only 6 equations remain, $Re\left({\alpha }_k\right)=0$ for $k=0,\cdots ,5$ , and from these 6 equations we can find the values of the parameters $a_0$ , $a_1$ , $b_0$ , $b_1$ , $b_2$ and $Q$ as functions of the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right\}$ . Then, substituting (45), (50) and (81) in Equation (4), we find that Equation (4) accepts a solution of the form:

$u\left(x,z\right)=\left[a_0+\frac{\sqrt{4b_0{b}_2-b_1^2}\tan \left(\sqrt{4b_0{b}_2-b_1^2}x/2\right)-b_1}{2b_2/a_1}\right]{\text{e}}^{iQz}$ (84)

and this function, for particular values of the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right\}$ , may reduce to a standard dark soliton of the form:

$u\left(x,z\right)=A_1{\text{e}}^{i\left(\pi /2-Q_{1}z\right)}\tanh \left(k_{1}x\right)$ (85)

where $A_1=0.344712$ , $Q_1=0.104707$ and $k_1=0.353$ .

5. Summary and Conclusions

In this article, we have studied the solutions of a generalized NLS equation [Equation (4)] that contains a weak nonlocality, second- and fourth-order diffraction terms, and cubic and quintic nonlinearities. We show that this equation has several exact analytical solutions. In particular, we find that this model has exact solutions of the following types:

- embedded solitons,

- standard (not embedded) bright solitons,

- standard dark solitons,

- asymmetric dark solitons,

- atypical kinks,

- solutions of the form (84).

Using the Vakhitov-Kolokolov criterion of stability, we found that, depending on the values of the coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right\}$ , the embedded solitons may be VK-stable or VK-unstable. The standard bright solitons that we have analyzed turned out to be VK-unstable. Moreover, using a renormalized form of the total energy, it was found that the dark solitons are VK-stable. It is worth mentioning that for a given set of coefficients $\left\{{\epsilon }_1,{\epsilon }_2,{\gamma }_1,{\gamma }_2,\mu \right\}$ , Equation (4) may have two different dark solitons, it may have only one, or none at all.

Using an auxiliary equation method (with a Riccati equation as auxiliary equation), it was found that Equation (4) has asymmetric dark solitons (see Figure 12), which is an interesting result, as asymmetric dark solitons of this type do not exist in the standard NLS equation. It should be observed that the existence of these asymmetric dark solitons is inextricably linked to the weak nonlocality present in Equation (4). This is an interesting result, as it is not obvious that a relationship should exist between nonlocality and asymmetric dark solitons.

Kinks were also obtained, but these kinks are different from the usual kinks that exist in the cubic-quintic NLS equation, as the kinks of Equation (4) do not fall to zero, neither when $x\to \infty$ , nor when $x\to -\infty$ (see Figure 13). Moreover, using Equation (81) as auxiliary equation, it was found that Equation (4) also has solutions of the form (84).

The response of embedded solitons to perturbations was also studied, and it was found that the embedded solitons of Equation (4) emit monochromatic radiation when they are perturbed, and the x-wavenumber of the radiation waves (the wavenumber along the x coordinate) coincides with the wavenumber $k$ defined by the dispersion relation (11), when the embedded soliton’s wavenumber $p$ is introduced in this relation.

The results found in this communication show that Equation (4) is an interesting model, as it possesses analytical solutions and solitons that neither the standard NLS, nor the cubic-quintic NLS equation do have.

Acknowledgements

The authors thank DGTIC-UNAM (Direc. Gral. De Cómputo y de Tecnologías de Información y Comunicación de la Universidad Nacional Autónoma de México) for their authorization to use their computer Miztli for this work, through the Project LANCAD-UNAM-DGTIC-164.

Conflicts of Interest

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

References