Stabilization of Two-Dimensional Optical Continuous-Wave States by a Potential Trough

1. Introduction and the Models

Spatial-domain propagation of light in bulk waveguides is determined by the interplay of diffraction and material nonlinearity. This setting gives rise to a variety of self-supporting modes, including solitons [1,2,3], vortices [4,5,6,7] (similar to their counterparts in atomic Bose-Einstein condensates (BEC) [8,9]), skyrmions [10,11], etc., which are adequately modeled by the two-dimensional (2D) nonlinear Schrödinger (NLS) equation, or a coupled system of such equations [12,13]. Spatiotemporal light propagation in planar waveguides is accurately modeled by 2D NLS equations which are essentially identical to their counterparts in the spatial domain [14]. While the ubiquitous nonlinearity of optical materials is represented by the Kerr cubic term, accurate experiments were also performed in chalcogenide glasses [15] and colloidal suspensions of metallic particles [16], which are adequately represented by combined nonlinearities – most typically, cubic-quintic (CQ) ones, which feature the competition of the cubic self-focusing and quintic defocusing.

While 1D solitons are typically stable states, the stability is a major challenge for studies of 2D nonlinear modes in optics, as well as in other areas of physics, especially in BEC [14,17]. The basic problem is the critical collapse which destabilizes 2D fundamental solitons under the action of the cubic self-focusing [12,13], while the solitons with embedded vorticity are still more vulnerable to the splitting instability [18]. Two major mechanisms which provide stabilization of both the fundamental (zero-vorticity) and vortex optical solitons in 2D rely upon the use of the CQ nonlinearity and/or effective potentials, which represent the transverse modulation of the refractive index in the waveguide. In particular, the CQ nonlinearity readily stabilizes a part of the vortex-soliton states, with different values of the topological charge (vorticity), against the splitting [19,20]. On the other hand, spatially periodic (lattice) potentials, both fully two-dimensional (2D) [21,22] and quasi-one-dimensional (Q1D) ones, which depend on a single coordinate [23,24]), also secure the stabilization of the fundamental and vortex solitons (although the lattice potentials destroy the underlying spatial isotropy, vortex solitons with a nonzero topological charge can be defined in this case too [21,22]).

The most basic 2D model combining the CQ nonlinearity and spatial inhomogeneity includes a Q1D potential trough, which depends on the single transverse coordinate, x, and does not depend on the other coordinate, y [25]. The corresponding scaled NLS equation for the slowly varying amplitude of the optical wave, $u\left(x,y;z\right)$, is

$$i{u}_z=-{\displaystyle \frac{1}{2}}\left(u_{xx}+u_{yy}\right)-{\left|u\right|}^{2}u+g{\left|u\right|}^{4}u-W_{0}exp\left(-x^2\right)u,$$

(1)

where z is the propagation distance, $g>0$ is the strength of the quintic self-defocusing nonlinearity (while the strength of the cubic self-focusing is set to be 1 by scaling), and $W_0>0$ is the depth of the potential trough with the Gaussian profile, whose width is also fixed to be 1 by means of scaling.

Equation (1) applies as well to the spatiotemporal light propagation in a planar waveguide with anomalous group-velocity dispersion, In that case, x remains the transverse coordinate, while $y=t-z/V_{\mathrm{gr}}$ is the temporal variable, with t and $V_{\mathrm{gr}}$ being the scaled time and group velocity of the carrier wave. In the spatiotemporal setting, the trough potential represents a stripe of a material with an enhanced refractive index, embedded into the planar waveguide. In the latter case, solely the Q1D effective potential is possible, because it cannot be made a function on time.

Fully localized states, predicted as solutions of Eq. (1), are characterized by the value of the 2D integral power in the spatial domain (or energy in the spatiotemporal one),

$$P_{2\mathrm{D}}={\int }_{-\infty }^{+\infty }dy{\int }_{-\infty }^{+\infty }dx{\left|u\left(x.y\right)\right|}^2,$$

(2)

which is a dynamical invariant of Eq. (1). Other conserved quantities are the Hamiltonian,

$$H={\int }_{-\infty }^{+\infty }dy{\int }_{-\infty }^{+\infty }dx\left[{\displaystyle \frac{1}{2}}\left({\left|u_x\right|}^2+{\left|u_y\right|}^2\right)-{\displaystyle \frac{1}{2}}{\left|u\right|}^4+{\displaystyle \frac{g}{3}}{\left|u\right|}^6-W_{0}exp\left(-x^2\right){\left|u\right|}^2\right],$$

