Dipole Solitons in a Nonlocal Nonlinear Medium with Self-Focusing and Self-Defocusing Quintic Nonlinear Responses

Abstract

Stability dynamics of dipole solitons have been numerically investigated in a nonlocal nonlinear medium with self-focusing and self-defocusing quintic nonlinearity by the squared-operator method. It has been demonstrated that solitons can stay nonlinearly stable for a wide range of each parameter, and two nonlinearly stable regions have been found for dipole solitons in the gap domain. Moreover, it has been observed that instability of dipole solitons can be improved or suppressed by modification of the potential depth and strong anisotropy coefficient.

1. Introduction

Many phenomena in nature are modeled mathematically using nonlinear differential equations. Traveling wave solutions of nonlinear partial differential equations play a significant role in nonlinear wave propagation problems that are observed in various fields such as nonlinear optics, fluid dynamics, plasma physics, elastic media, and biology [1]. Some of the solutions to such nonlinear wave propagation problems are called solitons, which are localized wave solutions.

Optical solitons are formed because of the balance between the medium’s diffraction and the self-phase modulation [2]. As a consequence of this, an optical field that does not change its shape occurs during propagation [3]. Recently, spatial solitons that can be used for optical switching and processing applications [4] have been extensively investigated in nonlinear optical systems with external optical lattices. There is a considerable amount of research about this subject in the literature. In 2003, Segev et al. experimentally observed spatial solitons in optically induced periodic potentials [5]. Fundamental and vortex solitons with real or complex lattices have been investigated in optical media with the cubic Kerr-type [6, 7, 8, 9, 10], the saturable [11], and competing nonlinearities [12]. Moreover, the existence of solitons has been observed in aperiodic or quasicrystal lattice structures [13, 14, 15, 16, 17, 18] and the lattices that possess defects [19, 20] and dislocations [21, 22].

The dynamics of solitons are governed by nonlinear Schrödinger (NLS) type equations in optical media with nonlinearities and/or external potentials as in the referred studies. Additionally, the cubic nonlinear NLS equation needs to be modified to describe nonlinear optical materials that have both cubic and quadratic nonlinear responses [23, 24, 25, 26, 27, 28, 29], such as potassium niobate (${\text{KNbO}}_3$) [30] or lithium niobate (${\text{LiNbO}}_3$) [31]. These dynamics in quadratically polarized media are governed by the NLS equation with coupling to a mean term (d.c. field), which are denoted as NLSM systems and sometimes referred to as Benney-Roskes or Davey-Stewartson systems [32, 33].

NLSM equations were first studied by Benney and Roskes for water of finite depth in the free surface conditions in 1969 [32]. Later, in 1974, Davey and Stewartson derived the limiting integrable case, which is a reduced case of the Benney–Roske’s system by studying the evolution of a 3D wave packet for water of finite depth [33]. In 1975, Ablowitz and Haberman [34] studied the integrability of NLSM systems in the shallow water limit. The effects of surface tension were included in the results of Benney and Roskes by Djordevic and Reddekopp [35] in 1977. From the first principles, Ablowitz et al. [23, 36, 37] discovered that NLSM-type equations describe the evolution of the electromagnetic field in a quadratic nonlinear media. The general NLSM system is given by [23, 36, 37]

$${\mathit{iu}}_z+\mathrm{\Delta }u+{\left|u\right|}^{2}u-\mathit{\rho u\varphi }=0,{\varphi }_{\mathit{xx}}+v{\varphi }_{\mathit{yy}}={\left({\left|u\right|}^2\right)}_{\mathit{xx}},$$E1

where u(x, y, z) corresponds to the normalized amplitude of the envelope of the static electric field propagating in the z direction, x and y are transverse spatial coordinates. $\mathrm{\Delta }u\equiv u_{\mathit{xx}}+u_{\mathit{yy}}$ corresponds to diffraction, the cubic term in u originates from the Kerr-type nonlinear change of the refractive index. The parameter ρ is a coupling constant that comes from the combined optical rectification and electro-optic effects modeled by the $\varphi \left(x,y\right)$ field, and v is the coefficient that comes from the anisotropy of the material [37]. Such systems of equations arise due to the growth and depletion of the fundamental and second-harmonic fields at the moment that the phase velocity of the fundamental and the second-harmonic wave are not equal during propagation [38]. When the phase-matching condition is not satisfied, the equation of the second-harmonic field can be solved directly and generates an additional self-phase modulation contribution as a result of cascaded nonlinearity. Similarly, the NLSM systems describe the nonlocal–nonlinear coupling between the first harmonic with the cascading effect from the second harmonic and a static field that is related to the mean term [36, 37].

