Discrete and Semi-discrete Multidimensional Solitons and Vortices: Established Results and Novel Findings

1. Introduction

1.1. Discrete Nonlinear Schrödinger (DNLS) Equations

1.1.1. The Basic Equation

Commonly adopted models of physical media are based on linear and nonlinear partial differential equations, such as the Gross-Pitaevskii (GP) equations for the mean-field wave function $\psi \left(x,y,z;t\right)$ in atomic Bose-Einstein condensates (BECs) [1], and nonlinear Schrödinger (NLS) equations for the amplitude of the electromagnetic field in optical waveguides [2,3]. The scaled form of the three-dimensional (3D) GP/NLS equation is

$$i{\psi }_t=-(1/2){\nabla }^2\psi +\sigma {\left|\psi \right|}^2\psi +U\left(x,y,z\right)\psi ,$$

(1)

where $\sigma =+1$ and $-1$ correspond to the self-defocusing and focusing signs of the cubic nonlinearity, and $U\left(x,y,z\right)$ is an external potential. In the application to optics, time t is replaced, as the evolution variable, by the propagation distance z, while the original coordinate z is then replaced by the temporal one, $\tau =t-z/V_{\mathrm{gr}}$, where $V_{\mathrm{gr}}$ is the group velocity of the carrier wave [2]. In optics, the effective potential may be two-dimensional (2D), $-U(x,y)$, which represents a local variation of the refractive index in the transverse plane.

In many cases the potential is spatially periodic, such as one induced by optical lattices (OLs) in BEC [4,5], or by photonic crystals which steer the propagation of light in optics [6],

$$U_{\mathrm{latt}}\left(x,y,z\right)=-\epsilon \left[cos(2\pi L/x)+cos(2\pi L/y)+cos(2\pi L/z)\right],$$

(2)

as well as its 2D and 1D reductions. A deep lattice potential, which corresponds to large amplitude $\epsilon$ in Eq. (2), splits the continuous wave function into a set of “droplets” trapped in local potential wells, which are linearly coupled by tunneling. Accordingly, in the framework of the tight-binding approximation, the NLS equation is replaced by a discrete NLS (DNLS) equation with the linear coupling between adjacent sites of the discrete lattice (nearest neighbors). This equation was derived, in the 1D form, for arrays of optical fibers [7,8,9,10] and plasmonic nanowires [11], as well as for BEC loaded in a deep OL potential [12]:

$$\begin{array}{c}i{\dot{\psi }}_{l,m,n}=-(1/2)\left[\left(u_{l+1,m,n}+{\psi }_{l-1,m,n}-2{\psi }_{l,m,n}\right)+\left({\psi }_{l,m+1,n}+{\psi }_{l,m-1,n}-2{\psi }_{l,m,n}\right)\right.\\ \left.+\left({\psi }_{l,m,n+1}+{\psi }_{l,m,n-1}-2{\psi }_{l,m,n}\right)\right]+\sigma {\left|{\psi }_{l,m,n}\right|}^2{\psi }_{l,m.n}+V_{l,m,n}{\psi }_{l,m,n},\end{array}$$

(3)

where the overdot stands for $d/dt$, and the set of integer indices, $\left(l,m,n\right)$, replaces continuous coordinates $\left(x,y,z\right)$ in Eq. (1). Here, potential $V_{l,m,n}$ is a possible smooth addition to the deep lattice potential which imposes the discretization. The full DNLS equation (3) is often reduced to 2D and 1D forms. 1D lattices are sometimes built in the form of zigzag chains, making it relevant to add couplings between the next-nearest neighbors [13,14]. 2D lattices with similar additional couplings were elaborated too [15].

Unlike the continuous NLS equation (1), the nonlinearity signs $\sigma =\pm 1$ in Eq. (3) are actually equivalent to each other. Indeed, the staggering transformation of the discrete wave function,

$${\psi }_{l,m,n}\left(t\right)\equiv {(-1)}^{l+m+n}exp\left(-6it\right){\tilde{\psi }}_{l,m,n}^{*}\left(t\right)$$

(4)

where * stands for the complex-conjugate, transforms Eq. (3) with $\sigma =-1$ into the same equation for ${\tilde{\psi }}_{l,m,n}$ with $\sigma =+1$. For the 2D and 1D DNLS equations, $exp\left(-6it\right)$ in the corresponding versions of staggering transformation (4) is replaced, respectively, by $exp\left(-4it\right)$ and $exp\left(-2it\right)$.

1.1.2. Extended Equations

The DNLS equation and its extensions, such as systems of coupled DNLS equations [16,17], constitute a class of models with a large number of physical realizations [10,18,20,21,22]. They have also drawn much interest as subjects of mathematical studies [23]. One of incentives for this interest is the fact that the discreteness arrests the development of the critical and supercritical collapse, which is driven by the self-focusing nonlinear term with $\sigma =-1$ in the 2D and 3D continuous NLS equations (1), respectively. The collapse leads to the emergence of singular solutions in the form of infinitely tall peaks, after a finite evolution time [3]. Naturally, the discreteness causes the arrest of the collapse, replacing it by a quasi-collapse [24], when the width of the shrinking peak becomes comparable to the spacing of the DNLS lattice.

The possibility of the collapse destabilizes formal 2D and 3D soliton solutions which are produced by Eq. (1), therefore a challenging problem is prediction and experimental realization of physical settings than make it possible to produce stable multidimensional solitons [25,26]. Thus, the discreteness provides a general method for the stabilization of 2D and 3D solitons.

The Gross-Pitaevskii (GP) Equations Amended by Effects of Quantum Fluctuations

Another promising possibility for the suppression of the collapse is offered by binary BEC, in which the cubic inter-component attraction creates 3D soliton-like states in the form of “quantum droplets" (QDs), while the development of the supercritical collapse is arrested by the self-repulsive quartic term that takes into account a correction to the mean-field nonlinearity produced by quantum fluctuations (known as the celebratedLee-Huang-Yang effect [28]). For the symmetric binary condensate, with identical wave functions $\psi$ of its components, the accordingly amended scaled GP equations (in the absence of the trapping potential) was derived by Petrov [29]

$$i{\psi }_t=-(1/2){\nabla }^2\psi -{\left|\psi \right|}^2\psi +{\left|\psi \right|}^3\psi .$$

(5)

Surprisingly quickly, the QD modes predicted by Eq. (5) have been created experimentally, in the quasi-2D [30,31] and full 3D [32] forms. The reduction of the spatial dimension to 2D and 1D replaces Eq. (5) by the following GP equations, respectively [33]:

$$\begin{array}{ccc}\hfill 2\mathrm{D}:& & i{\psi }_t=-(1/2){\nabla }^2\psi +{\left|\psi \right|}^2{ln\left(\right|\psi |}^2)·\psi ,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill 1\mathrm{D}:& & i{\psi }_t=-(1/2){\psi }_{xx}+{\left|\psi \right|}^2\psi -\left|\psi \right|\psi .\hfill \end{array}$$

(7)

Note, in particular, that in the 1D equation (7) the quantum correction is represented by the attractive term $-\left|\psi \right|\psi$, on the contrary to its repulsive counterpart $+{\left|\psi \right|}^3\psi$ in the 3D equation (5). For this reason, the usual mean-field cubic term is taken in Eq. (7) with the self-repulsion sign, to make it possible to study effects of competition of the quadratic self-attraction and cubic repulsion [34]. A semi-discrete version of Eq. (7) is considered below, see Eq. (66).

The Ablowitz-Ladik (AL) and Salerno-Model (SM) Equations

The 1D continuous NLS equation without an external potential is integrable by means of the inverse-scattering transform, with either sign of the nonlinearity, $\sigma =\pm 1$ [35,36,37,38], although 2D and 3D extensions of the NLS equation are nonintegrable. The straightforward discretization destroys the integrability of the 1D NLS equation [39,40]. Nevertheless, the NLS equation admits a specially designed 1D discretization, which leads to an integrable discrete model, viz., the Ablowitz-Ladik (AL) equation [41]:

$$i{\dot{\psi }}_n=-\left({\psi }_{n+1}+{\psi }_{n-1}\right)\left(1+\mu {\left|{\psi }_n\right|}^2\right),$$

(8)

where positive and negative values of the nonlinearity coefficient, $\mu$, correspond to the self-focusing and defocusing, respectively. Integrable discrete equations, such as the AL one, are exceptional models which provide exact solutions for discrete solitons [42].

Equation (8) gives rise to an exact solution for solitons in the case of $\mu >0$. Setting $\mu \equiv +1$ by means of rescaling, the solution is

$${\psi }_n\left(t\right)=\left(sinh\beta \right)\mathrm{sech}\left[\beta (n-\xi (t\left)\right)\right]exp\left[i\alpha \left(n-\xi \left(t\right)\right)-i\phi \left(t\right)\right],$$

(9)

where $\beta$ and $\alpha$ are arbitrary real parameters that determine the soliton’s amplitude, $A\equiv sinh\beta$, its velocity, $V\equiv \dot{\xi }=2{\beta }^{-1}\left(sinh\beta \right)sin\alpha$, and overall frequency $\Omega \equiv \dot{\phi }=-2\left[\left(cosh\beta \right)cos\alpha +(\alpha /\beta )\left(sinh\beta \right)sin\alpha \right]$. The existence of exact solutions for traveling solitons in the discrete system is a highly nontrivial property of the AL equation, which follows from its integrability. Generically, motion of a discrete soliton through a lattice is braked by emission of radiation, even if this effect may seem very weak in direct simulations [65]. Another integrable discrete model which admits exact solutions for moving solitons is the Toda-lattice equation for real coordinates $x_n\left(t\right)$ of particles with unit mass and exponential potential of interaction between adjacent ones [66]:

$${\ddot{x}}_n+exp\left(-\left(x_{n+1}-x_n\right)\right)-exp\left(-\left(x_n-x_{n-1}\right)\right)=0.$$

(10)

Considerable interest was also drawn to the nonintegrable combination of the AL and DNLS equations, in the form of the Salerno model (SM) [43], with an additional onsite cubic term, different from the intersite one in Eq. (8):

$$i{\dot{\psi }}_n=-\left({\psi }_{n+1}+{\psi }_{n-1}\right)\left(1+\mu {\left|{\psi }_n\right|}^2\right)-2{\left|{\psi }_n\right|}^2{\psi }_n.$$

(11)

Here, the signs and magnitude of the onsite nonlinearity coefficient are fixed by means of the staggering transformation (4) and rescaling. The SM finds a physical realization in the context of the Bose-Hubbard model, which represents the BEC loaded in a deep OL, taking into regard the nonlinearity of the intersite coupling [44].