and the y-component of the momentum, $M_y=i{\int }_{-\infty }^{+\infty }dy{\int }_{-\infty }^{+\infty }dx{u}_y^{*}u$, where * stands for the complex conjugate.

The simplest stationary states admitted by Eq. (1) are continuous waves (CWs), which are solutions trapped in the potential trough and thus localized in the x direction, while their y dependence is represented by a real wavenumber q:

$$u_{\mathrm{CW}}\left(x,y,z\right)=exp\left(ikz+iqy\right)U\left(x\right).$$

(3)

Here, k is a real propagation constant, and a real modal function $U\left(x\right)$ obeys the the transverse equation:

$$\left(k+{\displaystyle \frac{1}{2}}{q}^2\right)U={\displaystyle \frac{1}{2}}{\displaystyle \frac{d^{2}U}{d{x}^2}}+U^3-g{U}^5+W_{0}exp\left(-x^2\right)U.$$

(4)

Using the obvious Galilean invariance of Eq. (1) in the y direction, it is sufficient to consider the CW solutions with $q=0$, which is done below, the respective CW solution being independent of y. Then, Eq. (4) produces bound states (solutions localized at $\left|x\right|\to \infty$) for $k>0$.

Equation (4) is the 1D stationary NLS equation with the Gaussian trapping potential, Therefore, it is natural to classify its bound-state solutions as a nonlinear extension of eigenstates which are produced by the linear Schrödinger equation with the trapping potential, in the small-amplitude limit. Accordingly, the bound states are classified as the spatially even ($U(-x)=U\left(x\right)$) ground state (GS), spatially odd ($U(-x)=-U\left(x\right)$) first excited state (alias the dipole mode, DM), spatially even second excited state (alias the quadrupole mode), etc.

The basic problem for CW solutions is their modulational instability (MI, alias the Benjamin-Feir instability [26]) against small perturbations which impose periodic modulation of the CW state along the free coordinate y [27,28,29,30,31,32,33]. In particular, MI plays the crucial role in the generation of rogue waves [34,35,36,37,38]. For the CW solutions of Eq. (1) with the GS transverse structure, this problem was originally addressed in Ref. [25]. In the presents work, we aim to further elaborate this topic. First, in Section 2 we recapitulate the MI analysis performed in the framework of Eq. (1), including new findings for the CWs with the DM transverse structure. While all the CWs are subject to MI, the instability is very weak for the CWs of both the GS and DM types, with values of the propagation constant close to existence boundaries of these solution families, the respective interval of the “practical stability" being wider for the DMs than for GSs.

It is natural to expect that MI may split the CW into a periodic chain of 2D solitons along the y coordinate, trapped in the potential trough. In Section 3, we produce a numerical solution for the stationary soliton chain and demonstrate its stability. We also report results produced by the application of a kick to a particular soliton in the chain. The kick, that may be directed along x or y, eventually generates a complex dynamical picture.

The most essential novel results are reported in Section 4 for the 2D model with the CQ nonlinearity and the potential trough, $W\left(x\right)$, represented by the delta-functional profile:

$$i{u}_z=-{\displaystyle \frac{1}{2}}\left(u_{xx}+u_{yy}\right)-{\left|u\right|}^{2}u+g{\left|u\right|}^{4}u-\sqrt{\pi }{W}_0\delta \left(x\right)u.$$

(5)

Factor $\sqrt{\pi }$ is introduced in Eq. (5) to keep the same transverse area of the trough potential, ${\int }_{-\infty }^{+\infty }W\left(x\right)dx=\sqrt{\pi }{W}_0$, as in Eq. (1). We obtain two families of exact analytical solutions of Eq. (5), in parameter regions

$$W_0<W_{\mathrm{cr}}\left(\mathbf{Case}\mathbf{1}\right)\mathrm{and}{W}_0>W_{\mathrm{cr}}\left(\mathbf{Case}\mathbf{2}\right),$$

(6)

where the critical potential strength is

$$W_{\mathrm{cr}}=\sqrt{3/\left(8\pi g\right)}.$$

(7)

The numerical analysis demonstrates that, while the exact solutions are subject to MI in Case 1, they are completely stable in Case 2. The possibility to produce stable CW states in the form of Q1D bright solitons, supported by the trough potential, is a remarkable finding (in particular, because the stable solutions are obtained here in the analytical form), which was not reported for previously studied models. In this connection, it is relevant to mention that it was demonstrated in Ref. [39] that a Q1D (stripe-shaped) dark soliton may be stabilized against MI, in the framework of the two-dimensional NLS equation with the cubic self-defocusing term (the quintic one is not needed in this case), by means of a sufficiently strong Q1D potential ridge (rather than the trough), which corresponds to $W_0<0$ in Eq. (1).