Wave collapses play a significant role in various branches of science. The peak amplitude of the wave solutions tends to infinity (blow-up) in finite time or finite propagation distance when a singularity occurs. This phenomenon is often called wave collapse [39]. In the NLS equation, it was first observed numerically by Kelley in 1965 [40]. In fact, this wave collapse phenomenon is similar for the NLSM systems. Wave collapse in the NLSM systems occurs with a modulated profile [41]. Merle and Raphael [42] analyzed the collapse behavior of the NLS equation and other related equations in detail. Moreover, Moll et al. investigated experimental observations of optical wave collapse in cubic nonlinearity and showed that the amplitude of the wave increases as the spatial extent decreases in a self-similar profile [43]. In Ref [39], Ablowitz and coworkers studied wave collapse that occurs with a quasi-self-similar profile in the NLSM system and found that collapse can be arrested by the small nonlinear saturation. Furthermore, in Ref [30], NLSM collapse was arrested by wave self-rectification. In this aforementioned study, they considered only the nonlinear evolution of beams with an initial Gaussian beam profile with several values of input power and/or beam ellipticity and found that the wave collapse can be arrested by increasing the coupling constant ρ or for an initially highly elliptic beam. Recently, the NLSM system collapse was arrested by adding a real periodic [24] and partially parity-time-symmetric [44] and azimuthal [45] external lattices (potential) to the governing system, and it was shown numerically that modification of potential depth provides great controllability on the stability of soliton.

More recently, Bağcı et al. [46] have numerically investigated stability dynamics of fundamental lattice solitons that are solutions of extended NLSM system in a nonlocal nonlinear medium with self-focusing and self-defocusing quintic nonlinearity. It has been shown that as the absolute value of γ increases for both self-focusing and self-defocusing cases, the obtained fundamental solitons become nonlinearly unstable. However, the stability of unstable fundamental solitons can be improved by modification of potential depth [46].

Dipole (two-phased) and higher-phase vortex solitons in the presence of an induced lattice have been studied analytically and experimentally in Bose-Einstein condensates (BECs) [47, 48] and in optical Kerr media [49, 50, 51, 52, 53, 54]. In recent years, these types of solitons have attracted considerable interest because of their unique features and potential applications [55].

In this chapter, we numerically study the existence and stability of dipole soliton solutions of the NLSM system in a nonlocal nonlinear medium with the self-defocusing quintic nonlinear response by adding an external lattice. In fact, this study is about the dynamics of dipole solitons instead of fundamental solitons in the problem that Bağcı and coworkers have addressed in recent book chapter [46]. The purpose of this study is to numerically investigate the effects of the strength of quintic nonlinearity that specify characteristics of the model and variation of potential depth on the existence and stability of dipole solitons. In several applications, many optical materials such as chalcogenide glasses are required quintic and seventh-order effects in addition to cubic nonlinear effects [56], and effective higher-order nonlinearities can reveal with pure Kerr materials in an inhomogeneous propagation media [57, 58, 59].

The chapter is outlined as follows: In Sec. 2, we present the model equations, and the squared-operator method is explained so that it is modified for the model. The dipole solitons are computed by this numerical method. Nonlinear evolution of the dipole solitons is examined to perform stability analysis, In Sec. 3. Finally in Sec. 4, results of this study are outlined.

2. The model

In this chapter, we modify the NLSM system (1) as follows to describe the dynamics of lattice solitons in a nonlocal nonlinear medium with cubic and quintic nonlinearity

$$\begin{array}{l}{\mathit{iu}}_z+\frac{1}{2}\mathrm{\Delta }u+\beta {\left|u\right|}^{2}u-\mathit{\rho u\varphi }+\gamma {\left|u\right|}^{4}u-V\left(x,y\right)u=0,\\ {\varphi }_{\mathit{xx}}+v{\varphi }_{\mathit{yy}}={\left({\left|u\right|}^2\right)}_{\mathit{xx}}\end{array}$$E2

where γ is the coefficient of quintic nonlinearity and V(x, y) is the optical lattice. In this chapter, we consider lattices that can be written as the intensity of a sum of N phase-modulated plane waves [13]

$$V\left(x,y\right)=\frac{V_0}{N^2}{\left|\sum _{n=0}^{N-1}{e}^{i\left(k_x^{n}x+k_y^{n}y\right)}\right|}^2,$$E3

where V₀> 0 is the peak depth of the potential and the wave vector $\left(k_x^n,k_y^n\right)=\left[Kcos\left(2\mathit{\pi n}/N\right),Ksin\left(2\mathit{\pi n}/N\right)\right].$ The potential for N = 2, 3, 4, 6 yield crystal (periodic) lattices, while N = 5, 7 yield quasi-crystals (aperiodic) lattices. Contour image, contour plot, and diagonal cross-section of the lattice V(x, y) are plotted in Figure 1 for $V_0=12.5,$N = 4 and $k_x=k_y=2\pi .$ It can be seen that the lattice is periodic, and the center of lattice is a local maximum.

Figure 1.

(a) Contour image, (b) contour plot, and (c) diagonal cross-section of the lattice V(x, y) when $V_0=12.5,$ N = 4 and (x,y) ∈ [15,15].

2.1 Numerical solution for the dipole solitons

Yang and Lakoba developed an iterative numerical method called the squared-operator method (SOM) [60]. The idea of this method is to iterate a modified differential equation whose linearization operator is square of the original equation together with a preconditioning (or acceleration) operator. To obtain the soliton solution of the (2+1)D NLSM model, this method is modified as follows:

Soliton solutions are sought in the form $u\left(x,y,z\right)=U\left(x,y\right)e^{\mathit{i\mu z}}$ where U(x, y) is real-valued function and μ is the propagation constant (or eigenvalue). Substituting the ansatz u(x, y, z), we get the following expressions:

$$\begin{array}{l}{u}_z=\mathit{i\mu }{\mathit{Ue}}^{\mathit{i\mu z}},\\ u_{\mathit{xx}}=U_{\mathit{xx}}{e}^{\mathit{i\mu z}},\\ u_{\mathit{yy}}=U_{\mathit{yy}}{e}^{\mathit{i\mu z}},\\ {\left|u\right|}^2={\mathit{Ue}}^{\mathit{i\mu z}}{U}^{\ast }{e}^{-\mathit{i\mu z}}={\left|U\right|}^2,\\ {\left|u\right|}^4={\left|U\right|}^2{\left|U\right|}^2={\left|U\right|}^4\end{array}$$E4

where $U=U^{\ast }$ in our case. Substituting the set of the terms in Eq. (4) into the (2+1) NLSM model, the following nonlinear equations for U are obtained

$$\begin{array}{l}-\mathit{\mu U}+\frac{1}{2}\mathrm{\Delta }U+\beta {\left|U\right|}^{2}U-\mathit{\rho \varphi U}+\gamma {\left|U\right|}^{4}U-\mathit{VU}=0,\\ {\varphi }_{\mathit{xx}}+v{\varphi }_{\mathit{yy}}={\left({\left|U\right|}^2\right)}_{\mathit{xx}}.\end{array}$$E5

Applying the Fourier transform to the eigenequations system (5) yields

$$\begin{array}{l}-\mu \widehat{U}-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{U}+\mathcal{F}\left\{\beta {\left|U\right|}^{2}U-\mathit{\rho \varphi U}+\gamma {\left|U\right|}^{4}U-\mathit{VU}\right\}=0,\\ k_x^2\widehat{\varphi }+{\mathit{vk}}_y^2\widehat{\varphi }=k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\},\end{array}$$E6

where $\mathcal{F}$ denotes the Fourier transform, $\widehat{U}=\mathcal{F}\left\{U\right\},k_x$ and $k_y$ are the Fourier transform variables. Isolating $\widehat{\varphi }$ from the second equation of Eq. (6) gives

$$\widehat{\varphi }=\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}.$$E7

Taking the inverse Fourier transform of Eq. (7), we get

$$\varphi ={\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\},$$E8

where ${\mathcal{F}}^{-1}$ denotes the inverse Fourier transform and during iteration, the first element of $k_x^2+{\mathit{vk}}_y^2$ is set to 1 in order to avoid division by zero error. By applying the inverse Fourier transform to first equation of Eq. (6) and substituting Eq. (8) into the obtained equation, we get

$$-\mathit{\mu U}+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{U}\right\}+\beta {\left|U\right|}^{2}U-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U+\gamma {\left|U\right|}^{4}U-\mathit{VU}=0.$$E9

To obtain operator $L_0$, Eq. (9) can be written as

$$L_{0}U={\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{U}\right\}+T_{0}U=0,$$E10

where

$$T_0=-\mu +\beta {\left|U\right|}^2-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}+\gamma {\left|U\right|}^4-V.$$E11

Now, we should obtain operator $L_1$, which denotes the linearized operator of $L_{0}U=0$, with respect to the solution $\mathit{U}$, i.e., $L_0\left(U+\stackrel{\sim }{U}\right)=L_1\stackrel{\sim }{U}+O\left({\stackrel{\sim }{U}}^2\right)$, where $\stackrel{\sim }{U}\ll 1$. However, it should be noted that we have obtained the operator $L_0$ by substituting the mean field term $\varphi \left(x,y\right)$ into the governing equation. Therefore, at this point, we have to perturb $\varphi \left(x,y\right)$ function as well. In accordance with this purpose, the soliton solution and the mean-field term should be perturbed as follows, respectively,

$$\begin{array}{l}u\left(x,y,z\right)=\left[U\left(x,y\right)+\stackrel{\sim }{U}\left(x,y\right)\right]e^{\mathit{i\mu z}}\\ \varphi \left(x,y\right)=\varphi \left(x,y\right)+\stackrel{\sim }{\varphi }\left(x,y\right)\end{array},$$E12

where $\stackrel{\sim }{U}\ll 1$ and $\stackrel{\sim }{\varphi }\ll 1$. Firstly, consider that Eq. (5) can be written as a general type of nonlinearities

$$\begin{array}{l}-\mathit{\mu U}+\frac{1}{2}\mathrm{\Delta }U+F\left({\left|U\right|}^2\right)U-\mathit{\rho \varphi U}-\mathit{VU}=0,\\ {\varphi }_{\mathit{xx}}+v{\varphi }_{\mathit{yy}}={\left({\left|U\right|}^2\right)}_{\mathit{xx}},\end{array}$$E13

