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.

The basic 2D cubic DNLS equation is the 2D version of Eq. (3) with the self-attraction ($\sigma =-1$) and without the external potential: 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: 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].

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: 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.

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: where the horizontal and vertical coupling constants are modulated as follows: 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].

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: 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]: 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].

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] 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, and the asymmetry of the symmetry-broken state is determined by parameter

The system under the consideration can be analyzed by means of the VA, based on the 2D ansatz 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.

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: 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 where the stationary discrete wave function $u_{l,m,n}$ obeys the corresponding equation, 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 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]: Here $\beta >0$ is the relative strength of the inter-component attractive interaction with respect to the intra-component self-attraction.