The paper is concluded by Section 5, which summarizes the obtained results and puts forward possibilities for additional studies on the present topic.

2. MI (Modulational Instability) of the Q1D (Quasi-One-Dimensional) CW States

2.1. Linearized Equations for Modulational Perturbations

To analyze the modulational (in)stability of the CW states, a perturbed solution is looked for as

$$u\left(z,x,y\right)=exp\left(ikz\right)\left\{U\left(x\right)+\epsilon \left[a\left(x\right)exp\left(\gamma z+ipy\right)+b^{*}\left(x\right)exp\left({\gamma }^{*}z-ipy\right)\right]\right\},$$

(8)

where $U\left(x\right)$ is a real solution of Eq. (4) with propagation constant k (recall we set $q=0$ in Eqs. (3) and (4)), $a\left(x\right)$ and $b\left(x\right)$ are components of the eigenmode of small perturbations, with a real infinitesimal amplitude $\epsilon$ and wavenumber p, and $\gamma \left(p\right)$ is the corresponding eigenvalue, the stability condition being $\mathrm{Re}\left(\gamma \left(p\right)\right)=0$ for all real p. The substitution of ansatz (8) in Eq. (1) and linearization with respect to $\epsilon$ leads to the system of linear equations for the perturbation eigenmode (the Bogoliubov-de Gennes (BdG) equations, in terms of the BEC theory [40]):

$$\begin{array}{c}\left(-k-{\displaystyle \frac{p^2}{2}}+i\gamma \right)a=\end{array}$$

$$\begin{array}{c}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^{2}a}{d{x}^2}}-W_{0}exp\left(-x^2\right)a-2U^2\left(x\right)a-U^2\left(x\right)b+3g{U}^4\left(x\right)a+2g{U}^4\left(x\right)b,\\ \left(-k-{\displaystyle \frac{p^2}{2}}-i\gamma \right)b=\end{array}$$

(9)

$$\begin{array}{c}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^{2}b}{d{x}^2}}-W_{0}exp\left(-x^2\right)b-2U^2\left(x\right)b-U^2\left(x\right)a+3g{U}^4\left(x\right)b+2g{U}^4\left(x\right)a.\end{array}$$

(10)

The analysis of the CW states and their stability may be performed treating $W_0>0$ and $g>0$ as control parameters, while $k>0$ is a free parameter of the family of the CW states.

2.2. The MI of the CW States with the Ground-State (GS) Transverse Structure

In the uniform space, with $W_0=0$ in Eq. (1), the transverse profile of the CW state is given by the well-known Q1D soliton solution of Eq. (1) [41]:

$$u_{\mathrm{sol}}(x,z)=e^{ikz}\sqrt{{\displaystyle \frac{2k}{1+\sqrt{1-4gk/3}cosh\left(2\sqrt{k}x\right)}}},$$

(11)

in which the propagation constant k takes values in the interval

$$0<k<k_{max}\equiv {\displaystyle \frac{3}{4g}}.$$

(12)

The total 1D power of soliton (11),

$$P\left(k\right)={\int }_{-\infty }^{+\infty }{\left|u_{\mathrm{sol}}\left(x\right)\right|}^{2}dx={\displaystyle \frac{\sqrt{3}}{2\sqrt{g}}}ln\left({\displaystyle \frac{\sqrt{3}+2\sqrt{kg}}{\sqrt{3}-2\sqrt{kg}}}\right),$$

(13)

diverges in the limit of $k\to k_{max}$, when the soliton expands into a CW state in the x direction.

Typical examples of the stationary GSs, produced by the numerical solution of Eq. (4) with $q=0$ and $W_0>0$, are plotted in Figure 1. The solutions were produced by means of the modified squared-operator method (MSOM) [42] in the spatial domain $\left(-X/2,+X/2\right)$ of size $X=40$, discretized on a mesh of 128 points, with the mesh size $\Delta x=0.3125$. They may be construed as solitons trapped in the potential well. For fixed values of parameters $W_0$ and g, the CW family is characterized by the dependence of the 1D power P on k, cf. Eq. (13) for $W_0=0$. In particular, a set of $P\left(k\right)$ curves for $g=0.5$ and different values of $W_0$ is plotted in Figure 2. At $P=0$, each curve starts from the eigenvalue which, in terms of quantum mechanics, is determined by the GS energy, $E_{\mathrm{GS}}\equiv -k(P=0)$, of the potential well $W\left(x\right)=-W_{0}exp\left(-x^2\right)$. With the increase of of P, the curves with $W_0>0$ attain their turning points, $k=k_{\mathrm{turn}}$, and then move back, in the direction of $k<k_{\mathrm{turn}}$, as shown by the dashed lines in Figure 2. The extension of the curves towards $P\to \infty$ will bring them asymptotically close to the value $k=3/\left(4g\right)$, see Eq. (12).