where $F\left({\left|U\right|}^2\right)=\beta {\left|U\right|}^2+\gamma {\left|U\right|}^4$ for cubic-quintic nonlinearity. Then, substituting perturbation $U+\stackrel{\sim }{U}$ and using linear Taylor expansion yield

$$\begin{array}{l}F\left({\left|U\right|}^2\right)=F\left({\left|U+\stackrel{\sim }{U}\right|}^2\right)=F\left(U^2+2U\stackrel{\sim }{U}+{\stackrel{\sim }{U}}^2\right)\\ \approx F\left(U^2+2U\stackrel{\sim }{U}\right)\\ \approx F\left({\left|U\right|}^2\right)+2U\stackrel{\sim }{U}{F}_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)\end{array},$$E14

where $U=U^{\ast },\stackrel{\sim }{U}={\stackrel{\sim }{U}}^{\ast }$in our case and $F_{{\left|U\right|}^2}^{\prime }=\mathit{\partial F}/\mathit{\partial }{\left|U\right|}^2$

$$\begin{array}{l}F\left({\left|U\right|}^2\right)U=F\left({\left|U+\stackrel{\sim }{U}\right|}^2\right)\left(U+\stackrel{\sim }{U}\right)\\ \approx \left[F\left({\left|U\right|}^2\right)+2U\stackrel{\sim }{U}{F}_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)\right]\left(U+\stackrel{\sim }{U}\right)\\ \approx \left(U+\stackrel{\sim }{U}\right)F\left({\left|U\right|}^2\right)+2U^2\stackrel{\sim }{U}{F}_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)+2U{\stackrel{\sim }{U}}^2{F}_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)\\ \approx \left(U+\stackrel{\sim }{U}\right)F\left({\left|U\right|}^2\right)+2U^2\stackrel{\sim }{U}{F}_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)+O\left({\stackrel{\sim }{U}}^2\right).\end{array}$$E15

Substituting perturbations in Eq. (12) and Eq. (15) into Eq. (13) and only terms of $O\left(\stackrel{\sim }{U}\right)$ and $O\left(\stackrel{\sim }{\varphi }\right)$ are retained, we get

$$\begin{array}{l}-\mu \left(U+\stackrel{\sim }{U}\right)+\frac{1}{2}\mathrm{\Delta }\left(U+\stackrel{\sim }{U}\right)+\left(U+\stackrel{\sim }{U}\right)F\left({\left|U\right|}^2\right)+2U^2\stackrel{\sim }{U}{F}_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)\\ -\mathit{\rho \varphi }\left(U+\stackrel{\sim }{U}\right)-\rho \stackrel{\sim }{\varphi }U-V\left(U+\stackrel{\sim }{U}\right)=0,\\ {\varphi }_{\mathit{xx}}+v{\varphi }_{\mathit{yy}}+{\stackrel{\sim }{\varphi }}_{\mathit{xx}}+v{\stackrel{\sim }{\varphi }}_{\mathit{yy}}={\left({\left|U\right|}^2\right)}_{\mathit{xx}}+{\left(2U\stackrel{\sim }{U}\right)}_{\mathit{xx}}.\end{array}$$E16

Substituting $F\left({\left|U\right|}^2\right)=\beta {\left|U\right|}^2+\gamma {\left|U\right|}^4\Rightarrow F_{{\left|U\right|}^2}^{\prime }\left({\left|U\right|}^2\right)=\beta +2\gamma {\left|U\right|}^2$ into Eq. (16) yields

$$\begin{array}{l}-\mu \left(U+\stackrel{\sim }{U}\right)+\frac{1}{2}\mathrm{\Delta }\left(U+\stackrel{\sim }{U}\right)+\left(U+\stackrel{\sim }{U}\right)\left(\beta {\left|U\right|}^2+\gamma {\left|U\right|}^4\right)+2U^2\stackrel{\sim }{U}\left(\beta +2\gamma {\left|U\right|}^2\right)\\ -\mathit{\rho \varphi }\left(U+\stackrel{\sim }{U}\right)-\rho \stackrel{\sim }{\varphi }U-V\left(U+\stackrel{\sim }{U}\right)=0,\\ {\varphi }_{\mathit{xx}}+v{\varphi }_{\mathit{yy}}+{\stackrel{\sim }{\varphi }}_{\mathit{xx}}+v{\stackrel{\sim }{\varphi }}_{\mathit{yy}}={\left({\left|U\right|}^2\right)}_{\mathit{xx}}+{\left(2U\stackrel{\sim }{U}\right)}_{\mathit{xx}}.\end{array}$$E17

Applying the Fourier transform to the first perturbed equation in Eq. (17) and the inverse Fourier transform to the obtained equation, we get

$$\begin{array}{l}-\mu \left(U+\stackrel{\sim }{U}\right)+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\left(\widehat{U}+\widehat{\stackrel{\sim }{U}}\right)\right\}+\left(U+\stackrel{\sim }{U}\right)\left(\beta {\left|U\right|}^2+\gamma {\left|U\right|}^4\right)\\ +2U^2\stackrel{\sim }{U}\left(\beta +2\gamma {\left|U\right|}^2\right)-\mathit{\rho \varphi }\left(U+\stackrel{\sim }{U}\right)-\rho \stackrel{\sim }{\varphi }U-V\left(U+\stackrel{\sim }{U}\right)=0.\end{array}$$E18