The AL and SM equations (8) and (11) conserve the total norm N, whose definition is different from the straightforward one for the DNLS equation, given below by Eq. (27); namely,

$$N_{\mathrm{AL},\mathrm{SM}}=(1/\mu )\sum _{n}ln\left|1+\mu |{\psi }_n{|}^2\right|$$

(12)

[41,107]. The Hamiltonians of the AL and SM equations, which are dynamical invariants too, are

$$H_{\mathrm{AL}}=-\sum _n\left({\psi }_n{\psi }_{n+1}^{*}+{\psi }_{n+1}{\psi }_n^{*}\right),$$

(13)

$$H_{\mathrm{SM}}=-\sum _n\left[\left({\psi }_n{\psi }_{n+1}^{*}+{\psi }_{n+1}{\psi }_n^{*}\right)+(2/\mu ){\left|{\psi }_n\right|}^2\right]+(2/\mu )N_{\mathrm{SM}}.$$

(14)

In particular, the ostensible “simplicity” of Hamiltonian (13) is related to the complexity of the respective Poisson brackets (symplectic structure), which determine the evolution equations for ${\psi }_n$ as $d{\psi }_n/dt=\left\{H,{\psi }_n\right\}$. For the AL and SM models, the Poisson brackets, written for a pair of arbitrary functions $B\left({\psi }_n,{\psi }_n^{*}\right)$ and $C\left({\psi }_n,{\psi }_n^{*}\right)$, are [43,44]

$$\left\{B,C\right\}=i\sum _n\left(\frac{\partial B}{\partial {\psi }_n}\frac{\partial C}{\partial {\psi }_n^{*}}-\frac{\partial B}{\partial {\psi }_n^{*}}\frac{\partial C}{\partial {\psi }_n}\right)\left(1+\mu {\left|{\psi }_n\right|}^2\right).$$

(15)

It is also relevant to consider the continuum limit of the SM, which is introduced by approximating the intersite combination of the discrete fields by a truncated Taylor’s expansion,

$${\psi }_n\left(t\right)\equiv e^{2it}\Psi (x,t),\Psi \left(x=n\pm 1,t\right)\approx \Psi \left(x=n\right)\pm {\Psi }_x{|}_{x=n}+(1/2){\Psi }_{xx}{|}_{x=n},$$

(16)

where $\Psi \left(x\right)$ is considered as a function of continuous coordinate x, whose integer values coincide with the discrete lattice coordinate n. The substitution of approximation (16) in Eq. (11) leads to a generalized (nonintegrable) NLS equation [55]

$$i{\Psi }_t=-\left(1+\mu {\left|\Psi \right|}^2\right){\Psi }_{xx}-2\left(1+\mu \right){\left|\Psi \right|}^2\Psi .$$

(17)

Equation (17) conserves the total norm and Hamiltonian, which are continuum counterparts of expressions (12) and (14):

$$\begin{array}{ccc}\hfill {\left(N_{\mathrm{AL}}\right)}_{\mathrm{cont}}& =& \frac{1}{\mu }{\int }_{-\infty }^{+\infty }dxln\left|1+{\mu |\Psi |}^2\right|,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\left(H_{\mathrm{SM}}\right)}_{\mathrm{cont}}& =& {\int }_{-\infty }^{+\infty }dx\left[{\left|{\Psi }_x\right|}^2-2\left(\frac{1}{\mu }+1\right){|\Psi |}^2\right]+\frac{2}{\mu }{\left(N_{\mathrm{AL}}\right)}_{\mathrm{cont}}\hfill \end{array}$$

(19)

It is relevant to mention that the general approximation opposite to the continuum limit is the anti-continuum limit [18,56]. This approach starts with the limit form of the DNLS equation, in which the linear couplings between adjacent sites are dropped. Then, one can try to construct various states, including solitons, by introducing an input composed of simple solutions of the single-site equations corresponding to Eq. (3), viz., $i{\dot{\psi }}_{l,m,n}=\sigma {\left|{\psi }_{l,m,n}\right|}^2{\psi }_{l,m.n}+V_{l,m,n}{\psi }_{l,m,n}$, at a finite set of sites, and keeping the zero solution at all others. The single-site “simple solutions" are

$${\psi }_{l,m,n}=a_{l,m,n}exp\left(-i{\left|a_{l,m,n}\right|}^{2}t-i{V}_{l,m,n}t\right),$$

(20)

where $a_{l,m,n}$ is an arbitrary set of complex amplitudes. Next, one reintroduces weak linear intersite couplings and attempts to identify nontrivial solutions that may thus appear from the finite-set input composed of the single-site solutions (20).

Self-Trapping in Lattices with the Self-Repulsion Strength Growing from the Center to Periphery

DNLS equations with the onsite self-repulsive nonlinearity, corresponding to $\sigma >0$ in Eq. (3), may support discrete-soliton (self-trapped) states without the resort to the staggering transform (4) if the local self-repulsion strength is made a function of the lattice coordinates, growing fast enough from the center to periphery. Originally, this option was elaborated in the framework of the continuum NLS and GP equations in the space of dimension D, with the local self-defocusing (repulsion) coefficients growing at $r\to \infty$ (r is the radial coordinate) faster than $r^D$ [45]. In terms of the 1D and 2D DNLS equations, similar settings were introduced in Refs. [46] and [47], with the site-dependent self-attraction coefficients, ${\left({\sigma }_n\right)}_{1\mathrm{D}}={\sigma }_{0}exp\left(\alpha \left|n\right|\right)$ and ${\left({\sigma }_{m,n}\right)}_{2\mathrm{D}}={\sigma }_{0}exp\left(\alpha \left(\left|m\right|+\left|n\right|\right)\right)$, respectively, with positive constants ${\sigma }_0$ and $\alpha$. In the 2D model, solutions were constructed, and their stability analyzed, for fundamental, dipole, quadrupole, and vortical discrete solitons [45].

DNLS Equations with Long-Range Dipole-Dipole and Quadrupole-Quadrupole Intersite Interactions

It is well known that atomic BEC formed of ultracold atoms carrying permanent magnetic moments feature specific dynamical effects due to the long-range interactions between atomic moments [48]. This fact suggests to combine the dipole-dipole interactions and a deep OL potential, thus introducing DNLS equations with the nonlocal (long-range) coupling between lattice sites. In the 2D setting, this model gives rise to different forms of the DNLS equations. The simplest setup corresponds to the case when the atomic moments are polarized by external dc magnetic field perpendicular to the system’s plane. In this case, the dipole-dipole interactions amount to the isotropic nonlocal repulsion, accounted for by the respective interaction coefficient $\Gamma >0$ [49]

$$\begin{array}{ccc}\hfill i\frac{\partial {\psi }_{m,n}}{\partial t}& =& -\frac{1}{2}\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{m,n}\right)+\left[\sigma {\left|{\psi }_{m,n}\right|}^2\right.\hfill \\ & & \left.+\Gamma \sum _{\left\{m^{\prime },n^{\prime }\right\}\ne \left\{m,n\right\}}{\left|{\left(m-m^{\prime }\right)}^2+{\left(n-n^{\prime }\right)}^2\right|}^{-3/2}{\left|{\psi }_{m^{\prime },n^{\prime }}\right|}^2\right]{\psi }_{m,n},\hfill \end{array}$$

(21)

where $\sigma$ is the same coefficient of the onsite self-interaction as in Eq. (3). A more sophisticated setup corresponds to the atomic magnetic moments polarized parallel to the system’s plane. In the former case, the nonlocal term term in the respective DNLS equation is anisotropic, being attractive in one in-plane direction and repulsive in the other, cf. Ref. [51]:

$$\begin{array}{c}i\frac{\partial {\psi }_{m,n}}{\partial t}=-\frac{1}{2}\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{m,n}\right)+\left[\sigma {\left|{\psi }_{m,n}\right|}^2\right.\\ \left.+\Gamma \sum _{\left\{m^{\prime },n^{\prime }\right\}\ne \left\{m,n\right\}}\frac{{\left(n-n^{\prime }\right)}^2-2{\left(m-m^{\prime }\right)}^2}{{\left[{\left(m-m^{\prime }\right)}^2+{\left(n-n^{\prime }\right)}^2\right]}^{5/2}}{\left|{\psi }_{m^{\prime },n^{\prime }}\right|}^2\right]{\psi }_{m,n}.\end{array}$$

(22)

The analysis reported in Ref. [49] demonstrates that the nonlocal repulsion, added to Eq. (21), helps to stabilize discrete solitons with embedded vorticity. Solutions of Eq. (22) for anisotropic vortex solitons can be found too, but they are completely unstable [49].

A 2D DNLS model which combines the local onsite nonlinearity and long-range interaction between particles carrying permanent quadrupole electric moments was elaborated in Ref. [52]:

$$\begin{array}{c}i\frac{\partial {\psi }_{m,n}}{\partial t}=-\frac{1}{2}\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{mn}\right)+\left[\sigma |{\psi }_{mn}{|}^2\right.\\ \left.+\Gamma \sum _{\left\{m^{\prime },n^{\prime }\right\}\ne \left\{m,n\right\}}\frac{{(n-n^{\prime })}^2-4{(m-m^{\prime })}^2}{{[{(m-m^{\prime })}^2+{(n-n^{\prime })}^2]}^{7/2}}{\left|{\psi }_{m^{\prime }{n}^{\prime }}\right|}^2\right]{\psi }_{m,n},\end{array}$$

(23)

cf. Eq. (21). This model also gives rise to families of stable 2D discrete solitons [52].

The 2D Discrete Second-Harmonic-Generating (${\chi }^{\left(2\right)}$) system

The quadratic (alias ${\chi }^{\left(2\right)}$) nonlinearity is a fundamentally important effect which gives rise to coherent generation of the second harmonic in optics. In terms of the 2D spatial-domain propagation in a continuum material, the standard ${\chi }^{\left(2\right)}$ system for amplitudes $\psi \left(x,y,z\right)$ and $\varphi \left(x,y,z\right)$ of the fundamental-frequency (FF) and second-harmonic (SH) waves is [53]

$$\begin{array}{ccc}\hfill i{\psi }_z& =& -(1/2){\nabla }^2\psi -{\psi }^{*}\varphi ,\hfill \\ \hfill 2i{\varphi }_z& =& -(1/2){\nabla }^2\varphi -Q\varphi -\frac{1}{2}{\psi }^2,\hfill \end{array}$$

(24)

where z is the propagation distance, the paraxial-diffraction operator $(1/2){\nabla }^2$ acts on the transverse coordinates $\left(x,y\right)$, Q is a real mismatch parameter, and * stands for the complex conjugate. The discretized version of Eqs. (24), which represents, in the tightly-binding approximation, light propagation in a photonic crystal made of the ${\chi }^{\left(2\right)}$ material, is [54]