In agreement with Ref. [25], the numerical solution of the BdG system of Eqs. (9) and () demonstrates that all CW modes of the GS type are indeed subject to the MI, as the solution always produces eigenvalues with $\mathrm{Re}\left(\gamma \left(p\right)\right)\ne 0$. As an example, dependences of $\mathrm{Re}\left(\gamma \right)$ on p for $W_0=5.0$ and $g=0.5$ are plotted in Figure 3, for different values of the 1D power P.

The CW instability, predicted by the MI analysis, was confirmed by direct simulations of the perturbed evolution of the CW states, which were carried out in the framework of Eq. (1), by means of the usual split-step Fourier-transform method (here and below, the simulations were performed with stepsize $\Delta z=0.001$). In particular, results of the simulations for the CW state with the central GS profile from Figure 1, which is characterized by he largest MI gain, $\mathrm{Re}{\left(\gamma \right)}_{max}=0.8649$, for the parameters fixed in Figure 1, viz., $W_0=5.0$ and $g=0.5$, are displayed in Figure 4. It demonstrates that, in the $\left(x,y\right)$ plane, the instability splits the original Q1D CW state into a chain of localized 2D speckles, which performs irregular oscillations.

The CW states corresponding to other values of eigenvalue k in Eq. (4), which are taken closer to the linear limit, i.e., lower edge of the respective $P\left(k\right)$ curve (see Figure 2), demonstrate much weaker MI, so that, in direct simulations, they may seem as practically stable states. An example is presented in Figure 5, for $W_0=2.5$, $g=1.5$, and $P=1.3395$, $k=1.845$, It is obvious that this CW state, for which the numerical solution of the BdG equations (9) and () yields a very small peak value of the MI gain, $\mathrm{Re}{\left(\gamma \right)}_{max}=0.0343$, indeed exhibits virtually no instability. On the other hand, CWs with large values of P, such as $P=P=3.1453$ and $5.3634$ in Figure 3, may also demonstrate practical stability, explained by the domination of the self-defocusing quintic nonlinearity in this case.

2.3. The MI of the Dipole-Mode (DM) CWs

CW states of the DM type were not studied in previous works in the framework of Eqs. (1) and (4), or in similar models (cf. Ref. [25]). Note that, unlike the 1D soliton (11), Eq. (1) in the free space (with $W_0=0$) does not have solutions of the DM type.

Families of numerically found stationary CWs of the DM type are represented by the respective $P\left(k\right)$ curves in Figure 6(a) for $W_0=5.0$ and different values of g, where $P={\int }_{-\infty }^{+\infty }{\left|u\left(x\right)\right|}^{2}dx$ is the same 1D power as defined above. Like their GS counterparts (cf. Figure 2), the curves originate, at $P=0$, from eigenvalues k corresponding to the first excited state in the quantum-mechanical potential $W\left(x\right)=W_{0}exp\left(-x^2\right)$, reach turning points at $k=k_{\mathrm{turn}}$, and with further increase of P move back in the direction of $k<k_{\mathrm{turn}}$, as shown by dashed lines in Figure 6(a). Figure 5(b) demonstrates the DM profiles for two widely different values of the 1D power, $P=1.1011$ and $8.4537$, which corresponds to markers a and b in panel (a), for $g=1.0$. Lastly, Figure 5(c) exhibits spectra of the MI gain for these CW states.

In agreement with the predictions of the MI analysis, direct simulations of the perturbed evolution of the CW states of the DM type demonstrate strong instability when the respective MI gain, $\mathrm{Re}{\left(\gamma \right)}_{max}$, is relatively large, and “practical stability" when the MI gain is small, see typical examples displayed in Figure 7 and Figure 8, respectively. The strong MI splits the flat CW into a randomly oscillating chain of localized speckles, which keep the original dipole structure. It is relevant to stress that the region of the “practical stability" for the CW states of the DM type is essentially broader than for their GS counterparts. Similar to what is concluded above for the CW states of the GS type, the effective stabilization for the DM CWs may be provided, as in Figure 8, by a sufficiently large value of strength g of the defocusing term ($g=3$ in Figure 8), or by a large value of the 1D power, which also makes the defocusing nonlinear term a dominant one. For instance, in the same case of $W_0=5.0$ and $g=1.0$ as in Figure 7, but with the power larger by a factor $\sim 3$, the CW of the DM types seems to be fully stable in direct simulations (not shown here in detail)..