Using Eq. (5) and the second equation in Eq. (17), following equation is obtained

$${\stackrel{\sim }{\varphi }}_{\mathit{xx}}+v{\stackrel{\sim }{\varphi }}_{\mathit{yy}}={\left(2U\stackrel{\sim }{U}\right)}_{\mathit{xx}}.$$E19

Applying Fourier transform to Eq. (19) and isolating $\widehat{\stackrel{\sim }{\varphi }}$ from obtained equation, we get

$$\begin{array}{l}{k}_x^2\widehat{\stackrel{\sim }{\varphi }}+{\mathit{vk}}_y^2\widehat{\stackrel{\sim }{\varphi }}=k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}\\ \widehat{\stackrel{\sim }{\varphi }}=\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}.\end{array}$$E20

Taking the inverse Fourier transform of $\widehat{\stackrel{\sim }{\varphi }}$ yields

$$\stackrel{\sim }{\varphi }={\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}.$$E21

Substituting Eq. (8) and Eq. (21) into Eq. (18) yields

$$\begin{array}{l}-\mu \left(U+\stackrel{\sim }{U}\right)+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\left(\widehat{U}+\widehat{\stackrel{\sim }{U}}\right)\right\}+\left(U+\stackrel{\sim }{U}\right)\left(\beta {\left|U\right|}^2+\gamma {\left|U\right|}^4\right)\\ +2U^2\stackrel{\sim }{U}\left(\beta +2\gamma {\left|U\right|}^2\right)-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}\left(U+\stackrel{\sim }{U}\right)\\ -\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U-V\left(U+\stackrel{\sim }{U}\right)=0.\end{array}$$E22

After grouping the terms, Eq. (22) can be written as

$$\begin{array}{l}\left[-\mathit{\mu U}+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{U}\right\}+\beta {\left|U\right|}^{2}U-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U+\gamma {\left|U\right|}^{4}U-\mathit{VU}\right]\\ +\left[-\mu \stackrel{\sim }{U}+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{\stackrel{\sim }{U}}\right\}+\beta {\left|U\right|}^2\stackrel{\sim }{U}+\gamma {\left|U\right|}^4\stackrel{\sim }{U}+2\beta {\left|U\right|}^2\stackrel{\sim }{U}+4\gamma {\left|U\right|}^4\stackrel{\sim }{U}\right.\\ \left.-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}\stackrel{\sim }{U}-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U-V\stackrel{\sim }{U}\right]=0.\end{array}$$E23

Hence,

$$\begin{array}{l}\left[-\mathit{\mu U}+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{U}\right\}+\beta {\left|U\right|}^{2}U-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U+\gamma {\left|U\right|}^{4}U-\mathit{VU}\right]\\ +\left[-\mu \stackrel{\sim }{U}+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{\stackrel{\sim }{U}}\right\}+3\beta {\left|U\right|}^2\stackrel{\sim }{U}-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}\stackrel{\sim }{U}+5\gamma {\left|U\right|}^4\stackrel{\sim }{U}\right.\\ \left.-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U-V\stackrel{\sim }{U}\right]=0.\end{array}$$E24

From Eq. (9), we know that the first bracket is identically zero. Consequently, we obtain

$$\begin{array}{l}-\mu \stackrel{\sim }{U}+{\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{\stackrel{\sim }{U}}\right\}+3\beta {\left|U\right|}^2\stackrel{\sim }{U}-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}\stackrel{\sim }{U}+5\gamma {\left|U\right|}^4\stackrel{\sim }{U}\\ -\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U-V\stackrel{\sim }{U}=0.\end{array}$$E25

Moreover, Eq. (24) satisfied $L_0\left(U+\stackrel{\sim }{U}\right)=L_1\stackrel{\sim }{U}+O\left({\stackrel{\sim }{U}}^2\right)$. Therefore, to obtain a linearized operator $L_1{\mathrm{L}}_1$, Eq. (25) can be written as

$$L_1\stackrel{\sim }{U}={\mathcal{F}}^{-1}\left\{-\frac{1}{2}\left(k_x^2+k_y^2\right)\widehat{\stackrel{\sim }{U}}\right\}+T_1\stackrel{\sim }{U}-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{2U\stackrel{\sim }{U}\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}U=0,$$E26

where

$$T_1=-\mu +3\beta {\left|U\right|}^2-\rho {\mathcal{F}}^{-1}\left\{\frac{k_x^2\mathcal{F}\left\{{\left|U\right|}^2\right\}}{k_x^2+{\mathit{vk}}_y^2}\right\}+5\gamma {\left|U\right|}^4-V.$$E27

To obtain soliton solution U(x, y) in $L_{0}U=0$, we numerically integrate the following distance-dependent squared-operator evolution equation

$$U_z=-M^{-1}{L}_1^{†}{M}^{-1}{L}_{0}U,$$E28