$$\begin{array}{ccc}\hfill i\frac{d{\psi }_{m,n}}{dz}& =& -C_1\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{mn}\right)-{\psi }_{m,n}^{*}{\varphi }_{m,n},\hfill \\ \hfill 2i\frac{d{\varphi }_{m,n}}{dz}& =& -\frac{C_2}{2}\left({\varphi }_{m+1,n}+{\varphi }_{m-1,n}+{\varphi }_{m,n+1}+{\varphi }_{m,n-1}-4{\varphi }_{mn}\right)-Q{\varphi }_{m,n}-\frac{1}{2}{\psi }_{m,n}^2,\hfill \end{array}$$

(25)

where $C_1$ and $C_2$ are effective lattice-coupling constants for the FF and SH waves. The role of the conserved norm of the discrete ${\chi }^{\left(2\right)}$ system is played by the Manley-Rowe invariant, i.e., the total optical power, $I={\sum }_{m,n}\left({\left|{\psi }_{m,n}\right|}^2+4{\left|{\varphi }_{m,n}\right|}^2\right)$.

An essential property of 2D discrete solitons produced by Eqs. (25) is their mobility [54]. In this connection, it is relevant to mention that, while the development of the quasi-collapse in the 2D discrete NLS equation with the cubic self-attraction is arrested by the underlying lattice structure, the quasi-collapse strongly pins the 2D solitons to the same structure, and thus makes them immobile. On the other hand, the ${\chi }^{\left(2\right)}$ nonlinearity does not give rise to the collapse in the 2D (and 3D) space, therefore 2D ${\chi }^{\left(2\right)}$ solitons do not demonstrate a trend for strong pinning, remaining effectively mobile robust localized modes [54].

1.2. One-Dimensional DNLS Solitons

1.2.1. Fundamental Solitons

In the 1D setting, the model of basic interest is the DNLS equation with self-attraction, which corresponds to $\sigma =-1$ in the 1D version of Eq. (3), without the external potential ($V_{l,m,n}=0$):

$$i{\dot{\psi }}_n=-(1/2)\left({\psi }_{n+1}+{\psi }_{n-1}-2{\psi }_n\right)-{\left|{\psi }_n\right|}^2{\psi }_n.$$

(26)

This equation conserves two dynamical invariants, viz., the total norm,

$$N=\sum _{n=-\infty }^{+\infty }{\left|{\psi }_n\right|}^2,$$

(27)

and Hamiltonian (energy),

$$H=\sum _{n=-\infty }^{+\infty }\left[(1/2){\left|{\psi }_n-{\psi }_{n-1}\right|}^2-(1/4){\left|{\psi }_n\right|}^4\right].$$

(28)

Stationary solutions to Eq. (26) with real frequency $\omega$ are looked for as

$${\psi }_n\left(t\right)=e^{-i\omega t}{u}_n,$$

(29)

with real amplitudes $u_n$ satisfying the discrete equation,

$$\omega u_n=-(1/2)\left(u_{n+1}+u_{n-1}-2u_n\right)-u_n^3.$$

(30)

Note that Eq. (30) can be derived by varying its Lagrangian,

$$L=\sum _{n=-\infty }^{+\infty }\left\{(1/4)\left[{\left(u_n-u_{n-1}\right)}^2-u_n^4\right]-\omega u_n^2\right\},$$

(31)

with respect to the discrete real field $u_n$.

A fundamental property of the DNLS equation (26) with the self-attractive onsite nonlinearity is the modulational instability (MI) of its spatially homogeneous continuous-wave (CW) state [57], ${\psi }_n=aexp\left(i{a}^{2}t\right)$, with an arbitrary amplitude a [cf. Eq. (29)]. MI breaks the CW state into a chain of discrete solitons [10]. Analytical solutions for these solitons are not available, as the DNLS equation is not integrable. The solitons can be readily found in a numerical form, and studied in the framework of the variational approximation (VA) [23]. The VA is based on a particular ansatz, i.e. an analytical expression which aims to approximate the solution [58]. The only discrete ansatz for which analytical calculations are feasible is represented by the exponential function [59,60,61,62], namely,

$${\left(u_n\right)}_{\mathrm{onsite}}=Aexp\left(-a\left|n\right|\right),$$

(32)

with $a>0$. The corresponding norm, calculated as per Eq. (27), is

$$N_{\mathrm{ansatz}}=A^{2}cotha.$$

Note that ansatz (32) is appropriate for strongly and moderately discrete solitons, as shown in Figure 1, but not for broad (quasi-continuum) ones, which may be approximated by the commonly known soliton solution of the NLS equation (the 1D version of (1) with $U=0$),

$$\psi \left(x,t\right)=\eta \mathrm{sech}\left(\eta \left(x-\xi \right)\right)exp\left(i{\eta }^{2}t\right),$$

(33)

with width ${\eta }^{-1}$ which must be large in comparison with the discreteness spacing, ${\eta }^{-1}\gg 1$, and central coordinate $\xi$.

The substitution of ansatz (32) in Lagrangian (31) produces the corresponding VA Lagrangian:

$$L_{\mathrm{VA}}=(A^2/2)coth(a/2)-(A^4/4)coth\left(2a\right)-\omega A^{2}cotha.$$

(34)

Then, for given $\omega <0$ (solitons do not exist for $\omega >0$), the squared amplitude, $A^2$, and inverse width, a, of the discrete soliton are predicted by the Euler-Lagrange equations,

$$\frac{\partial L_{\mathrm{eff}}}{\partial \left(A^2\right)}=\frac{\partial L_{\mathrm{eff}}}{\partial a}=0.$$

(35)

This corresponding system of algebraic equations for $A^2$ and a can be easily solved numerically. The VA produces an accurate predictions for the solitons, as shown in Figure 1 and Ref. [63].

Rigorous justification of the VA was elaborated in Ref. [64]. Furthermore, the VA and a full numerical solution of Eq. (26) demonstrate that the entire family of the discrete solitons is stable [23].

1.2.2. Higher-Order One-Dimensional Modes: Twisted Discrete Solitons and Bound States

In addition to the fundamental (single-peak) solitons outlined above, Eq. (30) admits stable second-order states in the form of twisted modes, which are subject to the antisymmetry condition, $u_n=-u_{1-n}$ [67]. Such states exist and are stable only in a strongly discrete form, vanishing in the continuum limit. In particular, the above-mentioned anti-continuum approximation is appropriate for the construction of the twisted modes.

Stable 1D DNLS solitons may form bound states, which also represent higher-order modes of the DNLS equation. They are stable in the out-of-phase form, i.e., with opposite signs of the constituent solitons [68,69], which resembles the structure of the twisted modes (however, the tight antisymmetric structure of the twisted modes cannot be considered as a bound state of fundamental solitons), the same being true for 2D DNLS solitons [70]. Note that stationary bound states of fundamental solitons do not exist in the continuum limit.

1.2.3. 1D Solitons in the Salerno Model (SM)

The substitution of the usual ansatz (29) in Eq. (11) produces a stationary discrete equation for real $u_n$:

$$\omega u_n=-\left(u_{n+1}+u_{n-1}\right)\left(1+\mu u_n^2\right)-2u_n^3,$$

(36)

cf. (30). Discrete solitons produced by the SM equation (11) with $\mu >0$, i.e. with noncompeting intersite and onsite self-focusing nonlinearities, were investigated by means of numerical methods [107,108,109]. It was found that the SM gives rise to the 1D solitons at all positive values of $\mu$.

Another option is to consider the SM with $\mu <0$, which features competing nonlinearities, as the intersite cubic term, with coefficient $\mu <0$ in Eq. (11), which accounts for the nonlinear coupling between adjacent sites of the lattice, and the onsite cubic term in Eq. (11) represent, respectively, repulsive and attractive nonlinear interactions. This version of the SM gives rise to families of discrete solitons, in the usual form (29), with $\omega <0$ and real $u_n$, of two different types. One family represents ordinary discrete solitons, similar to those generated by the DNLS equation. Another family represents cuspons, featuring higher curvature of their profile at the center. A small subfamily of ordinary solitons produced by the SM with the competing nonlinearities is unstable, while all cuspons are stable.

As mentioned above, antisymmetric bound states of DNLS solitons are stable, while symmetric bound states are unstable [68,69]. The same is true for bound states of ordinary discrete solitons in the SM [55]. However, in the framework of the SM with the competing nonlinearities, the situation is exactly opposite for the cuspons: their symmetric and antisymmetric bound states are stable and unstable, respectively [55].

1.3. The Subject and Structure of the Present Article

The above-mentioned reviews [10,18,20,21,22] and [23] produce a comprehensive survey of theoretical and experimental results for discrete solitons in various 1D systems. The objective of this article is to produce a relatively brief summary of results for multidimensional (chiefly, two-dimensional) discrete and semi-discrete solitons, which were considered in less detail in previous reviews and, on the other hand, draw growing interest in the context of the current work with 2D and 3D solitons in diverse physical contexts [25,26]. In this context, the presence of the two or three coordinates makes it also possible to define semi-discrete states as ones which are discrete in one direction and continuous in the perpendicular one [71,72,73,74,77]. The article chiefly represents theoretical results, but some experimental findings for quasi-discrete 2D solitons in photonic lattices [79,80] are included too.

The review presented below does not claim to be comprehensive. The comprises results that are produced by conservative models of the DNLS types (including, in particular, the 2D SM). Discrete models of other types – in particular, those similar to the Toda lattice, see Eq. (10), Fermi-Pasta-Ulam-Tsingou lattices [81], and Frenkel-Kontorova systems [82] – are not considered here. Dissipative systems are not considered either, except for a 2D model with the parity-time ($\mathcal{PT}$) symmetry [83], see Eq. (70) below.

The rest of the article is arranged as follows. Basic results for fundamental (zero-vorticity) and vortex solitons, as well as bound states of such solitons, produced by the 2D DNLS equation and its generalizations, are summarized in Section II, which is followed, in Section III, by brief consideration of fundamental and vortex solitons in the 2D SM (Salerno model). Section IV addresses discrete solitons of the semi-vortex and mixed-mode types in the 2D spin-orbit-coupled (SOC) system of GP equations for a two-component BEC. Basic results for discrete self-trapped modes produced by 3D DNLS equations, including fundamental and vortex solitons, along with skyrmions, are presented in Section V. The findings for 2D semi-discrete systems, again including fundamental and vortex solitons, supported by combined quadratic-cubic and cubic-quintic nonlinearities (that are relevant for BEC and optics, respectively), are summarized in Section VI. This section also addresses transverse mobility of confined spatiotemporal modes in an array of optical fibers with the intrinsic cubic self-focusing (Kerr nonlinearity). Fundamental and vortex solitons produced by a $\mathcal{PT}$-symmetric discrete 2D coupler with the cubic nonlinearity are considered in Section VII. The article is concluded by Section VIII, which, in particular, suggests directions for the further work in this area, and mentions particular topics which are not included in the present review.