We have also investigated CW states which correspond, in the linear limit, to the second-order (quadrupole) bound state in the trapping potential $W\left(x\right)=W_{0}exp\left(-x^2\right)$. They always demonstrate strong MI, without any interval of the “practical stability", therefore they are not considered in detail here.

3. Stationary Chains of 2D Solitons

3.1. The Existence and Stability of the Chains

While, as shown in Figure 4, the MI splits the unstable CW states with the GS transverse structure into a chain of localized speckles in the state of random oscillations, it is natural to expect that the same Eq. (1) may produce a stationary chain of 2D solitons, arranged periodically along the y axis. Such a chain may be considered as a manifestation of the MI in terms of stationary states.

Here, we report numerical results for the chains, systematically collected for parameters $W_0=3$ and $g=1$ in Eq. (1), which makes it possible to exhibit generic findings. The stationary-chain solutions with propagation constant k were looked for as

$$u\left(x,y,z\right)=exp\left(ikz\right)U\left(x,y\right),$$

(14)

with real function $U(x,y)$ satisfying the equation (written with the above-mentioned values $W_0=3$, $g=1$)

$$kU={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\partial }^{2}U}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}U}{\partial y^2}}\right)+U^3-U^5+3exp\left(-x^2\right)U,$$

(15)

cf. Eq. (4). The chain solutions were obtained by means of MSOM applied to Eq. (15), starting with the y-periodic input,

$$U_0(x,y)=A_{0}exp\left(-{\displaystyle \frac{x^2}{4}}\right)cos\left({\displaystyle \frac{\pi y}{2}}\right).$$

(16)

The numerical solutions were constructed in the computational domain $-32\le x,y\le +32$ with $256\times 256$ grid points, periodic boundary conditions in the $y$ direction, and zero boundary conditions at $x=\pm 32$. The same boundary conditions were adopted in dynamical simulations of Eq. (1).

The soliton-chain family, produced by the numerical solution of Eq. (15), is presented in Figure 9 by means of the dependence of the 2D integral power (2) on the propagation constant k. The curve features a turning point at $k\approx 0.921$, which connects two branches with opposite slopes. The bottom branch originates, at $P=0$, from $k=k_{\max }\approx 0.735$, which corresponds to the eigenstate of the respective 2D linear Schrödinger equation. Stability of the soliton chains belonging to the continuous blue part of the $P\left(k\right)$ curve in Figure 9 was established by the computation of the respective eigenvalues $\gamma$, as per the BdG system of Eqs. (9) and (), and verified by direct simulations of the perturbed propagation.

Two representative chain solutions, belonging to the bottom and top branches of the $P\left(k\right)$ curve, which are marked (a) and (b) in Figure 9, are displayed in Figure 10, by means of the respective distributions of $\left|U\left(x,y\right)\right|$, with amplitudes ${\left(U_{max}\right)}_{\mathrm{a}}=0.5405$ and ${\left(U_{max}\right)}_{\mathrm{b}}=1.2666$ (in agreement with the fact that the power at point (b) in Figure 9 is essentially higher than at point (a)). Adjacent peaks in both chains keep alternating signs, inherited from the input given by Eq. (16).

3.2. Response of the Soliton Chains to Localized Kicks

To further examine the robustness of the solitons chains and their intrinsic dynamics, we applied a localized phase kick to the single soliton, whose center is located at $\left(x,y\right)=\left(0,0\right)$ in the stable chain marked (b) in Figure 9 (the stationary shape of the chain is plotted in Figure 10(b)). The evolution excited by the kick was observed in direct simulations of Eq. (1).

First, the kick is applied in the y direction, which is represented by the initial condition

$$u(x,y;z=0)=U(x,y)\left[1+\left(e^{iQy}-1\right)\xi (x,y)\right],$$

(17)

where $U(x,y)$ is the stationary soliton-chain solution,

$$\xi (x,y)=exp\left(-{\displaystyle \frac{x^2}{{\sigma }_x^2}}-{\displaystyle \frac{y^2}{{\sigma }_y^2}}\right)$$

(18)

is the 2D Gaussian mask, with widths ${\sigma }_x$ and ${\sigma }_y$ matched to the transverse dimensions of the individual soliton in the chain, and real Q is the strength of the kick. Here, we display the results for $Q=0.9$, which represents a generic case.