where ${\left(\cdot \right)}^{†}$ denotes the Hermitian of the operator and M is a real-valued positive definite Hermitian preconditioning operator that is introduced to accelerate the convergence. Since it is easily invertible to take the Fourier transform, we take the preconditioning operator M to be in the form of the following

$$M=c-\left({\mathit{\partial }}_{\mathit{xx}}+{\mathit{\partial }}_{\mathit{yy}}\right),$$E29

where c > 0 is a parameter for parametrizing the numerical scheme. Applying the Fourier transform to Eq. (29) yields

$$\mathcal{F}\left\{M\right\}=c+k_x^2+k_y^2.$$E30

Consequently, in Eq. (28)

$$M^{-1}{L}_1^{†}{M}^{-1}{L}_{0}U={\mathcal{F}}^{-1}\left\{\frac{\mathcal{F}\left\{L_1{M}^{-1}{L}_{0}U\right\}}{c+k_x^2+k_y^2}\right\}.$$E31

Using the forward Euler method, steady-state solution U is computed by an iterative scheme as follows

$$U_{n+1}=U_n-{\left[M^{-1}{L}_1^{†}{M}^{-1}{L}_{0}U\right]}_{U=U_n}\mathrm{\Delta }z,$$E32

where Δz is an auxiliary distance-step parameter. It has been demonstrated that the SOM algorithm converges to a soliton solution for a wide range of nonlinear PDEs if the initial condition is sufficiently close to the exact solution and the distance-step Δz in the iteration scheme is less than a specific threshold value [60, 61]. To obtain a convergent soliton solution, c and Δz are chosen heuristically as positive real numbers. Moreover, our convergence criterion is that the obtained solution satisfies Eq. (10) with an absolute error less than ${10}^{-5}$.

In this chapter, to obtain dipole solitons, the initial condition of the SOM algorithm is chosen as a multi-humped Gaussian function which is given by

$$U_0\left(x,y,0\right)=\sum _{n=0}^{H-1}{e}^{-A\left[{\left(x+x_n\right)}^2+{\left(y+y_n\right)}^2\right]+i{\theta }_n},$$E33

where $x_n$ and $y_n$ represent the location of the solitons on the lattice, H corresponds to the number of humps, A is a positive integer, and ${\theta }_n$ is the phase difference. Since we numerically investigate the dipole solitons, H is set to 2, thus Eq. (33) takes the following form:

$$U_0\left(x,y,0\right)=e^{-A\left[{\left(x+x_0\right)}^2+{\left(y+y_0\right)}^2\right]+i{\theta }_0}+e^{-A\left[{\left(x+x_1\right)}^2+{\left(y+y_1\right)}^2\right]+i{\theta }_1},$$E34

where ($x_0,y_0$) and ($x_1,y_1$) represent the locations of dipole solitons, ${\theta }_0$ and ${\theta }_1$ are the phase differences of dipole solitons. It was shown that the solitons located at the maximum of the lattices are unstable [13, 21, 24], due to this fact we will investigate the dipole solitons located on minima of the considered square lattice. A dipole (two-phased) localized soliton numerically found by

$$A=1,x_n=r\mathit{cos}{\theta }_n,y_n=r\mathit{sin}{\theta }_n,n=0,1.$$E35

Here r is set to be π and ${\theta }_n=\mathit{n\pi }$, so that the humps of the initial condition are located at the local minima of the lattice where $\left(x_0,y_0\right)=\left(\pi ,0\right)$ and $\left(x_1,y_1\right)=\left(-\pi ,0\right)$.

Unless otherwise stated, parameters in the NLSM model (2) are fixed to

$$\left(\mu ,\rho ,v,\beta ,\gamma ,V_0\right)=\left(-0.1,0.5,1.5,2,-0.1,12.5\right).$$E36

It is noted that $\rho =0.5$ and $v=1.5$ are especially chosen to simulate quadratic optical effects in potassium niobate (${\text{KNbO}}_3$) [30].

Dipole solitons of the NLSM model (2) are calculated by the SOM method. In Figure 2, 3D views (first column), phase structures (second column), and contour plots of dipole solitons on the underlying lattice (third column) are displayed for self-defocusing ($\gamma <0$) and self-focusing ($\gamma >0$) quintic nonlinearities. γ is set to be −0.3, −0.1 and 0.3 in Figure 2(a)–(c), respectively, and all other parameters are fixed to the values given in Eq. (36). Figure 2 shows that the dipole solitons can be generated on the lattice minima (see the third column), and the amplitudes of dipole solitons are decreased as γ increased (from −0.3 to 0.3) (see the first column).

Figure 2.

3D dipole profiles centered at the lattice minima (first column), the phase structures of the dipole (second column), and the contour plot of the dipole solitons superimposed on the underlying lattice (third column), which are obtained for (a) $\gamma =-0.3,c=2.1,\mathrm{\Delta }z=0.4$ and error is order of $10^{-6}$, (b) $\gamma =-0.1,c=2.5,\mathrm{\Delta }z=0.4$ and error is order of $10^{-8}$, and (c) $\gamma =0.3,c=2.5,\mathrm{\Delta }z=0.4$ and error is order of $10^{-8}$. All other parameters are fixed to the values in Eq. (36).