2. Two-Dimensional (2D) Nonlinear-Schrödinger Lattices: Fundamental and Vortex Solitons, and Their Bound States

2.1. Vortex Solitons: Theoretical and Experimental Results

The basic 2D cubic DNLS equation is the 2D version of Eq. (3) with the self-attraction ($\sigma =-1$) and without the external potential:

$$i{\dot{\psi }}_{m,n}=-(1/2)\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{m,n}\right)-{\left|{\psi }_{m,n}\right|}^2{\psi }_{m.n},$$

(37)

cf. Eq. (26). The substitution of ${\psi }_{m,n}=exp\left(-i\omega t\right)u_{m,n}$ in Eq. (37) produces the stationary equation, where the stationary discrete wave function, $u_{m,n}$, may be complex:

$$\omega u_{m,n}=-(1/2)\left(u_{m+1,n}+u_{m-1,n}+u_{m,n+1}+u_{m,n-1}-4u_{m,n}\right)-{\left|u_{m,n}\right|}^2{u}_{m,n},$$

(38)

cf. (30). Fundamental-soliton solutions to Eq. (38) can be predicted by means of VA [84,85], using an exponential ansatz, see (47) below (cf. Eq. (32) for the 1D soliton). More interesting in the 2D setting are discrete solitons with embedded vorticity, which were introduced in Ref. [86] (see also Ref. [87]). Vorticity, alias the topological charge, or winding number, is defined as $\Delta \phi /\left(2\pi \right)$, where $\Delta \phi$ is a total change of the phase of the complex discrete function $u_{m,n}$ along a contour surrounding the vortex’ pivot. Stability is an important issue for 2D discrete solitons, because, in the continuum limit, the 2D NLS equation gives rise to the well-known Townes solitons [88], which are unstable against the onset of the critical collapse [3]. In the same limit, the Townes solitons with embedded vorticity (vortex rings [89]) are subject to much stronger instability against spontaneous splitting of the ring in fragments [90].

The lattice structure of the DNLS equation provides for stabilization of both fundamental (zero-vorticity) and vortex solitons [86]. A typical example of a stable 2D vortex soliton with topological charge $S=1$ is displayed in Figure 2. 2D fundamental and vortex solitons, with topological charges $S=0$ and 1, are stable in regions $-\omega >|{\omega }_{\mathrm{cr}}^{(S=0)}|\approx 0.50$ and $-\omega >|{\omega }_{\mathrm{cr}}^{(S=1)}|\approx 1.23$, respectively [86], while the higher-order discrete vortices with charges $S=2$ and 4 are unstable, being replaced by stable modes in the form of quadrupoles and octupoles [91]. The vortex solitons with $S=3$ may be stable, but only in an extremely discrete form, viz., at $-\omega >|{\omega }_{\mathrm{cr}}^{(S=2)}|\approx 4.94$. In agreement with what is said above, these results imply that all the solitons are unstable in the continuum limit, corresponding, in the present notation, to $\omega \to 0$.

The experimentally relevant lattice structure may be anisotropic, with the linear combination $\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{m,n}\right)$ in Eq. (37) replaced by

$\left(\alpha \left({\psi }_{m+1,n}+{\psi }_{m-1,n}\right)+\left({\psi }_{m,n+1}+{\psi }_{m,n-1}\right)-2\left(1+\alpha \right){\psi }_{m,n}\right),$ with anisotropy parameter $\alpha \ne 1$. Effects of the anisotropy on the structure and stability of the fundamental and vortical discrete solitons were explored in Ref. [102].

The theoretically predicted 2D discrete solitons with vorticity $S=1$ were experimentally created in Refs. [79] and [80], using a photorefractive crystal. Unlike uniform media of this type, where delocalized (“dark") optical vortices were originally demonstrated [92,93], those works made use of a deep virtual photonic lattice as a quasi-discrete structure supporting the self-trapping of nonlinear modes in the optical field with extraordinary polarization (while the photonic lattice was induced as the interference pattern of quasi-linear beams in the ordinary polarization). Intensity distributions observed in vortex solitons of the onsite- and intersite-centered types (i.e., with the vortex’ pivot coinciding with a site of the underlying lattice, or set between the sites, respectively), are displayed in Figure 3.

Another interesting result demonstrated (and theoretically explained) in deep virtual photonic lattices is a possibility of periodic flipping of the topological charge of a vortex soliton initially created with topological charge $S=2$ [94]. Stable vortex solitons with $S=2$ were created using a hexagonal virtual photonic lattice [95].

2.2. Bound States of 2D Discrete Solitons and Solitary Vortices

As mentioned above, stable 2D discrete solitons may form stable bound states, composed of two or several items. Vortex solitons may form bound states as well. This possibility and stability of the resulting bound states are determined by an effective potential of interaction between two identcial discrete solitons with intrinsic vorticity S, which are separated by large distance L. The potential can be derived from an asymptotic expression for exponentially decaying tails of the soliton. In the quasi-continuum approximation, it is $u_{m,n}\sim {\left(m^2+n^2\right)}^{-1/4}exp\left(-\sqrt{-2\omega \left(m^2+n^2\right)}\right)$ (recall soliton solutions to Eq. (38) exist for $\omega <0$). Then, the overlap of the tail of each soliton with the central body of the other one gives rise to the following interaction potential:

$$U_{\mathrm{int}}\left(L\right)\approx \mathrm{const}·cos\delta ·{(-1)}^S{L}^{-1/2}exp\left(\sqrt{-2\omega }L\right),$$

(39)

with $\mathrm{const}>0$, where $\delta$ is the phase shift between the solitons. Thus, for the fundamental solitons with $S=0$, Eq. (39) predicts the attractive interaction between the in-phase solitons ($\delta =0$), and repulsion between out-of-phase ones ($\delta =\pi$). Accordingly, the interplay of the repulsive interaction with the effective Peierls-Nabarro potential, which is pinning the soliton to the underlying lattice [23], produces stable bound states of two or several mutually out-of-phase solitons, while the in-phase bound states are unstable. These predictions were confirmed by numerical results [70]. As an example, Figure 4 displays a numerically found stable bound state in the form of a string of three fundamental solitons with alternating phases.

For the pair of identical vortex solitons with $S=1$, Eq. (39) predicts the opposite result, viz., the repulsive interaction and stability of the ensuing bound states for in-phase vortices ($\delta =0$), and the attraction leading to instability of the bound state in the case of $\delta =\pi$. These predictions were also corroborated by numerical findings [70], see an example of a stable bound state of two identical vortex solitons in Figure 5.

2.3. 2D Discrete Solitons in Mini-Gaps of a Spatially Modulated Lattice

A specific class of self-trapped modes are gap solitons which may populate finite bandgaps in linear spectra of various nonlinear systems originating in optics and BEC [96,97,98]. While in most cases gap solitons are predicted theoretically and created experimentally in the context of matter waves [99] and optical pulses [100] in the continuum, they may naturally appear as discrete modes in mini-gaps, which are induced in the linear spectrum of lattice media by superimposed periodic spatial modulations (superlattice).

Such a 2D lattice model was introduced in Ref. [101], based on the following DNLS equation:

$$i{\dot{\psi }}_{m,n}+C_{m,n}\left({\psi }_{m+1,n}+{\psi }_{m-1,n}\right)+K_{m,n}\left({\psi }_{m,n+1}+{\psi }_{m,n-1}\right)+{\left|{\psi }_{m,n}\right|}^2{\psi }_{m,n}=0,$$

(40)

where the horizontal and vertical coupling constants are modulated as follows:

$$C_{m,n}=1+\Delta cos\left(Qm\right),K_{m,n}=1+\Delta cos\left(Qn\right),$$

(41)

cf. Eq. (3). The superlattice represented by Eqs. (40) and (41) can be created by means of the technique used for making OLs in experiments with BEC [101].

Looking for solutions in the usual form, ${\psi }_{m,n}\left(t\right)=exp\left(-i\omega t\right)u_{m,n}$, with real frequency (chemical potential) $\omega$, the numerical analysis produces the linear spectrum of this system, including the usual semi-infinite bandgap and a pair of additional narrow mini-gaps. Further, a family of fundamental 2D discrete solitons populating the mini-gaps was furnished by the numerical solution of the full nonlinear system, being stable in a small section of the mini-gap, as shown in Figure 5. The stable 2D soliton displayed in panel Figure 5(a) features a typical shape of gaps solitons, with a number of small satellite peaks surrounding the tall central one [97]. Bound states of two and four fundamental solitons were found too, featuring weak instability [101].

2.4. 2D Discrete Solitons in a Rotating Lattice

Dynamics of BEC loaded in OLs rotating at angular velocity $\Omega$, as well as the propagation of light in a twisted nonlinear photonic crystal with pitch $\Omega$, is modeled by the 2D version of Eq. (1) including the lattice potential, with depth $\epsilon$ and period $2\pi /k$, written in the rotating reference frame:

$$i{\psi }_t=-\left[(1/2){\nabla }^2+\Omega {\widehat{L}}_z\right]\psi -\epsilon \left[cos\left(kx\right)+cos\left(ky\right)\right]\psi +\sigma {\left|\psi \right|}^2\psi ,$$

(42)

where ${\widehat{L}}_z=i(x{\partial }_y-y{\partial }_x)\equiv i{\partial }_{\theta }$ is the operator of the z-component of the orbital momentum, $\theta$ being the angular coordinate in the $\left(x,y\right)$ plane. In the tight-binding approximation, Eq. (42) is replaced by the following variant of the DNLS equation [103]:

$$\begin{array}{c}i{\dot{\psi }}_{m,n}=-(C/2)\left\{\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{m,n}\right)\right.\\ \left.-i\Omega \left[m\left({\psi }_{m,n+1}-{\psi }_{m,n-1}\right)-n\left({\psi }_{m+1,n}-{\psi }_{m-1,n}\right)\right]\right\}+\sigma {\left|{\psi }_{m,n}\right|}^2{\psi }_{m,n},\end{array}$$

(43)