The results of the simulations, up to $z=120$, are displayed in Figure 11. The chain structure remains essentially undisturbed at $z\lesssim 40$, which is followed by the emission of a disturbance wave from the kicked soliton in the $\pm y$ directions. Due to the periodic boundary conditions, the disturbances re-enter from the opposite edges of the computational domain, making the cumulative effect visible at $z\approx 60$. As a result, an irregular pattern appears at $z\gtrsim 80$. This effect depends on the size (period) of the domain in the y direction, as confirmed by simulations performed in the domain with the double size, $\left|y\right|\le 64$: the onset of the complex pattern is delayed, and the disturbance waves re-enter from the opposite edge at $z\approx 120$ (not shown here in detail).

Figure 12 displays the results produced by the application of the kick in the x direction, defined by the initial condition

$$u(x,y;z=0)=U(x,y)\left[1+\left(e^{iQx}-1\right)\xi (x,y)\right],$$

(19)

again with $Q=0.9$. The ensuing picture is qualitatively similar to that in Figure 11, which was initiated by the kick applied in the y direction.

4. CWs Maintained by the Delta-Functional Trapping Potential

4.1. Analytical Results

The substitution of the usual ansatz for the stationary CW state, $u\left(x,y,z\right)=exp\left(ikz\right)U\left(x\right)$, in Eq. (5) leads to the real equation for $U\left(x\right)$, cf. Eq. (4):

$$kU={\displaystyle \frac{1}{2}}{\displaystyle \frac{d^{2}U}{d{x}^2}}+U^3-g{U}^5+\sqrt{\pi }{W}_0\delta \left(x\right)U.$$

(20)

The GS solution to Eq. (20) for an even function $U\left(x\right)$ must be continuous at $x=0$, with the jump of the first derivative, imposed by the delta-functional potential:

$${\displaystyle \frac{dU}{dx}}\left(x=+0\right)-{\displaystyle \frac{dU}{dx}}\left(x=-0\right)=-2\sqrt{\pi }{W}_{0}U(x=0).$$

(21)

A commonly known fact is that the linearized version of Eq. (20) supports precisely one bound state, in the form of

$$U\left(x\right)=U_{0}exp\left(-\sqrt{\pi }{W}_0\left|x\right|\right),\mathrm{with}k=(\pi /2)W_0^2.$$

(22)

It is possible to construct two exact solutions for the soliton pinned to the delta-functional potential. First, in Case 1, defined as per Eqs. (6) and (7), i.e., for $W_0<\sqrt{3/\left(8\pi g\right)}$, the exact solution can be constructed in the following interval of values of the propagation constant,

$$(\pi /2)W_0^2\le k<3/\left(16g\right),$$

(23)

where the left edge naturally coincides with the eigenvalue (22) of the linearized equation (20). If condition (23) holds, one easily finds a solution, combining the usual exact soliton (11) in the free space (with $W_0=0$) and the jump condition (21):

$$U_{\mathrm{sol}}^{\left(1\right)}\left(x,{\xi }_1\right)=\sqrt{{\displaystyle \frac{4k}{1+\sqrt{1-\frac{16}{3}gk}cosh\left(2\sqrt{2k}\left(\left|x\right|+{\xi }_1\right)\right)}}},$$

(24)

with the spatial shift ${\xi }_1>0$ determined as a root of equation

$$\sqrt{2k\left(1-{\displaystyle \frac{16}{3}}gk\right)}sinh\left(2\sqrt{2k}\xi \right)=\sqrt{\pi }{W}_0\left[1+\sqrt{1-{\displaystyle \frac{16}{3}}gk}cosh\left(2\sqrt{2k}\xi \right)\right].$$

(25)

An explicit solution of Eq. (25) is given by a cumbersome expression:

$${\xi }_1={\displaystyle \frac{1}{2\sqrt{2k}}}ln\left[{\displaystyle \frac{\sqrt{1+\left(1-\frac{16}{3}gk\right)\left(\frac{2k}{\pi W_0^2}-1\right)}+1}{\sqrt{1-\frac{16}{3}gk}\left(\frac{\sqrt{2k}}{\sqrt{\pi }{W}_0}-1\right)}}\right].$$

(26)

In the opposite Case 2, also defined as per Eqs. (6) and (7), i.e., as $W_0>\sqrt{3/\left(8\pi g\right)}$, another exact solution for the soliton pinned to the delta-functional potential can be found in the following interval of values of the propagation constant:

$$3/\left(16g\right)<k\le (\pi /2)W_0^2,$$

(27)