3. Stability analysis

The stability dynamics of dipole solitons obtained by the SOM method are studied by the power analysis and direct simulation of the nonlinear evolution.

The power of solitons plays an important role in the stability analysis and it is calculated by

$$P\left(\mu \right)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{\left|U(x,y;\mu )\right|}^2\mathit{\text{dxdy}}.$$E37

Vakhitov and Kolokolov proved a necessary condition for the linear stability of solitons in Ref [62]. They demonstrated that a soliton is linearly stable only if its power increases as propagation constant (or eigenvalue) μ increases. In other words, a necessary condition for the stability of solitons is

$$\frac{\mathit{dP}}{\mathit{d\mu }}>0.$$E38

Moreover, Weinstein and Rose [63, 64] proved that a necessary condition for the nonlinear stability of solitons is also the slope condition given in Eq. (38).

To analyze nonlinear stability of the NLSM model (2), we examine the direct simulation of dipole solitons obtained by the SOM method. A finite-difference discretization scheme is used in the spatial domain (x, y) and the dipole solitons are advanced in the z direction with a fourth-order Runge-Kutta method. The initial condition of the nonlinear evolution is taken to be a dipole soliton, and 1% random noise is inserted into the amplitude of the initial condition.

The power diagrams of dipole solitons are displayed for varied μ, γ, β and ρ values in Figure 3(a)–(d), respectively. It is noted that the domain of existence for the varied parameter is shown on the x-axis of each panel when other parameters are fixed to the values in Eq. (36). Figure 3 shows that the power of dipole solitons increases as μ and ρ increase, whereas the power of dipole solitons decreases as γ and β increase. Moreover, the stability (solid blue) and instability (red dotted) regions of parameters are determined by the nonlinear evolution of dipole solitons for each point on the power curves.

Figure 3.

Power of dipole solitons (a) for varying eigenvalue μ, (b) for varying quintic nonlinearity coefficient γ, (c) for varying cubic nonlinearity coefficient β, and (d) for varying quadratic nonlinear response ρ. The nonlinear stability and instability regions are shown by solid blue and red dotted lines, respectively.

The dipole solitons are found to be nonlinearly stable for self-defocusing quintic nonlinearity ($\gamma =-0.1$) when the power P ∈ [0.99, 1.87] and propagation constant $\mu \in \left[-0.75,-0.6\right]$. Also, the dipole solitons are nonlinearly stable when P ∈ [3.12, 4.26] and $\mu \in \left[-0.35,-0.06\right]$_, which is the second nonlinearly stable gap (see Figure 3(a)). These results are consistent with key analytical results on nonlinear stability, which Weinstein and Rose proved in Ref [63, 64], since slope of the power-eigenvalue (P − μ) diagram is positive. As can be seen from Figure 3(b), the dipole solitons are obtained for $\gamma \in \left[-0.7,25\right]$, when other parameters are fixed, and dipole solitons are nonlinearly stable for $\gamma \in \left[-0.21,0.25\right]$. Zoom-in view of this stability domain is depicted in Figure 3(b). Furthermore, it is observed that dipole solitons are nonlinearly stable for $\beta \in \left[1.6,18.9\right]$ (see Figure 3(c)), and dipole solitons are stable for $\rho \in \left[0,0.8\right]$ (see Figure 3(d)) in their existence domains when other parameters are fixed.

In Figure 4, nonlinear evolution of peak amplitudes, 3D views of the evolved dipole solitons, and the phase structures of evolved dipole solitons are plotted for the dipole solitons that are shown in Figure 2. The effect of quintic nonlinearity (γ) on nonlinear stability is investigated by fixing other parameters as in Eq. (36).

Figure 4.

Nonlinear evolution of maximum amplitudes as a function of propagation distance z (first column), 3D views of the dipole solitons after evolution (second column), and the phase structures of dipole solitons after evolution (third column) for (a) $\gamma =-0.3$, (b) $\gamma =-0.1$, and (c) $\gamma =0.3$. All other parameters are taken as in Eq. (36).

Figure 4(b) shows that peak amplitudes of dipole solitons oscillate mildly (first column), and the 3D profile (second column) and phase structure (third column) of dipole solitons are preserved for $\gamma =-0.1$. Thus, the dipole solitons are nonlinearly stable for the defocusing quintic nonlinearity for the considered parameter regime. On the other hand, as shown in Figure 4(a) and (c), when the quintic nonlinearity is strong ($\gamma =-0.3$ and $\gamma =+0.3$), peak amplitudes of dipole solitons increase significantly in a short propagation distance z (first column), dipole profiles (second column) cannot be preserved, and phase structures of dipole solitons (third column) break up after evolution. Comparing Figure 4(a) and (c), it is observed that the propagation distance of dipole solitons in a medium with strong self-focusing nonlinearity ($\gamma =0.3$) is longer than that of a medium with strong self-defocusing nonlinearity ($\gamma =-0.3$). Considering these evolution results in Figure 4 and the existing domain for γ in Figure 3(b), it is demonstrated that both strong self-focusing and self-defocusing quintic nonlinearities have a negative effect on the nonlinear stability of dipole solitons.