where C is the intersite coupling constant. In Ref. [103], stationary solutions to Eq. (43) were looked for in the usual form (29), fixing $\omega \equiv -1$ and varying C in (43) as a control parameter. Two species of localized states were thus constructed: off-axis fundamental discrete solitons, placed at distance R from the origin, and on-axis ($R=0$) vortex solitons, with topological numbers $S=1$ and 2. At a fixed value of rotation frequency $\Omega$, a stability interval for the fundamental soliton, $0<C<C_{max}^{\left(\mathrm{fund}\right)}\left(R\right)$, monotonously shrinks with the increase of R, i.e., most stable are the discrete fundamental solitons with the center placed at the rotation pivot. Vortices with $S=1$ are gradually destabilized with the increase of $\Omega$ (i.e., their stability interval, $0<C<C_{max}^{\left(\mathrm{vort}\right)}(\Omega )$, shrinks). On the contrary, a remarkable finding is that vortex solitons with $S=2$, which, as said above, are completely unstable in the usual DNLS equation with $\Omega =0$, are stabilized by the rotation, in an interval $0<C<C_{\mathrm{cr}}^{(S=2)}(\Omega )$, with $C_{\mathrm{cr}}^{(S=2)}(\Omega )$ growing as a function of $\Omega$. In particular, $C_{\mathrm{cr}}^{(S=2)}(\Omega )\approx 2.5\Omega$ at small $\Omega$ [103].

2.5. Spontaneous Symmetry Breaking of the 2D Discrete Solitons in Linearly-Coupled Lattices

A characteristic feature of many nonlinear dual-core systems, built of two identical linearly-coupled waveguides with intrinsic self-attractive nonlinearity, is a spontaneous-symmetry-breaking (SSB) bifurcation, which destabilizes the symmetric ground state, with equal components of the wave function in the coupled cores, and creates stable asymmetric states. The SSB bifurcation takes place at a critical value of the nonlinearity strength, the asymmetric state existing above this value [105]. A system of linearly-coupled DNLS equations is a basic model for SSB in discrete settings. Its 2D form is [17]

$$\begin{array}{ccc}\hfill i{\dot{\varphi }}_n& =& -(1/4)\left({\varphi }_{m+1,n}+{\varphi }_{m-1,n}+{\varphi }_{m,n+1}+{\varphi }_{m,n-1}-4{\varphi }_{m,n}\right)-{\left|{\varphi }_{m,n}\right|}^2{\varphi }_{m,n}-K{\psi }_{m,n},\hfill \\ \hfill i{\dot{\psi }}_n& =& -(1/4)\left({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1}-4{\psi }_{m,n}\right)-{\left|{\psi }_{m,n}\right|}^2{\psi }_{m,n}-K{\varphi }_{m,n},\hfill \end{array}$$

(44)

where ${\varphi }_{m,n}$ and ${\psi }_{m,n}$ are wave functions of the discrete coordinates m and n, and $K>0$ represents the linear coupling between the cores. Stationary states with frequency $\omega$ are looked for as $\left({\varphi }_{m,n},{\psi }_{m,n}\right)=exp\left(-i\omega t\right)\left(u_{m,n},v_{m.n}\right)$. Real stationary fields in the two components are characterized by their norms,

$$E_{u,v}=\sum _{m,n=-\infty }^{+\infty }\left(u_{m,n}^2,v_{m,n}^2\right),$$

(45)

and the asymmetry of the symmetry-broken state is determined by parameter

$$r=\left(E_u-E_v\right)/\left(E_u+E_v\right).$$

(46)

The system under the consideration can be analyzed by means of the VA, based on the 2D ansatz

$$\left(u_{m,n},v_{m,n}\right)=\left(A,B\right)exp\left[-a\left(\left|m\right|+\left|n\right|\right)\right],$$

(47)

with inverse width a and amplitudes, A and B, of the two components (cf. the 1D ansatz (32)). The SSB is represented by solutions with $A\ne B$. An example of a stable 2D discrete soliton is displayed in Figure 7(a), which corroborates accuracy of the VA. In Figure 7(b), the families of symmetric and asymmetric 2D discrete solitons is characterized by the dependence of asymmetry parameter r, defined as per Eq. (46), on the total norm, $E\equiv E_u+E_v$, see Eq. (45). Figure 7(b) demonstrates the SSB bifurcation of the subcritical type [106], with the two branches of broken-symmetry states originally going backward in the E direction, as unstable ones; they become stable after passing the turning point. Accordingly, Figure 7(b) demonstrates a bistability area, where symmetric and asymmetric states coexist as stable ones.

3. 2D Discrete Solitons in the Salerno Model (SM)

The 2D version of the SM was introduced in Ref. [111]:

$$\begin{array}{ccc}\hfill i{\dot{\psi }}_{n,m}& =& -\left[\left({\psi }_{n+1,m}+{\psi }_{n-1,m}\right)+\left({\psi }_{n,m+1}+{\psi }_{n,m-1}\right)\right]\hfill \\ & \times & \left(1+\mu {\left|{\psi }_{n,m}\right|}^2\right)-2{\left|{\psi }_{n,m}\right|}^2{\psi }_{n,m},\hfill \end{array}$$

(48)

Similar to its 1D version (11), Eq. (48) conserves the norm and Hamiltonian, cf. Eqs. (12) and (14),

$${\left(N_{\mathrm{SM}}\right)}_{2\mathrm{D}}=(1/\mu )\sum _{m,n}ln\left|1+\mu |{\psi }_{n,m}{|}^2\right|,$$

(49)

$$\begin{array}{c}{\left(H_{\mathrm{SM}}\right)}_{2\mathrm{D}}=-\sum _{n,m}\left[\left({\psi }_{n,m}{\psi }_{n+1,m}^{*}+{\psi }_{n+1,m}{\psi }_{n,m}^{*}\right)+\left({\psi }_{n,m}{\psi }_{n,m+1}^{*}+{\psi }_{n,m+1}{\psi }_{n,m}^{*}\right)\right.\\ \left.+(2/\mu )|{\psi }_{n,m}{|}^2\right]+(2/\mu ){\left(N_{\mathrm{SM}}\right)}_{2\mathrm{D}}.\end{array}$$

(50)

The continuum limit of this model is the 2D equation which is an extension of its 1D counterpart (17):

$$i{\Psi }_t+\left(1+\mu {\left|\Psi \right|}^2\right)\left({\Psi }_{xx}+{\Psi }_{yy}\right)+2\left(2\mu +1\right){|\Psi |}^2\Psi =0.$$

(51)

Note that term in $\mu {\left|\Psi \right|}^2\left({\Psi }_{xx}+{\Psi }_{yy}\right)$ prevents the onset of the collapse in Eq. (51), making it possible to create stable 2D solitons.

2D discrete solitons are looked for as solutions to Eq. (51) in the usual form, ${\psi }_{m,n}\left(t\right)=e^{-i\omega t}{\Phi }_{m,n}$. In the most interesting case of the competing nonlinearities, $\mu <0$, the situation is similar to that outlined above for the one-dimensional SM: there are ordinary discrete solitons, which have their stability and instability regions, and 2D cuspons, which are entirely stable in their existence region. Typical 2D solitons of both types are displayed in Figure 8. Antisymmetric bound states of ordinary 2D discrete solitons, and symmetric complexes built of 2D cuspons, are stable, while the bound states of cuspons with opposite parities are unstable, also like in the 1D model.

Along with the fundamental solitons, the 2D SM with the competing nonlinearities gives rise to vortex-soliton modes which may be stable in narrow parameter regions [111]. Examples of onsite- and intersite-centered vortex solitons (alias vortex cross and vortex square, respectively) are presented in Figure 9. In the 2D SM with non-competing nonlinearities ($\mu >0$ in Eq. (51)), vortex solitons are unstable, spontaneously transforming into fundamental ones and losing their vorticity. This transition is possible because the angular momentum is not conserved in the lattice system. The situation is different in the 2D SM with competing nonlinearities ($\mu <0$), where unstable vortex modes transform into vortical breathers, i.e., persistently oscillating localized modes that keep the original vorticity.

4. Solitons of the Semi-Vortex (SV) and Mixed-Mode (MM) Types in the Discrete 2D Spin-Orbit-Coupling (SOC) System

Recently, much interest has been drawn to emulation of the spin-orbit-coupling (SOC) effect, which is well known in physics of semiconductors, in binary BEC [112,113,114,115]. While SOC is a linear effect, its interplay with the intrinsic mean-field nonlinearity of atomic BEC gives rise to predictions of new species of 1D, 2D, and 3D solitons [116]. In particular, the effectively 2D binary BEC with SOC of the Rashba type is modeled by the following system of coupled GP equations for two components ${\varphi }^{(\pm )}$ of the pseudo-spinor wave function [117],

$$\begin{array}{ccc}\hfill i\frac{\partial {\varphi }^{(+)}}{\partial t}& =& -\frac{1}{2}{\nabla }^2{\varphi }^{(+)}-\left(\right|{\varphi }^{(+)}{|}^2+\gamma |{\varphi }^{(-)}{|}^2){\varphi }^{(+)}+2\lambda \left(\frac{\partial {\varphi }^{(-)}}{\partial x}-i\frac{\partial {\varphi }^{(-)}}{\partial y}\right),\hfill \\ \hfill i\frac{\partial {\varphi }^{(-)}}{\partial t}& =& -\frac{1}{2}{\nabla }_{-}^2{\varphi }^{(-)}-\left(\right|{\varphi }^{(-)}{|}^2+\gamma |{\varphi }^{(+)}{|}^2){\varphi }^{(-)}+2\lambda \left(-\frac{\partial {\varphi }^{(+)}}{\partial x}-i\frac{\partial {\varphi }^{(+)}}{\partial y}\right).\hfill \end{array}$$

(52)

In this system, SOC of the Rashba type is represented in by terms with coefficient $\lambda$, which couple the two equations through the first-order spatial derivatives. The system also includes the self- and cross-attractive nonlinearities, with scaled coefficients 1 and $\gamma$, respectively.

The system of coupled GP equations (52) maintains 2D solitons of two different types, namely, semi-vortices (SVs) and mixed modes (MMs) [117], The SV solitons, written in polar coordinates $\left(r,\theta \right)$, have zero vorticity in one component, and vorticity $+1$ or $-1$ in the other:

$$\begin{array}{ccc}\hfill {\varphi }_1^{(+)}& =& e^{-i\omega t}{f}_1\left(r^2\right),{\varphi }_1^{(-)}=e^{-i\omega t+i\theta }r{f}_2\left(r^2\right),\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill {\varphi }_2^{(+)}& =& e^{-i\omega t-i\theta }r{f}_2\left(r^2\right),{\varphi }_2^{(-)}=e^{-i\omega t}{f}_1\left(r^2\right),\hfill \end{array}$$

(54)