cf. interval (23). In this case, the right edge of the existence interval coincides with the eigenvalue (22) of the linearized equation (20). The solution existing in interval (27) is based on the following singular exact solution of Eq. (20) in the free space (with $W_0=0$),

$$U_{\mathrm{sing}}\left(x\right)=\sqrt{{\displaystyle \frac{4k}{1+\sqrt{\frac{16}{3}gk-1}sinh\left(2\sqrt{2k}x\right)}}}.$$

(28)

In the free space, the singular solution is irrelevant. However, in the presence of the delta-functional potential, an exact solution for the pinned soliton can be obtained from expression (28), replacing x by $\left|x\right|+{\xi }_2$:

$$U_{\mathrm{sol}}^{\left(2\right)}\left(x,{\xi }_2\right)=\sqrt{{\displaystyle \frac{4k}{1+\sqrt{\frac{16}{3}gk-1}sinh\left(2\sqrt{2k}\left(\left|x\right|+{\xi }_2\right)\right)}}}$$

(29)

(cf. expression (24) for the first soliton solution), with ${\xi }_2$ defined by the equation

$$\sqrt{2k\left({\displaystyle \frac{16}{3}}gk-1\right)}cosh\left(2\sqrt{2k}{\xi }_2\right)=\sqrt{\pi }{W}_0\left[1+\sqrt{{\displaystyle \frac{16}{3}}gk-1}sinh\left(2\sqrt{2k}{\xi }_2\right)\right],$$

(30)

cf. Eq. (25) for the first soliton solution. The explicit solution of Eq. (30) is

$${\xi }_2={\displaystyle \frac{1}{2\sqrt{2k}}}ln\left[{\displaystyle \frac{\sqrt{1+\left(\frac{16}{3}gk-1\right)\left(1-\frac{2k}{\pi W_0^2}\right)}-1}{\sqrt{\frac{16}{3}gk-1}\left(1-\frac{\sqrt{2k}}{\sqrt{\pi }{W}_0}\right)}}\right],$$

(31)

cf. Eq. (26). Note that, unlike positive ${\xi }_1$ given by Eq. (26) for the first soliton solution, Eq. (31) may produce ${\xi }_2<0$.

4.2. Numerical Results

The accuracy and stability of the above analytical solutions were verified by comparison with numerical solutions obtained from Eqs. (5) and (20) with the delta-function approximated by a narrow Gaussian,

$$\delta \left(x\right)\to \tilde{\delta }\left(x\right)={\displaystyle \frac{1}{\sqrt{\pi }{x}_0}}exp\left(-{\displaystyle \frac{x^2}{x_0}}\right),$$

(32)

with a small regularization scale $x_0={\left(2\Delta x\right)}^2$, where $\Delta x=0.3125$. is the above-mentioned mesh size adopted in realization of the MSOM algorithm. This value is definitely much smaller than characteristic length scales, in the x direction, of all solutions considered here. However, this approximation cannot exactly reproduce the cusp at the center of exact analytical solutions (24) and (29) for the pinned solitons, see Figure 14 below.

Figure 13. An example of the MI destroying the CW based on the stationary profile (24) (corresponding to Case 1, see Eq. (6)). The parameters are $W_0=0.1$, $g=2$, $k=0.05$. (a) The initial intensity profile ${\left|u(x,y;z=0)\right|}^2$. (b) The MI gain, $\mathrm{Re}\left(\gamma \right)$, versus p, showing the instability band with the peak value $\mathrm{Re}{\left(\gamma \right)}_{max}\approx 0.037$, attained at $p\approx 0.27$. (c) The evolution of the local intensity ${\left|u(x,y=0,z)\right|}^2$ in cross-section $y=0$, initiated by the perturbed input (33) with $p=0.2688$ and $\epsilon ={10}^{-3}$. The CW breakup commences at $z\approx 150$.

Figure 13. An example of the MI destroying the CW based on the stationary profile (24) (corresponding to Case 1, see Eq. (6)). The parameters are $W_0=0.1$, $g=2$, $k=0.05$. (a) The initial intensity profile ${\left|u(x,y;z=0)\right|}^2$. (b) The MI gain, $\mathrm{Re}\left(\gamma \right)$, versus p, showing the instability band with the peak value $\mathrm{Re}{\left(\gamma \right)}_{max}\approx 0.037$, attained at $p\approx 0.27$. (c) The evolution of the local intensity ${\left|u(x,y=0,z)\right|}^2$ in cross-section $y=0$, initiated by the perturbed input (33) with $p=0.2688$ and $\epsilon ={10}^{-3}$. The CW breakup commences at $z\approx 150$.