In previous studies [6, 14, 24], it is found that modification of the depth of potential can suppress nonlinear instabilities. More recently, Bağcı and coworkers [46] have demonstrated that nonlinear stability of fundamental solitons in an NLSM system (2) with quintic nonlinearity can be improved by the modification of lattice depth $V_0$. They showed that increased lattice depth supports the stability of fundamental solitons in a medium with strong self-focusing ($\gamma =0.3$) quintic nonlinearity, and the stability of solitons in a medium with strong self-defocusing ($\gamma =-0.3$) quintic nonlinearity can be improved by decreasing lattice depth. For the dipole solitons, evolution of peak amplitudes is depicted for varying potential depths, when $\gamma =-0.3$ and $\gamma =0.3$ in Figure 5(a) and (b), respectively. Figure 5(a) shows that the stability of dipole solitons is improved by decreasing lattice depth (from 25 to 5) for strong self-defocusing nonlinearity ($\gamma =-0.3$), and collapse can be arrested when $V_0=5$. In contrast, as shown in Figure 5(b), the propagation distance of dipole solitons in a medium with strong self-focusing nonlinearity ($\gamma =0.3$) is extended by increasing lattice depth (from 5 to 50). It should be noted that these results are in agreement with the findings of the aforementioned studies. Thus modification of the lattice depth can be utilized to improve the nonlinear stability of dipole solitons.

Figure 5.

Maximum amplitudes of the evolved dipole solitons for varying depth of potential V₀, when the dipole soliton is obtained for (a) $\gamma =-0.3$ and (b) $\gamma =0.3$.

It is also known that when the quadratic [15, 24, 44] and quintic [46] electro-optic effects are strong, the instability of fundamental solitons can be improved by increasing the anisotropy parameter. To examine the effect of anisotropy coefficient v on the nonlinear stability of dipole solitons in a medium with strong quintic nonlinearity ($\gamma =-0.3$ and $\gamma =0.3$), evolution of the peak amplitudes is displayed for varied v values in Figure 6. Figure 6 shows that increasing the anisotropy coefficient v from 0.001 to 10 stabilizes the dipole solitons in a medium with strong self-defocusing nonlinearity ($\gamma =-0.3$), and increasing v from 0.001 to 1000 extends the propagation distance of dipole solitons in a medium with strong self-focusing nonlinearity ($\gamma =0.3$). Thus, larger anisotropy coefficient supports the nonlinear stability of the dipole solitons, and this result complies with the results of the previous studies [46]. It is important to note that the parameters ρ and v are predetermined coefficients that depend on the type of optical materials; larger values of v cannot be applied to real optical systems. In this chapter, the effect of extremely large v values on the stability of dipole solitons is explored numerically.

Figure 6.

Maximum amplitudes of the evolved dipole solitons for varying anisotropy coefficients v, when the dipole soliton is obtained for (a) $\gamma =-0.3$, and (b) $\gamma =0.3$.

4. Conclusions

In this chapter, the existence and nonlinear stability dynamics of dipole solitons have been investigated for a nonlocal nonlinear medium with quintic nonlinear response. This medium was characterized by the (2+1)D NLSM system with a periodic external lattice. Dipole solitons were obtained for self-defocusing ($\gamma <0$) and self-focusing ($\gamma >0$) quintic nonlinearities by the SOM method, and the nonlinear stability of these dipole structures has been investigated by the direct simulation of the model equations. Power of dipole solitons was determined for varying μ, γ, β, and ρ parameters and it has shown that the power of dipole solitons increases as the eigenvalue μ and quadratic nonlinear response ρ increase, whereas the power of dipole solitons decreases as quintic nonlinearity coefficient γ and cubic nonlinearity coefficient β increase.

Nonlinear evolution of the dipole solitons showed that the dipole solitons are stable for the weak self-focusing and self-defocusing quintic nonlinearity. In other words, as an absolute value of γ increases, the obtained dipole solitons become nonlinearly unstable in both self-focusing and self-defocusing media. It has been demonstrated that the collapse of dipole solitons can be arrested by decreased potential depth in a medium with strong self-defocusing quintic nonlinearity ($\gamma =-0.3$), while the deeper lattice extends the propagation distance of dipole solitons in a medium with strong self-focusing quintic nonlinearity ($\gamma =0.3$). Furthermore, it has been observed that increasing the anisotropy coefficient (v) extends the propagation distance of the dipole solitons for strong self-focusing quintic nonlinearity, and it stabilizes the dipole solitons for strong self-defocusing quintic nonlinearity.

In conclusion, the existence and stability properties of dipole solitons have been numerically explored in a nonlocal nonlinear medium with quintic nonlinear response, and it has been demonstrated that the instability of dipole solitons can be suppressed by modification of the lattice depth and increased anisotropy coefficient.

Conflict of interest

The authors declare no conflict of interest.