where $\omega$ is the chemical potential, and $f_{1,2}\left(r^2\right)$ are real functions which take finite values at $r=0$ and exponentially decay $\sim \left(sin\left(2\lambda r\right),cos\left(2\lambda r\right)\right)exp\left(-\sqrt{2\left(-\omega -2{\lambda }^2\right)}r\right)$ at $r\to \infty$. The two SV solutions (53) and (), which are mirror images of each other, exist in the semi-infinite bandgap, $\omega <-2{\lambda }^2$.

The combination of zero and nonzero vorticities in the SV solutions (53) and () is exactly compatible with the structure of the coupled GP equations (52). On the contrary to this, MM solitons cannot be represented by an exact ansatz similar to Eqs. (53) and (), but they may be approximated by a linear combination of both types of the SVs, $\left\{{\varphi }_1^{(+)}+{\varphi }_2^{(+)},{\varphi }_1^{(-)\left(1\right)}+{\varphi }_2^{(+)}\right\}$. An essential result is that the SVs and MMs are stable and represent the system’s ground state in the cases of $\gamma <1$ and $\gamma >1$, respectively, i.e., when the self-attraction is stronger or weaker than the cross-attraction in Eqs. (52) [117]. On the other hand, the SVs and MMs are unstable, as excited states, in the opposite cases, i.e., $\gamma >1$ and $\gamma <1$, respectively.

The discretized version of the SOC GP system (52), which corresponds to the spin-orbit-coupled binary BEC trapped in a deep OL potential, with discrete coordinates $(m,n)$, was introduced in Ref. [118]:

$$\begin{array}{c}i\frac{\partial {\varphi }_{m,n}^{(+)}}{\partial t}=-\frac{1}{2}\left({\varphi }_{m+1,n}^{(+)}+{\varphi }_{m-1,n}^{(+)}+{\varphi }_{m,n+1}^{(+)}+{\varphi }_{m,n-1}^{(+)}-4{\varphi }_{m,n}^{(+)}\right)\\ -\left(\right|{\varphi }_{m,n}^{(+)}{|}^2+\gamma |{\varphi }_{-}{|}^2){\varphi }_{m,n}^{(+)}+\lambda \left[{\varphi }_{m+1,n}^{(-)}-{\varphi }_{m-1,n}^{(-)}-i\left({\varphi }_{m,n+1}^{(-)}-{\varphi }_{m,n-1}^{(-)}\right)\right],\\ i\frac{\partial {\varphi }_{m,n}^{(-)}}{\partial t}=-\frac{1}{2}\left({\varphi }_{m+1,n}^{(-)}+{\varphi }_{m-1,n}^{(-)}+{\varphi }_{m,n+1}^{(-)}+{\varphi }_{m,n-1}^{(-)}-4{\varphi }_{m,n}^{(-)}\right)\\ -\left(\right|{\varphi }_{m,n}^{(-)}{|}^2+\gamma |{\varphi }_{m,n}^{(+)}{|}^2){\varphi }_{-}+\lambda \left[-\left({\varphi }_{m+1,n}^{(+)}-{\varphi }_{m-1,n}^{(+)}\right)-i\left({\varphi }_{m,n+1}^{(+)}-{\varphi }_{m,n-1}^{(+)}\right)\right].\end{array}$$

(55)

The linearized version of this system gives rise to the following dispersion relation for the plane waves, ${\varphi }_{m,n}^{(\pm )}\sim exp\left(-i\omega t+ipx+iqy\right)$, with wavenumbers taking values in the first Brillouin zone, $0<p,q<2\pi$:

$$\omega =2\left({sin}^2\frac{p}{2}+{sin}^2\frac{q}{2}\right)\pm 2\lambda \sqrt{{sin}^{2}p+{sin}^{2}q}.$$

(56)

The numerical solution of Eq. (55) has produced 2D modes which are discrete counterparts of the SV and MM solitons of the continuum system (52), see examples in Figure 10. As concerns the stability, the discreteness extends the stability of the SV and MM solitons towards $\gamma >1$ and $\gamma <1$, respectively.

A drastic difference of the discrete solitons of both the SV and MM types from their counterparts in the continuum is that they suddenly suffer delocalization (decay) when the SOC strength $\lambda$ in Eq. (55) exceeds a certain critical value, ${\lambda }_{\mathrm{cr}}$. The dependence of ${\lambda }_{\mathrm{cr}}$ on the soliton’s norm,

$$N=\sum _{m,n}\left({\left|{\varphi }_{m,n}^{(+)}\right|}^2+{\left|{\varphi }_{m,n}^{(-)}\right|}^2\right),$$

(57)

for the SV and MM solitons is displayed in Figure 11. The onset of the delocalization may be explained as a transition of the solution from the spectral bandgap to the band populated by the small-amplitude plane-wave states in the system’s linear spectrum, which is produced by Eq. (56).

5. Stable Soliton Species in the 3D DNLS Equation

5.1. The 3D Setting

A natural development of the analysis of the solitons and solitary vortices, and their bound states, produced by the 2D discrete DNLS equation and its extensions, which is outlined above in Section 2 – 4, is to construct self-trapped states (solitons) in the framework of the 3D equation:

$$\begin{array}{c}i{\dot{\varphi }}_{l,m,n}+C\left({\varphi }_{l+1,m,n}+{\varphi }_{l,m+1,n}+{\varphi }_{l,m,n+1}+{\varphi }_{l-1,m,n}+{\varphi }_{l,m-1,n-1}+{\varphi }_{l,m,n-1}-6{\varphi }_{l,m,n}\right)\\ +{\left|{\varphi }_{l,m,n}\right|}^2{\varphi }_{l,m,n}=0,\end{array}$$

(58)

where, as above, the overdot stands for the time derivative, $\left(l,m,n\right)$ is the set of the 3D discrete coordinates, and $C>0$ is the coefficient of the intersite coupling. The 3D DNLS equation cannot be realized in optics, but it admits natural implementation for BEC loaded in a deep 3D OL potential [97,98]. In that case, ${\varphi }_{l,m,n}\left(t\right)$ is the the respective effectively discretized BEC wave function.

As above, stationary soliton solutions to Eq. (58) with chemical potential $\omega$ are looked for as

$${\varphi }_{l,m,n}=exp(-i\omega t)u_{l,m,n},$$

(59)

where the stationary discrete wave function $u_{l,m,n}$ obeys the corresponding equation,

$$\begin{array}{c}\omega u_{l,m,n}+C\left(u_{l+1,m,n}+u_{l,m+1,n}+u_{l,m,n+1}+u_{l-1,m,n}+u_{l,m-1,n-1}+u_{l,m,n-1}-6u_{l,m,n}\right)\\ +{\left|u_{l,m,n}\right|}^2{u}_{l,m,n}=0.\end{array}$$

(60)

In particular, numerical solutions of Eq. (60) for 3D discrete solitons with embedded vorticity $S=0,1,2,...$ ($S=0$ corresponds to the fundamental solitons, for which the wave function $u_{l,m,n}$ is real) can be obtained, starting from the natural input

$$u_{l,m,n}^{\left(\mathrm{init}\right)}=A{\left(l+im\right)}^S\mathrm{sech}\left(\eta \sqrt{l^2+m^2}\right)exp\left(-\left|n\right|\right),$$

(61)

where $\eta$ is a real scale parameter, and it is implied that the vorticity axis is directed along coordinate n [119].

It is also relevant to consider a two-component system of nonlinearly-coupled 3D DNLS equations, for wave functions ${\varphi }_{l,m,n}\left(t\right)$ and ${\psi }_{l.m.n}\left(t\right)$ of two interacting BEC species (most typically, these are different hyperfine states of the same atom) [119]:

$$\begin{array}{c}i{\dot{\varphi }}_{l,m,n}+C\left({\varphi }_{l+1,m,n}+{\varphi }_{l,m+1,n}+{\varphi }_{l,m,n+1}+{\varphi }_{l-1,m,n}+{\varphi }_{l,m-1,n-1}+{\varphi }_{l,m,n-1}-6{\varphi }_{l,m,n}\right)\\ +\left({\left|{\varphi }_{l,m,n}\right|}^2+\beta {\left|{\psi }_{l,m,n}\right|}^2\right){\varphi }_{l,m,n}=0,\\ i{\dot{\psi }}_{l,m,n}+C\left({\psi }_{l+1,m,n}+{\psi }_{l,m+1,n}+{\psi }_{l,m,n+1}+{\psi }_{l-1,m,n}+{\psi }_{l,m-1,n-1}+{\psi }_{l,m,n-1}-6{\psi }_{l,m,n}\right)\\ +\left({\left|{\psi }_{l,m,n}\right|}^2+\beta {\left|{\varphi }_{l,m,n}\right|}^2\right){\psi }_{l,m,n}=0,\end{array}$$

(62)

Here $\beta >0$ is the relative strength of the inter-component attractive interaction with respect to the intra-component self-attraction.

5.2. Results

5.2.1. Single-Component 3D Solitons

The numerical analysis, starting from input (61), has provided families of fundamental and vortex solitons. Here, following Ref. [119], the results are displayed for a fixed value of the chemical potential, $\omega =-2$ in Eq. (59), while varying coupling constant C in Eqs. (58) and (60). In particular, the discrete fundamental solitons with $S=0$ are stable at $C<$$C_{\mathrm{cr}}^{\left(0\right)}\approx 2$, and the vortex modes with $S=1$ are stable at $C<$$C_{\mathrm{cr}}^{\left(1\right)}\approx 0.65$. Note that the limit of $C\to \infty$ corresponds to the 3D NLS equation in the continuum, in which all solitons are definitely unstable; therefore, all discrete solitons become unstable at sufficiently large values of C. At $C>C_{\mathrm{cr}}^{\left(1\right)}$ the simulations demonstrate that the development of the instability destroys the vortical structure and, eventually, transforms the soliton into a fundamental one, with $S=0$ (not shown here in detail). The change of the topological charge is possible, as the angular momentum is not a dynamical invariant of the lattice dynamics.

The vortex solitons with $S=2$ are completely unstable, but an unusual feature of these states is that, at sufficiently small values of C (in particular, at $C=0.01$), the instability spontaneously rebuilds them into stable discrete solitons with a larger vorticity, $S=3$ [119]. An example of a stable soliton with $S=3$ is displayed in Figure 12.

In addition to the fundamental and vortex solitons, Eqs. (58) and (60) produce diverse multimode species of stable discrete 3D solitons in the form of dipoles, quadrupoles and octupoles [120]. Examples of such states are presented in Figure 13 for $C=0.1$. This figure displays tightly-bound dipoles with different orientations with respect to the lattice, $viz$., straight, 2D-diagonal, and 3D-diagonal ones (panels (a,b,c)), quadrupole (panel (d)), and octupole (panel (f)), in which the field fills adjacent sites of the lattice (with the lattice distance between them $d=1$). Also displayed are loosely-bound quadrupole and octupole (in panels (e) and (g), respectively), with distance $d=2$ between the filled sites. Similar multimode states with still larger separations d between the filled sites were found too. The results are summarized in Figure 13(h), which shows stability boundaries $C_{cr}^{(3D,d)}$ for dipoles, quadrupoles, and octupoles vs. d. Naturally, the stability region, $C<C_{cr}^{(3D,d)}$, increases with the increase of d, as the interaction between the filled sites, which leads to the possible dynamical instability of the multipole states, is weaker for larger d.