Figure 14. An example of the transverse profile of a modulationally stable CW. The analytical profile is produced by Eqs. (29) and (31), with $W_0=g=1$ and $k=1.3$. Its numerically found counterpart cannot exactly reproduce the cusp of the analytical solution, determined by the jump condition (21).

Figure 14. An example of the transverse profile of a modulationally stable CW. The analytical profile is produced by Eqs. (29) and (31), with $W_0=g=1$ and $k=1.3$. Its numerically found counterpart cannot exactly reproduce the cusp of the analytical solution, determined by the jump condition (21).

First, both the numerical computation of the stability eigenvalues, as per the BdG system of Eqs. (9) and (), and direct simulations of the perturbed evolution easily demonstrate that, in terms of the 1D reduction of Eq. (5) (excluding coordinate y), all pinned states of both types, (24) and (29), are completely stable. However, in the framework of the full 2D equation (5), all the CW states based on the 1D solitons of type (24) are subject to MI, see a typical example in Figure 13. In particular, the direct test of the perturbed propagation, displayed in Figure 13(c), was performed by simulating Eq. (20) with the initial condition

$$u(x,y;z=0)=U\left(x\right)\left[1+2\epsilon cos\left(py\right)\right],$$

(33)

where $U\left(x\right)$ is the stationary profile, p and $\epsilon$ being the wavenumber and amplitude of the modulational perturbation (cf. Eq. (8)).

In view of the full instability of the CWs with the stationary profile (24), we do not consider this case in further detail, and focus on the CWs with the transverse profile (29), which turn out to be completely stable. First, an example of this profile for $W_0=1$, $g=1$ and $k=1.3$ is plotted in Figure 14.

Further, a family of such profiles, with the same parameters $W_0=1$, $g=1$ and varying values of k, is displayed in Figure 15, by means of the respective dependence of the 1D power P on the propagation constant k. In fact, the power of the analytical solution diverges in the limit of $k\to 3/\left(16g\right)$, which represents the left edge of the existence region in Figure 15.

Note that the curve $P\left(k\right)$ in Figure 15 obeys the anti-Vakhitov-Kolokolov(anti-VK) criterion, $dP/dk<0$, which is known as the necessary stability condition for self-trapped states dominated by self-repulsive nonlinearity [43], which is the case here, because the pinned solitons of type (29) do not exist in the absence of the self-defocusing quintic term in Eq. (5) (the VK criterion proper, $dP/dk>0$, is the necessary stability condition for solitons supported by self-attractive nonlinearity [44,45])

An illustration of the stability of the CW with the transverse profile provided by Eqs. (29) and (31), with $W_0=1.5$, $g=3$, $k=0.5$, is provided by Figure 16, in which the simulated propagation of the CW state, including initial perturbation (33) with $\epsilon =0.001$ and $p=0.3$, is displayed. To additionally test the CW robustness, a similar simulation is displayed in Figure 17 for the perturbation with a much larger amplitude, $\epsilon =0.1$. It is clearly observed that the CW is not destroyed even by strong perturbations.

5. Conclusion

The objective of this work is to consider the dynamics of Q1D (quasi-one-dimensional) CW (continuous-wave) states and related patterns (such as chains of 2D solitons) in the optical medium with the CQ (cubic-quintic) nonlinearity and the trapping potential trough. In the model with the usual smooth trapping trough, we have briefly reproduced the known MI (modulational instability) of CW states with the transverse structure corresponding to the GS (ground state) of the trapping potential, and produced new results for the CW states with the transverse DM (dipole-mode) structure, which corresponds to the first excited state of the potential trap. In particular, it is found that the CWs of both types, GS and DM, are “practically stable" near edges of their existence domains. Then, stable periodic chains of 2D solitons with alternating signs, maintained by the potential trough, have been found as the stationary manifestation of the MI, and their dynamics, initiated by the application of localized kicks, was explored. An essentially novel finding is represented by two families of exact CW solutions in the model with the delta-functional trapping potential, a remarkable result being the complete stability of one family, which does not exist without the trapping potential.

As an extension of the analysis, it may be interesting to consider chains of 2D solitons with embedded vorticities, identical or alternating (in the latter case, these are vortex-antivortex chains), trapped in the Q1D potential. It is known that individual vortex solitons with unitary and multiple topological charges have their stability regions in the free-space NLS equation with the CQ nonlinearity [19,20], therefore the trapped chain may be stable too. As a dynamical problem, one may address collisions between counterpropagating 2D solitons moving in the potential trough.