In addition to the above-mentioned states, Eq. (60) admits more sophisticated stable composite states, such as “vortex cubes", built as a pair of identical parallel quasi-planar vortices with topological numbers $S_1=S_2=1$, with opposite signs (phase shift $\pi$). set in parallel planes, as shown in Figure 14(a). Stationary solutions representing vortex-antivortex cubes, in the form of parallel quasi-planar vortices with opposite topological charges, $S_1=-S_2=1$, can be found too, as shown in Figure 14(b), but they are completely unstable.

The same Eq. (58) gives rise to other stable self-trapped modes, such as vortex solitons with the axis directed along the 2D diagonal, cf. Figure 13(b). Vortex modes with the axis parallel to the 3D diagonal exist too, but they are unstable, see further details in Ref. [120].

Two-Component 3D Solitons (Including Skyrmions)

The system of coupled 3D DNLS equations (62) produces stable soliton states which are specific to the two-component nonlinear lattice medium. A noteworthy example is a composite mode built as a bound state of vortex solitons in the two components with mutually orthogonal orientations, see an example in Figure 15. These bound states are stable for sufficiently small values of the coupling constant, such as $C=0.01$ in Figure 15, and for $\beta <1$ in Eq. (62), i.e., under the condition that the self-attraction in each component is stronger than the inter-component attraction.

The system of coupled GP equations with the repulsive sign of the nonlinearity may be used as the simplest model producing skyrmions in the binary BEC [121,122,123]. The discretization of the GP system leads to Eqs. (62) with the opposite sign in front of the nonlinear terms [124]. Then, these equations are reduced to stationary ones by the usual substitution with chemical potential $\omega$, $\left\{\varphi ,\psi \right\}=exp\left(-i\omega t\right)\left\{u_{l,m,n},v_{l,m,n}\right\}$:

$$\begin{array}{c}C\left(u_{l+1,m,n}+u_{l,m+1,n}+u_{l,m,n+1}+u_{l-1,m,n}+u_{l,m-1,n-1}+u_{l,m,n-1}-6u_{l,m,n}\right)\\ -\left({\left|u_{l,m,n}\right|}^2+\beta {\left|v_{l,m,n}\right|}^2\right)u_{l,m,n}=-\omega u_{l,m,n},\\ C\left(v_{l+1,m,n}+v_{l,m+1,n}+v_{l,m,n+1}+v_{l-1,m,n}+v_{l,m-1,n-1}+v_{l,m,n-1}-6v_{l,m,n}\right)\\ -\left({\left|v_{l,m,n}\right|}^2+\beta {\left|u_{l,m,n}\right|}^2\right)v_{l,m,n}=-\omega v_{l,m,n},\end{array}$$

(63)

where the relative strength $\beta$ of the inter-component repulsion with respect to the self-repulsion remain a positive coefficient. For $\omega >0$, skyrmions can be constructed by choosing field $u_{l,m,n}$ as a complex one, representing a quasi-flat vortex soliton with topological charge $S=1$, and real field $v_{l,m,n}$ as a bubble into which the vortex soliton is embedded, with a nonzero background value at $\left(\left|l\right|,\left|m\right|,\left|n\right|\right)\to \infty$, viz., $v_{\mathrm{background}}^2=\omega$ [124]. An example of a numerically found skyrmion solution of this type is displayed in Figure 16.

The same work [124] reported solutions for 2D discrete “baby skyrmions", which are lattice counterparts of the modes produced by the 2D reduction of the Skyrme model [125,126]. The have a simple structure similar to its 3D counterpart displayed in Figure 16, i.e., a complex 2D vortex soliton in one component, embedded into a bubble of the delocalized field in the other real component.

6. 2D solitons and Solitary Vortices in Semi-Discrete Systems

6.1. Spatiotemporal Optical Solitons in Arrayed Waveguides

The consideration of 2D and 3D settings suggests a natural option to introduce 2D semi-discrete systems, with a continuous coordinate in one direction and a discrete coordinate in another, as well as 3D systems, where one or two coordinates are continuous, while the remaining two or one coordinates are discrete. In optics, a well-known 2D setting belonging to this class represents spatiotemporal propagation of light in an array of optical fibers [127]. It is modeled by the system of linearly-coupled NLS equations for complex amplitudes $u_n\left(z,\tau \right)$ of optical fields in individual fibers:

$$i{\partial }_z{u}_n+(1/2)D{\partial }_{\tau }^2{u}_n+(\kappa /2)\left(u_{n+1}+u_{n-1}-2u_n\right)+{\left|u_n\right|}^2{u}_n=0,$$

(64)

where z is the propagation distance, $\tau \equiv t-x/V_{\mathrm{gr}}$ (with time t and carrier group velocity $V_{\mathrm{gr}}$) is the usual temporal variable, real D is the group-velocity-dispersion coefficient in each fiber, $\kappa >0$ is the coefficient of coupling between adjacent fibers in the array, and the nonlinearity coefficient is normalized to be 1. It is commonly known that optical solitons (semi-discrete ones, in the present case) can be supported in the case of anomalous dispersion, i.e., $D>0$.

A remarkable counter-intuitive property of semi-discrete localized modes generated by Eq. (64) is their ability to stably move across the array, under the action of a kick applied at $z=0$ [71]:

$$u_n\left(\tau \right)\to exp\left(ian\right)u_n\left(\tau \right),$$

(65)

with real a, in spite of the presence of the respective quasi-1D Peierls-Nabarro potential, An example of such a moving mode is displayed in Figure 17. This property may be compared to the above-mentioned mobility of 1D discrete solitons in the DNLS equation [65], and of 2D discrete solitons in the framework of the ${\chi }^{\left(2\right)}$ system (25).

Similarly, quasi-discrete settings modeled by an extension of (64) with two transverse spatial coordinates were used for the creation for spatiotemporal optical solitons (“light bullets") [128], as well as soliton-like transient modes with embedded vorticity [129]. Waveguides employed in those experiments feature a transverse hexagonal-lattice structure, written in bulk silica by means of an optical technology. A spatiotemporal vortex state observed in the bundle-like structure (in the experiment, it is actually a transient one) is represented by Figure 18, which displays both numerically predicted and experimentally observed distributions of intensity of light in the transverse plane, together with a phase plate used in the experiment to embed the vorticity into the incident spatiotemporal pulse which was used to create the mode.

6.2. Semi-Discrete Quantum and Photonic Droplets

A new type of semi-discrete solitons was recently elaborated in Ref. [73], in the framework of an array of linearly coupled 1D GP equations, including the above-mentioned Lee-Huang-Yang correction, which represents an effect of quantum fluctuations around the mean-field states of a binary BEC [29,33]. The system is

$$i{\partial }_t{\psi }_j=-(1/2){\partial }_{zz}{\psi }_j-\left(C/2\right)\left({\psi }_{j+1}-2{\psi }_j+{\psi }_{j-1}\right)+g|{\psi }_j{|}^2{\psi }_j-\left|{\psi }_j\right|{\psi }_j,$$

(66)

where ${\psi }_j\left(z\right)$ is the mean-field wave function in the j-th core with coordinate z, C is the effective inter-core coupling constant, the self-attractive quadratic term represents the Lee-Hung-Yang correction in the 1D limit, cf. Eq. (7), and $g>0$ accounts for the mean-field self-repulsion.

A semi-discrete system similar to the one modeled by Eq. (66), but with the cubic-quintic nonlinearity instead of the combination of the quadratic and cubic terms in Eq. (66), was derived in the context of nonlinear optics [74]:

$$i{\partial }_z{u}_n=-(1/2){\partial }_{xx}{u}_n-(C/2)\left(u_{n+1}-2u_n+u_{n-1}\right)-|u_n{|}^2{u}_n+{\left|u_n\right|}^4{u}_n.$$

(67)

It corresponds to the array of parallel-coupled planar waveguides, as shown in Figure 19. In this case, $u_n(x,z)$ is the complex local amplitude of the optical wave in the n-th waveguide, z is the propagation distance, and x is the transverse coordinate in each waveguide, while C is the effective coupling constant, similar to Eq. (66). By analogy with the quantum droplets, semi-discrete solitons produced by Eq. (67) may be called “photonic droplets".

The droplets produced by Eqs. (66) and (67) are characterized by the total norm, which is proportional to the number of atoms in BEC,

$$N=\sum _j{\int }_{-\infty }^{+\infty }{\left|{\psi }_j\left(z\right)\right|}^{2}dz,$$

(68)

or the total power of the photonic droplet,

$$P=\sum _n{\int }_{-\infty }^{+\infty }{\left|u_n\left(x\right)\right|}^{2}dx.$$

(69)

For solitons produced by Eqs. (66) and (67) sets of control parameters are, respectively, $\left(C,g\right)$ for fixed N, or $\left(C,P\right)$.

The models based on Eqs. (66) and (67) give rise to many families of semi-discrete solitons, including a novel species of semi-discrete vortex solitons. Typical examples of the onsite- and intersite vortices with topological charge $S=1$, produced by Eq. (66), are displayed in Figure 20.

An example of a stable semi-discrete vortex soliton produced by Eq. (67) with $S=2$ is displayed in Figure 21.

Getting back to the semi-discrete system (66), a chart in the plane of $\left(C,g\right)$ which displays stability areas for the semi-discrete vortex solitons with topological charges $S=2$, 3, 4, and 5, is plotted in Figure 22. The chart demonstrates abundant multistability – for instance, the stable solitons with $S=2$ coexist with all higher-order ones (with $S=3$, 4, and 5).

Self-trapped solutions of a continuum model, which are similar to semi-discrete vortex solitons outlined above, were recently reported for a photonic crystal built in a ${\chi }^{\left(2\right)}$ material with a checkerboard structure representing quasi-phase matching [75].

Semi-discreteness of another type is possible in two-component systems, where one component is governed by a discrete equation, and the other one by a continuous equation. This type of two-component systems was proposed in [76]. It introduced a ${\chi }^{\left(2\right)}$ model, assuming that the continuous second-harmonic wave propagates in a slab with a continuous transverse coordinate, while the fundamental-harmonic field is trapped in a discrete waveguiding array built on top of the slab.

7. 2D Fundamental and Vortical Discrete Solitons in a Two-Component $\mathcal{PT}$ (parity-time) symmetric lattice

While the above presentation deals solely with conservative discrete systems, many properties of conservative settings are shared by a very special type of dissipative ones, viz., systems with the parity-time ($\mathcal{PT}$) symmetry. They include mutually symmetric spatially separated elements carrying linear gain and loss [138,139,140]. The experimental realization of such systems in optics [140] suggests one to include the Kerr nonlinearity, thus opening the way to the prediction an creation of $\mathcal{PT}$-symmetric solitons [141,142]. In particular, exact solutions for 1D $\mathcal{PT}$-symmetric solitons and exact results for their stability boundaries were found in the model of the nonlinear $\mathcal{PT}$-symmetric coupler (dual-core waveguide), with mutually symmetric linear gain and loss carried by the linearly coupled cores [143,144]. Stability limits for 2D fundamental solitons in the 2D $\mathcal{PT}$-symmetric coupler with the cubic-quintic nonlinearity in each core (essentially the same as in Eqs. (67), chosen to prevent the critical-collapse instability) were identified in Ref. [145].

The definition of the $\mathcal{PT}$ symmetry makes it also natural to consider discrete $\mathcal{PT}$-symmetric systems. Various species of stable discrete solitons were predicted in chains of $\mathcal{PT}$-symmetric elements [83,146,147,148,149,150], and the existence of such solitons was demonstrated experimentally [151].

A natural model for the creation of $\mathcal{PT}$-symmetric discrete 2D solitons is a generalization of the 2D discrete nonlinear coupler, based on Eqs. (44), by adding the linear gain and loss terms with strength $\gamma >0$ to the coupled equations [83]:

$$\begin{array}{ccc}& & i\frac{d{\psi }_{m,n}}{dz}=-\frac{C}{2}({\psi }_{m,n+1}+{\psi }_{m,n-1}+{\psi }_{m-1,n}+{\psi }_{m+1,n}-4{\psi }_{m,n})\hfill \\ & & -|{\psi }_{m,n}{|}^2{\psi }_{m,n}-{\phi }_{m,n}+i\gamma {\psi }_{m,n},\hfill \\ & & i\frac{d{\phi }_{m,n}}{dz}=-\frac{C}{2}({\phi }_{m,n+1}+{\phi }_{m,n-1}+{\phi }_{m-1,n}+{\phi }_{m+1,n}-4{\phi }_{m,n})\hfill \\ & & -|{\phi }_{m,n}{|}^2{\phi }_{m,n}-{\psi }_{m,n}-i\gamma {\phi }_{m,n}.\hfill \end{array}$$

(70)

Here, in terms of the optical realization, the evolution variable z is the propagation distance, the inter-core coupling coefficient is scaled to be 1, and $C>0$ is constant of the intra-core coupling between adjacent sites of the lattice. The dispersion relations for plane-wave solutions to the linearized version of Eqs. (70), $\left\{{\psi }_{m,n}\left(z\right),{\phi }_{m,n}\left(z\right)\right\}\sim exp\left(ipm+iqn+ikz\right)$, is

$$k=-2C\left({sin}^2\frac{p}{2}+{sin}^2\frac{q}{2}\right)\pm \sqrt{1-{\gamma }^2}.$$

(71)

As it follows from Eq. (71), the $\mathcal{PT}$ symmetry holds under condition $\gamma <{\gamma }_{max}\equiv 1$, i.e., the gain-loss strength $\gamma$ must be smaller than the linear-coupling coefficient, that is 1 in the present notation, which is a generic property of $\mathcal{PT}$-symmetric couplers [143,144].

Stationary modes with real propagation constant k are looked for as solutions to the full nonlinear system of Eqs. (70) in the usual form, $\left\{{\psi }_{m,n}\left(z\right),{\phi }_{m,n}\left(z\right)\right\}=e^{ikz}\left\{u_{m,n},v_{m,n}\right\}$, with stationary amplitudes obeying equations

$$\begin{array}{ccc}\hfill \left(k+i\gamma \right)u_{m,n}& =& \frac{C}{2}(u_{m,n+1}+u_{m,n-1}+u_{m-1,n}+u_{m+1,n}-4u_{m,n})+{\left|u_{m,n}\right|}^2{u}_{m,n}+v_{m,n},\hfill \\ \hfill \left(k-i\gamma \right)v_{m,n}& =& \frac{C}{2}(v_{m,n+1}+v_{m,n-1}+v_{m-1,n}+v_{m+1,n}-4v_{m,n})+{\left|v_{m,n}\right|}^2{v}_{m,n}+u_{m,n}.\hfill \end{array}$$

(72)

Localized states produced by Eqs. (72) are characterized, as above, by the total power,

$$P={\sum }_{m,n}\left(|u_{m,n}{|}^2+{\left|v_{m,n}\right|}^2\right).$$

(73)

Straightforward analysis of Eqs. (72) demonstrates that the system produces $\mathcal{PT}$-symmetric fundamental-soliton solutions, which must be subject to the relation $v_{m.n}=u_{m.n}^{*}$ (with * standing for the complex conjugate), in the form of

$$\left\{u_{m.n},v_{mn}\right\}=w_{m,n}exp\left(\pm (i/2)arcsin\gamma \right),$$

(74)

where real discrete distribution $w_{m,n}$ should be found as a solution of the usual stationary equation for 2D discrete solitons,

$$\left(-k+\sqrt{1-{\gamma }^2}\right)w_{m,n}+\frac{C}{2}\left(w_{m,n+1}+w_{m,n-1}+w_{m-1,n}+w_{m+1,n}-4w_{m,n}\right)+w_{m,n}^3=0,$$

(75)

cf. Eq. (38). In agreement with the linear spectrum (71), Eq. (75) may produce soliton solutions for $k>\sqrt{1-{\gamma }^2}$. An example of a stable fundamental $\mathcal{PT}$-symmetric soliton is displayed, by means of its cross-section shapes, in Figure 23.

The existence and stability of the $\mathcal{PT}$-symmetric fundamental discrete solitons is summarized in the plane of $\left(\gamma ,P\right)$ for the fundamental solitons in Figure 24. It is seen that, naturally, the stability area shrinks as the gain-loss coefficient $\gamma$ is approaching its limit value, ${\gamma }_{max}=1$ (cf. the 1D situation considered in Refs. [143] and [144]). The existence boundary, i.e., the minimum value of P, below which no solitons are found (in the white area), corresponds to the limit of very broad small-amplitude solitons. In this limit, the discrete soliton may be approximated by its counterpart in the continuum NLS equation, i.e., the above-mentioned Townes soliton, whose power takes the unique value, which thus determines the existence boundary in Figure 24.

The stability boundary in Figure 24 may be understood as the one at which the symmetric soliton is destabilized by the spontaneous symmetry breaking (as described in detail above for 2D solitons produced by the linearly-coupled conservative DNLS equations (44), see also Ref. [17]), which is here modified by the presence of the linear gain and loss. Because asymmetric solitons cannot exist in the system with the balanced gain and loss, the symmetry breaking always leads to either blowup or decay of the soliton [83]. In their stability region, the $\mathcal{PT}$-symmetric fundamental discrete solitons actually represent the system’s ground state [83].

Alongside the fundamental discrete $\mathcal{PT}$-symmetric solitons, the same system of Eqs. (72) produces $\mathcal{PT}$-symmetric vortex solitons, which also have their stability area, see details in Ref. [83]. An example of a stable $\mathcal{PT}$-symmetric vortex soliton is presented in Figure 25.

In addition to the $\mathcal{PT}$-symmetric solitons, Eqs. (72) give rise to anti-$\mathcal{PT}$-symmetric ones, defined by relation $v_{m.n}=-u_{m.n}^{*}$. They, as well as anti-$\mathcal{PT}$-symmetric vortex solitons, are stable in some parameter areas (see details in Ref. [83]), but those areas are essentially smaller than their counterparts for the $\mathcal{PT}$-symmetric modes. The reduced stability area for the anti-$\mathcal{PT}$-symmetric fundamental solitons is explained by the fact that they cannot be the system’s ground state.

8. Conclusion

The interplay of the discreteness and intrinsic nonlinearity in various physical media – chiefly, in nonlinear optics and BEC – gives rise to a great variety of self-trapped localized states, in the form of discrete solitons. This article aims to produce a concise review, starting from the brief survey of basic theoretical models combining the discreteness in 1D, 2D, and 3D geometries and various nonlinearities, such as cubic, quadratic, and quintic. The main subject addressed in the article is a summary of basic results for 2D and 3D discrete solitons produced by such models. Unlike the topic of 1D discrete solitons, the multidimensional ones were not previously reviewed in a systematic form. Along with the fundamental solitons, topologically organized ones, in the form of solitary vortices and discrete skyrmions, are considered too. Some experimental findings are also included, such as the observation of 2D discrete optical solitons with embedded vorticity.

In many cases, the discreteness helps to produce states which either do not exist or are definitely unstable in continuum analogs of the discrete settings. In particular, these are 2D fundamental and vortex solitons, which may be stable in the discrete form, while their continuum counterparts are completely unstable in the free space. On the other hand, mobility of solitons, which is their obvious property in the continuum, is a nontrivial issue for the lattice (discrete) solitons.

The work in this area remains a subject of ongoing theoretical and experimental work, promising new findings. A perspective direction is to produce 2D and 3D self-trapped states with intrinsic topological structures. Some results obtained in this direction are presented in this article, such as discrete solitons in the system with spin-orbit coupling [118] (see also Ref. [130]), sophisticated 3D discrete modes with embedded vorticity [119,120], and discrete skyrmions [124]. A challenging task is experimental realization of these states which, thus far, were only predicted in the theoretical form.

It is relevant to mention some topics which may be relevant in the present context but are not included here, to keep a reasonable size of the review. In particular, these are interactions of discrete solitons with local defects in the underlying lattice, as well as with interfaces and edges. It is known that defects and surfaces may often help to create and stabilize localized modes which do not exist or are unstable in uniform lattices, such as Tamm [131] and topological-insulator [132,133] states. Another vast area of studies, which is not considered here, deals with dissipative discrete nonlinear systems. In this article, only the very special case of $\mathcal{PT}$-symmetric systems is addressed. Basic nonlinear dissipative models are represented by discrete complex Ginzburg-Landau equations, i.e., DNLS equations with complex coefficients in front of the onsite linear and nonlinear terms, which account for losses and compensating gain [134]. Unlike conservative and $\mathcal{PT}$-symmetric models, the dissipative ones may only give rise to stable discrete solitons which do not exist in continuous families, but rather as isolated attractors [135,136,137].