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

This article presents a concise survey of basic discrete and semi-discrete nonlinear models, which produce two- and three-dimensional (2D and 3D) solitons, and a summary of the main theoretical and experimental results obtained for such solitons. The models are based on the discrete nonlinear Schrödinger (DNLS) equations and their generalizations, such as a system of discrete Gross–Pitaevskii (GP) equations with the Lee–Huang–Yang corrections, the 2D Salerno model (SM), DNLS equations with long-range dipole–dipole and quadrupole–quadrupole interactions, a system of coupled discrete equations for the second-harmonic generation with the quadratic (χ(2)) nonlinearity, a 2D DNLS equation with a superlattice modulation opening mini-gaps, a discretized NLS equation with rotation, a DNLS coupler and its PT-symmetric version, a system of DNLS equations for the spin–orbit-coupled (SOC) binary Bose–Einstein condensate, and others. The article presents a review of the basic species of multidimensional discrete modes, including fundamental (zero-vorticity) and vortex solitons, their bound states, gap solitons populating mini-gaps, symmetric and asymmetric solitons in the conservative and PT-symmetric couplers, cuspons in the 2D SM, discrete SOC solitons of the semi-vortex and mixed-mode types, 3D discrete skyrmions, and some others.


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 ψ (x, y, z; t) 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 where σ = +1 and −1 correspond to the self-defocusing and focusing signs of the cubic nonlinearity, and U (x, y, z) 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, τ = t − z/V gr , where V 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], as well as its 2D and 1D reductions.A deep lattice potential, which corresponds to large amplitude ε 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]: where the overdot stands for d/dt, and the set of integer indices, (l, m, n), replaces continuous coordinates (x, y, z) in Eq. (1).In Eq. ( 3), potential V l,m,n is a possible smooth addition to the deep lattice potential which imposes the discretization.
The rigorous mathematical derivation of the DNLS equation by the discretization of the underlying continuum NLS equation with the deep spatially periodic potential is based on the expansion of the continuous wave field over the set of Wannier functions [13].These are linear combinations of the quasiperiodic Bloch wave functions which feature shapes localized around potential minima [14], thus offering a natural basis for the transition to the discrete limit.
The full DNLS equation ( 3) is often reduced to its 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 [15,16].2D lattices with similar additional couplings were elaborated too [17].

Extended equations
The DNLS equation and its extensions, such as systems of coupled DNLS equations [18,19], constitute a class of models with a large number of physical realizations [10,20,[22][23][24].They have also drawn much interest as subjects of mathematical studies [25].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 σ = −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 [26], 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 [27,28].Thus, the discreteness provides a general method for the stabilization of 2D and 3D solitons.
a.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 celebrated Lee-Huang-Yang effect [30]).For the symmetric binary condensate, with identical wave functions ψ of its components, the accordingly amended scaled GP equations (in the absence of the trapping potential) was derived by Petrov [31] iψ t = −(1/2)∇ 2 ψ − |ψ| 2 ψ + |ψ| 3 ψ. ( Surprisingly quickly, the QD modes predicted by Eq. ( 5) have been created experimentally, in the quasi-2D [32,33] and full 3D [34] forms.The reduction of the spatial dimension to 2D and 1D replaces Eq. ( 5) by the following GP equations, respectively [35]: 2D: 1D: Note, in particular, that in the 1D equation (7) the quantum correction is represented by the attractive term −|ψ|ψ, on the contrary to its repulsive counterpart +|ψ| 3 ψ 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 [36].A semi-discrete version of Eq. ( 7) is considered below, see Eq. (66).b.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, σ = ±1 [37][38][39][40], although 2D and 3D extensions of the NLS equation are nonintegrable.The straightforward discretization destroys the integrability of the 1D NLS equation [41,42].Nevertheless, the NLS equation admits a specially designed 1D discretization, which leads to an integrable discrete model, viz., the Ablowitz-Ladik (AL) equation [43]: where positive and negative values of the nonlinearity coefficient, µ, 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 [44].Equation (8) gives rise to an exact solution for solitons in the case of µ > 0. Setting µ ≡ +1 by means of rescaling, the solution is where β and α are arbitrary real parameters that determine the soliton's amplitude, A ≡ sinh β, its velocity, V ≡ ξ = 2β −1 (sinh β) sin α, and overall frequency Ω ≡ φ = −2 [(cosh β) cos α + (α/β) (sinh β) sin α].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 [72].Another integrable discrete model which admits exact solutions for moving solitons is the Toda-lattice equation for real coordinates x n (t) of particles with unit mass and exponential potential of interaction between adjacent ones [73]: Considerable interest was also drawn to the nonintegrable combination of the AL and DNLS equations, in the form of the Salerno model (SM) [45], with an additional onsite cubic term, different from the intersite one in Eq. ( 8): 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 [46].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, [ 43,114].The Hamiltonians of the AL and SM equations, which are dynamical invariants too, are 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 ψ n as dψ n /dt = {H, ψ n }.For the AL and SM models, the Poisson brackets, written for a pair of arbitrary functions B (ψ n , ψ * n ) and C (ψ n , ψ * n ), are [45,46] {B, 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, where Ψ(x) 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 [57] Equation ( 17) conserves the total norm and Hamiltonian, which are continuum counterparts of expressions ( 12) and ( 14): It is relevant to mention that the general approximation opposite to the continuum limit is the anti-continuum limit [20,58].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 ψl,m,n = σ |ψ l,m,n | 2 ψ l,m,n +V l,m,n ψ l,m,n , at a finite set of sites, and keeping the zero solution at all others.The single-site "simple solutions" are 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).c.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 σ > 0 in Eq. ( 3), may support discrete-soliton (selftrapped) 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 → ∞ (r is the radial coordinate) faster than r D [47].In terms of the 1D and 2D DNLS equations, similar settings were introduced in Refs.[48] and [49], with the site-dependent self-attraction coefficients, (σ n ) 1D = σ 0 exp (α|n|) and (σ m,n ) 2D = σ 0 exp (α (|m| + |n|)), respectively, with positive constants σ 0 and α.In the 2D model, solutions were constructed, and their stability analyzed, for fundamental, dipole, quadrupole, and vortical discrete solitons [47].
d. 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 [50].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 Γ > 0 [51] where σ 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. [53]: The analysis reported in Ref. [51] 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 [51].
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. [54]: cf. Eq. ( 21).This model also gives rise to families of stable 2D discrete solitons [54].e.The 2D discrete second-harmonic-generating (χ (2) ) system The quadratic (alias χ (2) ) 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 χ (2) system for amplitudes ψ (x, y, z) and ϕ (x, y, z) of the fundamental-frequency (FF) and second-harmonic (SH) waves is [55] where z is the propagation distance, the paraxial-diffraction operator (1/2)∇ 2 acts on the transverse coordinates (x, y), 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 χ (2)  material, is [56] 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 χ (2) system is played by the Manley-Rowe invariant, i.e., the total optical power, An essential property of 2D discrete solitons produced by Eqs. ( 25) is their mobility [56].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 χ (2) nonlinearity does not give rise to the collapse in the 2D (and 3D) space, therefore 2D χ (2) solitons do not demonstrate a trend for strong pinning, remaining effectively mobile robust localized modes [56].

Fundamental solitons
In the 1D setting, the model of basic interest is the DNLS equation with self-attraction, which corresponds to σ = −1 in the 1D version of Eq. ( 3), without the external potential (V l,m,n = 0): This equation conserves two dynamical invariants, viz., the total norm, and Hamiltonian (energy), Stationary solutions to Eq. ( 26) with real frequency ω are looked for as with real amplitudes u n satisfying the discrete equation, Note that Eq. ( 30) can be derived by varying its Lagrangian, 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 [59], ψ n = a exp ia 2 t , 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) [25].The VA is based on a particular ansatz, i.e. an analytical expression which aims to approximate the solution [60].The only discrete ansatz for which analytical calculations are feasible is represented by the exponential function [61][62][63][64], namely, with a > 0. The corresponding norm, calculated as per Eq. ( 27), is Note that ansatz (32) is appropriate for strongly and moderately discrete solitons, as shown in Fig. 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), with width η −1 which must be large in comparison with the discreteness spacing, η −1 ≫ 1, and central coordinate ξ.
The substitution of ansatz (32) in Lagrangian (31) produces the corresponding VA Lagrangian: Then, for given ω < 0 (solitons do not exist for ω > 0), the squared amplitude, A 2 , and inverse width, a, of the discrete soliton are predicted by the Euler-Lagrange equations, 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 Fig. 1 and Ref. [65].Rigorous justification of the VA was elaborated in Ref. [71].Furthermore, the VA and a full numerical solution of Eq. ( 26) demonstrate that the entire family of the discrete solitons is stable [25].
In addition to the bright solitons considered here, the DNLS equation also gives rise to discrete dark solitons, which have been studied in detail theoretically and experimentally [66][67][68][69].As concerns the topic of the present review, two-dimensional discrete dark modes, such as delocalized lattice vortices, were studied too [70].However, the consideration of dark modes is not included in this article.

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 [74].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 [75,76], 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 [77].Note that stationary bound states of fundamental solitons do not exist in the continuum limit.

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 : cf. (30).Discrete solitons produced by the SM equation ( 11) with µ > 0, i.e. with noncompeting intersite and onsite self-focusing nonlinearities, were investigated by means of numerical methods [114][115][116].It was found that the SM gives rise to the 1D solitons at all positive values of µ.
Another option is to consider the SM with µ < 0, which features competing nonlinearities, as the intersite cubic term, with coefficient µ < 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 ω < 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 [75,76].The same is true for bound states of ordinary discrete solitons in the SM [57].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 [57].

C. The subject and structure of the present article
The above-mentioned reviews [10,20,[22][23][24] and [25] 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 [27,28].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 [78][79][80][81]84].The article chiefly represents theoretical results, but some experimental findings for quasi-discrete 2D solitons in photonic lattices [86,87] are included too.
The review presented below does not claim to be comprehensive.It 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 [88], and Frenkel-Kontorova systems [89] -are not considered here.Dissipative systems are not considered either, except for a 2D model with the parity-time (PT ) symmetry [90], 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 spinorbit-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 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.

II. TWO-DIMENSIONAL (2D) NONLINEAR-SCHR ÖDINGER LATTICES: FUNDAMENTAL AND VORTEX SOLITONS, AND THEIR BOUND STATES
A. Vortex solitons: theoretical and experimental results The basic 2D cubic DNLS equation is the 2D version of Eq. ( 3) with the self-attraction (σ = −1) and without the external potential: cf. Eq. ( 26).The substitution of ψ m,n = exp (−iωt) 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 [91,92], 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. [93] (see also Ref. [94]).Vorticity, alias the topological charge, or winding number, is defined as ∆φ/ (2π), where ∆φ 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 [95], which are unstable against the onset of the critical collapse [3].In the same limit, the Townes solitons with embedded vorticity (vortex rings [96]) are subject to much stronger instability against spontaneous splitting of the ring in fragments [97].
The lattice structure of the DNLS equation provides for stabilization of both fundamental (zero-vorticity) and vortex solitons [93].A typical example of a stable 2D vortex soliton with topological charge S = 1 is displayed in Fig. 2. 2D fundamental and vortex solitons, with topological charges S = 0 and 1, are stable in regions −ω > |ω 23, respectively [93], 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 [98].The vortex solitons with S = 3 may be stable, but only in an extremely discrete form, viz., at −ω > |ω 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 ω → 0.
The experimentally relevant lattice structure may be anisotropic, with the linear combination ( ) , with anisotropy parameter α ̸ = 1.Effects of the anisotropy on the structure and stability of the fundamental and vortical discrete solitons were explored in Ref. [109].
The theoretically predicted 2D discrete solitons with vorticity S = 1 were experimentally created in Refs.[86] and [87], using a photorefractive crystal.Unlike uniform media of this type, where delocalized ("dark") optical vortices were originally demonstrated [99,100], 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 Fig. 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 [101].Stable vortex solitons with S = 2 were created using a hexagonal virtual photonic lattice (while, as mentioned above, the localized modes with S = 2 are completely unstable in the case of the square lattice) [102].

B. 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 (recall soliton solutions to Eq. ( 38) exist for ω < 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 const > 0, where δ 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 (δ = 0), and repulsion between out-of-phase ones (δ = π).Accordingly, the interplay of the repulsive interaction with the effective Peierls-Nabarro potential, which is pinning the soliton to the underlying lattice [25], 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 [77].
As an example, Fig. 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 (δ = 0), and the attraction leading to instability of the bound state in the case of δ = π.These predictions were also corroborated by numerical findings [77], see an example of a stable bound state of two identical vortex solitons in Fig. 5.

C. 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 [103][104][105].While in most cases gap solitons are predicted theoretically and created experimentally in the context of matter waves [106] and optical pulses [107] 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. [108], 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 [108].
Looking for solutions in the usual form, ψ m,n (t) = exp (−iωt) u m,n , with real frequency (chemical potential) ω, 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 Fig. 5.The stable 2D soliton displayed in panel 5(a) features a typical shape of gaps solitons, with a number of small satellite peaks surrounding the tall central one [104].Bound states of two and four fundamental solitons were found too, featuring weak instability [108].

D. 2D discrete solitons in a rotating lattice
Dynamics of BEC loaded in OLs rotating at angular velocity Ω, as well as the propagation of light in a twisted nonlinear photonic crystal with pitch Ω, is modeled by the 2D version of Eq. ( 1) including the lattice potential, with depth ε and period 2π/k, written in the rotating reference frame: where Lz = i(x∂ y − y∂ x ) ≡ i∂ θ is the operator of the z-component of the orbital momentum, θ being the angular coordinate in the (x, y) plane.In the tight-binding approximation, Eq. ( 42) is replaced by the following variant of the DNLS equation [110]: where C is the intersite coupling constant.In Ref. [110], stationary solutions to Eq. ( 43) were looked for in the usual form (29), fixing ω ≡ −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 Ω, a stability interval for the fundamental soliton, 0 < C < C (fund) max (R), 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 Ω (i.e., their stability interval, 0 < C < C (vort) max (Ω), 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 Ω = 0, are stabilized by the rotation, in an interval 0 < C < C (Ω) ≈ 2.5Ω at small Ω [110].
E. 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 [112].A system of linearly-coupled DNLS equations is a basic model for SSB in discrete settings.Its 2D form is [19] i where ϕ m,n and ψ 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 ω are looked for as (ϕ m,n , ψ m,n ) = exp (−iωt) (u m,n , v m,n ).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 ̸ = B.An example of a stable 2D discrete soliton is displayed in Fig. 7(a), which corroborates accuracy of the VA.In Fig. 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 ≡ E u + E v , see Eq. ( 45). Figure 7(b) demonstrates the SSB bifurcation of the subcritical type [113], 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, Fig. 7(b) demonstrates a bistability area, where symmetric and asymmetric states coexist as stable ones.

III. 2D DISCRETE SOLITONS IN THE SALERNO MODEL (SM)
The 2D version of the SM was introduced in Ref. [118]: Similar to its 1D version (11), Eq. ( 48) conserves the norm and Hamiltonian, cf.Eqs. ( 12) and ( 14), The continuum limit of this model is the 2D equation which is an extension of its 1D counterpart (17): Note that the effective nonlinear-diffraction term µ |Ψ| 2 (Ψ xx + Ψ yy ) in Eq. ( 51) prevents the onset of the collapse because, in the limit of the catastrophic self-compression, this term becomes dominant, giving a positive contribution to the energy.Thus, this term makes it possible to construct stable 2D solitons [118].
2D discrete solitons are looked for as solutions to Eq. ( 51) in the usual form, ψ m,n (t) = e −iωt Φ m,n .In the most interesting case of the competing nonlinearities, µ < 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 Fig. 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 [118].Examples of onsite-and intersite-centered vortex solitons (alias vortex cross and vortex square, respectively) are presented in Fig. 9.In the 2D SM with non-competing nonlinearities (µ > 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 (µ < 0), where unstable vortex modes transform into vortical breathers, i.e., persistently oscillating localized modes that keep the original vorticity.

IV. 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 [119][120][121][122].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 [123].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 ϕ (±) of the pseudo-spinor wave function [124], In this system, SOC of the Rashba type is represented in by terms with coefficient λ, 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 γ, respectively.The system of coupled GP equations (52) maintains 2D solitons of two different types, namely, semi-vortices (SVs) and mixed modes (MMs) [124], The SV solitons, written in polar coordinates (r, θ), have zero vorticity in one component, and vorticity +1 or −1 in the other: where ω is the chemical potential, and f 1,2 r 2 are real functions which take finite values at r = 0 and exponentially decay ∼ (sin (2λr) , cos (2λr)) exp − 2 (−ω − 2λ 2 )r at r → ∞.The two SV solutions ( 53) and ( 54), which are mirror images of each other, exist in the semi-infinite bandgap, ω < −2λ 2 .
The combination of zero and nonzero vorticities in the SV solutions ( 53) and ( 54) 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 ( 54), but they may be approximated by a linear combination of both types of the SVs, ϕ . An essential result is that the SVs and MMs are stable and represent the system's ground state in the cases of γ < 1 and γ > 1, respectively, i.e., when the self-attraction is stronger or weaker than the cross-attraction in Eqs. ( 52) [124].On the other hand, the SVs and MMs are unstable, as excited states, in the opposite cases, i.e., γ > 1 and γ < 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. [125]: The linearized version of this system gives rise to the following dispersion relation for the discrete plane waves, ϕ m,n ∼ exp (−iωt + ipx + iqy), with wavenumbers taking values in the first Brillouin zone, 0 < p, q < 2π: 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 Fig. 10.As concerns the stability, the discreteness extends the stability of the SV and MM solitons towards γ > 1 and γ < 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 λ in Eq. ( 55) exceeds a certain critical value, λ cr .The dependence of λ cr on the soliton's norm, FIG. 11.Left and right: The dependence of the critical value of the SOC strength, λcr, above which the 2D discrete solitons of the SV and MM types, produced by the numerical solution of Eq. ( 55) with γ = 0 and 2, respectively, suffer the delocalization, on the total soliton's norm (57).The figure is borrowed from Ref. [125].
for the SV and MM solitons is displayed in Fig. 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).

V. STABLE SOLITON SPECIES IN THE 3D DNLS EQUATION
A. 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 Sections 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, (l, m, n) 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 [104,105].In that case, ϕ l,m,n (t) is the the respective effectively discretized BEC wave function.
As above, stationary soliton solutions to Eq. ( 58) with chemical potential ω 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 η is a real scale parameter, and it is implied that the vorticity axis is directed along coordinate n [126].
It is also relevant to consider a two-component system of nonlinearly-coupled 3D DNLS equations, for wave functions ϕ l,m,n (t) and ψ l,m,n (t) of two interacting BEC species (most typically, these are different hyperfine states of the same atom) [126]: i φl,m,n + C (ϕ l+1,m,n + ϕ l,m+1,n + ϕ l,m,n+1 Here β > 0 is the relative strength of the inter-component attractive interaction with respect to the intra-component self-attraction.
B. Results

Single-component 3D solitons
The numerical analysis, starting from input (61), has provided families of fundamental and vortex solitons.Here, following Ref.[126], the results are displayed for a fixed value of the chemical potential, ω = −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 (0) cr ≈ 2, and the vortex modes with S = 1 are stable at C < C (1) cr ≈ 0.65.Note that the limit of C → ∞ 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 (1) cr 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 [126].An example of a stable soliton with S = 3 is displayed in Fig. 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 [127].Examples of such states are presented in Fig. 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 Fig. 13( , 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 π).set in parallel planes, as shown in Fig. 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 Fig. 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.Fig. 13(b).Vortex modes with the axis parallel to the 3D diagonal exist too, but they are unstable, see further details in Ref. [127].a. 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 Fig. 15.These bound states are stable for sufficiently small values of the coupling constant, such as C = 0.01 in Fig. 15, and for β < 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 [128][129][130].The discretization of the GP system leads to Eqs. (62) with the opposite sign in front of the nonlinear terms [131].Then, these equations are reduced to stationary ones by the usual substitution with chemical potential ω, {ϕ, ψ} = exp (−iωt) {u l,m,n , v l,m,n }: where the relative strength β of the inter-component repulsion with respect to the self-repulsion remain a positive coefficient.For ω > 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 (|l|, |m|, |n|) → ∞, viz., v 2 background = ω [131].An example of a numerically found skyrmion solution of this type is displayed in Fig. 16.
The same work [131] reported solutions for 2D discrete "baby skyrmions", which are lattice counterparts of the modes produced by the 2D reduction of the Skyrme model [132,133].The have a simple structure similar to its 3D counterpart displayed in Fig. 16, i.e., a complex 2D vortex soliton in one component, embedded into a bubble of the delocalized field in the other real component.

VI. 2D SOLITONS AND SOLITARY VORTICES IN SEMI-DISCRETE SYSTEMS A. 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 [134].It is modeled by the system of linearly-coupled NLS equations for complex amplitudes u n (z, τ ) of optical fields in individual fibers: where z is the propagation distance, τ ≡ t − x/V gr (with time t and carrier group velocity V gr ) is the usual temporal variable, real D is the group-velocity-dispersion coefficient in each fiber, κ > 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 [78]: 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 Fig. 17.This property may be compared to the above-mentioned mobility of 1D discrete solitons in the DNLS equation [72], and of 2D discrete solitons in the framework of the χ (2) 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") [135], as well as soliton-like transient modes with embedded vorticity [136].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 Fig. 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.

B. Semi-discrete quantum and photonic droplets
A new type of semi-discrete solitons was recently elaborated in Ref. [80], 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 [31,35].The system is where ψ j (z) 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-Huang-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 [81]: It corresponds to the array of parallel-coupled planar waveguides, as shown in Fig. 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, or the total power of the photonic droplet, For solitons produced by Eqs. ( 66) and ( 67) sets of control parameters are, respectively, (C, g) for fixed N , or (C, P ).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 Fig. 20.An example of a stable semi-discrete vortex soliton produced by Eq. ( 67) with S = 2 is displayed in Fig. 21.
Getting back to the semi-discrete system (66), a chart in the plane of (C, g) which displays stability areas for the semi-discrete vortex solitons with topological charges S = 2, 3, 4, and 5, is plotted in Fig. 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 χ (2) material with a checkerboard structure representing quasi-phase matching [82].
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 [83].It introduced a χ (2) 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.

VII. 2D FUNDAMENTAL AND VORTICAL DISCRETE SOLITONS IN A TWO-COMPONENT 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 (PT ) symmetry.They include mutually symmetric spatially separated elements carrying linear gain and loss [145][146][147].The experimental FIG.22. Stability areas in the parameter plane (C, g), produced by the numerical solution of Eq. ( 66) for onsite-centered semi-discrete vortex solitons with S = 2 (all colored regions), 3 (orange + brown + green), 4 (brown + green + blue), and 5 (green + blue + dark gray).For the convenience of plotting, the normalizations for S = 2, 3, 4, and 5 are fixed as N = 400, 900, 2500, and 4500, respectively.The figure is borrowed from Ref. [80].
realization of such systems in optics [147] suggests one to include the Kerr nonlinearity, thus opening the way to the prediction an creation of PT -symmetric solitons [148,149].In particular, exact solutions for 1D PT -symmetric solitons and exact results for their stability boundaries were found in the model of the nonlinear PT -symmetric coupler (dual-core waveguide), with mutually symmetric linear gain and loss carried by the linearly coupled cores [150,151].Stability limits for 2D fundamental solitons in the 2D 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. [152].
The definition of the PT symmetry makes it also natural to consider discrete PT -symmetric systems.Various species of stable discrete solitons were predicted in chains of PT -symmetric elements [90,[153][154][155][156][157], and the existence of such solitons was demonstrated experimentally [158].
A natural model for the creation of 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 γ > 0 to the coupled equations [90]: 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 discrete plane-wave solutions to the linearized version of Eqs.(70), As it follows from Eq. ( 71), the PT symmetry holds under condition γ < γ max ≡ 1, i.e., the gain-loss strength γ must be smaller than the linear-coupling coefficient, that is 1 in the present notation, which is a generic property of PT -symmetric couplers [150,151].
Localized states produced by Eqs. ( 72) are characterized, as above, by the total power, Straightforward analysis of Eqs. ( 72) demonstrates that the system produces 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 {u m,n , v mn } = w m,n exp (±(i/2) arcsin γ) , (74) where real discrete distribution w m,n should be found as a solution of the usual stationary equation for 2D discrete solitons, cf. Eq. (38).In agreement with the linear spectrum (71), Eq. ( 75) may produce soliton solutions for k > 1 − γ 2 .An example of a stable fundamental PT -symmetric soliton is displayed, by means of its cross-section shapes, in Fig.

23.
The existence and stability of the PT -symmetric fundamental discrete solitons is summarized in the plane of (γ, P ) for the fundamental solitons in Fig. 24.It is seen that, naturally, the stability area shrinks as the gain-loss coefficient γ is approaching its limit value, γ max = 1 (cf. the 1D situation considered in Refs.[150] and [151]).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 Fig. 24.
The stability boundary in Fig. 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. [19]), 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 [90].In their stability region, the PT -symmetric fundamental discrete solitons actually represent the system's ground state [90].
Alongside the fundamental discrete PT -symmetric solitons, the same system of Eqs.(72) produces PT -symmetric vortex solitons, which also have their stability area, see details in Ref. [90].An example of a stable PT -symmetric vortex soliton is presented in 25.  75) and ( 74), in the plane of the gain-loss coefficient, γ, and total power P , which is defined as per Eq. ( 73).The soliton solutions do not exist in the white area.The figure is borrowed from Ref. [90].72) for (C, γ) = (0.06, 0.4), with propagation constant k = 1 and total power P = 1.65, defined as per Eq. ( 73).The figure is borrowed from Ref. [90].
In addition to the PT -symmetric solitons, Eqs.(72) give rise to anti-PT -symmetric ones, defined by relation v m,n = −u * m,n .They, as well as anti-PT -symmetric vortex solitons, are stable in some parameter areas (see details in Ref. [90]), but those areas are essentially smaller than their counterparts for the PT -symmetric modes.The reduced stability area for the anti-PT -symmetric fundamental solitons is explained by the fact that they cannot be the system's ground state.

VIII. 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 [125] (see also Ref. [137]), sophisticated 3D discrete modes with embedded vorticity [126,127], and discrete skyrmions [131].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 [138] and topological-insulator [139,140] 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 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 [141].Unlike conservative and 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 [142][143][144].

FIG. 1 .
FIG. 1.Comparison of a typical discrete soliton, obtained as the numerical solution of Eq. (30), shown by chains of blue dots, and its counterpart produced by the VA (shown by red open circles).In this figure, ω = −1, see Eq. (29), the corresponding parameters of ansatz(32) being A ≈ 1.31, a ≈ 1.15.The figure is borrowed from Ref.[24].

FIG. 2 .
FIG.2.A stable discrete vortex soliton with topological charge S = 1, produced by Eq. (38) with ω = −3.2.The left and right panels show, respectively, distributions of the absolute value and phase of the complex wave function um,n in the plane with coordinates (m, n).The figure is borrowed from Ref.[93].

FIG. 4 .
FIG. 4.An example of a stable bound state of three 2D fundamental (zero-vorticity) solitons, which are mutually out of phase.The solution is produced by the numerical solution of Eq. (38) with ω = −1.The figure is borrowed from Ref. [77].[This low-quality figure will be replaced in the published version of the paper.]

FIG. 6 .
FIG. 6.(a) An example of a stable 2D discrete soliton with chemical potential µ = 4.05 (in this figure, the notation for the chemical potential is µ, instead of ω, adopted in the text), which corresponds to the right vertical line in (b), found in the mini-gap of the system based on Eqs.(40) and (41) with ∆ = 0.5 and Q = π/3.(b) The dependence P (µ) of the norm of the solitons populating the mini-gap, which is identical to the interval of values of µ represented in the panel.The solitons are stable in the narrow shaded interval.The figure is borrowed from Ref. [108].

FIG. 7 .FIG. 8 .
FIG.7.Left: A stable 2D two-component discrete soliton with spontaneously broken symmetry between the components, generated by system(44).The 2D soliton, with total norm E ≡ Eu + Ev = 1.435, is displayed by means of its cross section.Symbols labelled (UN, VN) and (UA, VA) stand, respectively, for the components of the numerically constructed soliton and its analytical counterpart predicted by the VA based on ansatz(47).Right: Families of 2D onsite-centered discrete solitons, generated by system(44), are shown by means of curves r(E), where r is the asymmetry parameter(46).The dashed-dotted curve shows the VA prediction, while the solid and dashed ones represent stable and unstable solitons produced by the numerical solution.The figure is borrowed from Ref.[19].

FIG. 9 .
FIG. 9. Examples of discrete vortex solitons with topological charge S = 1, produced by the 2D SM, based on Eq. (48).Profiles of the real part of the stationary wave function Φm,n for the vortices of the onsite-centered (stable vortex cross) and intersite-centered (unstable vortex square) types are displayed in the top and bottom panels, respectively.Both solutions are obtained for µ = −0.4 and ω = 7.0.The figure is borrowed from Ref. [118].

FIG. 10 .
FIG. 10.Left: Juxtaposed profiles of ϕ (+) m,n and ϕ (−) m,n (solid and dashed lines, respectively) of a stable 2D discrete soliton of the semi-vortex type, in the central cross section, produced by a numerical solution of Eq. (55) with λ = 0.53 and γ = 0.The soliton's norm (see Eq. (57)) is N = 3.5.Right: The same for a stable discrete soliton of the mixed-mode type, with λ = 0.58, γ = 2, and N = 2. Values of the discrete fileds at lattice sites are connected by lines, for better visualization.The figure is borrowed from Ref. [125].

FIG. 13 .
FIG. 13.Stable 3D multipole solutions of Eq. (60) with C = 0.1 and ω = −2.The top row depicts stable tightly-bound dipoles (with intersite separation d = 1): (a) straight, (b) 2D-diagonal, and (c) 3D-diagonal ones.(d) and (e): Quadrupoles set in the n = 0 plane, with intrinsic separation d = 1 and d = 2, respectively.(f) and (g): Octupoles with d = 1 and d = 2. Panel (h) displays the stability boundary C (3D,d) cr as a function of the intrinsic separation d for diagonal, oblique, and straight dipoles, octupoles, and quadrupoles, from top to bottom.The horizontal dashed line designates the stability threshold for the fundamental discrete soliton.Note that, for the quadrupoles (the bottom boundary), C (3D,d) quad is a linear function of d at d ≤ 4 (see the dashed straight line with slope 0.325, included for the guidance).In panels (a)-(g), level contour corresponding to Re(u l,m,n ) = ±0.5 are shown by blue and red (colors, respectively.The figure is borrowed from Ref. [127].

FIG. 14 .
FIG. 14. Vortex cubes produced by the numerical solution of Eq. (60) with Λ = 2 and C = 0.1.Panel (a) shows a stable composite mode, built of two parallel identical quasi-planar vortices with topological numbers S1 = S2 = 1 and a phase shift of π.Panel (b) shows an unstable vortex-antivortex cube, formed by vortices with opposite topological charges, S1 = −S2.Level contours corresponding to Re(u l,m,n ) = ±0.5 are shown by blue and red colors, and the contours corresponding to Im(u l,m,n ) = ±0.5 are shown by green and yellow colors, respectively.The figure is borrowed from Ref. [127].

FIG. 17 .
FIG.17.An example of a semi-discrete localized spatiotemporal mode, generated by Eq. (64), which performs stable transverse motion under the action of the kick, defined according to(65), with a = 1.5.The cross section of the plot at any fixed z shows the distribution of power |un(τ )| 2 for each n.The figure is borrowed from Ref.[78].

FIG. 19 .
FIG.19.The realization of the semi-discrete system 67: the array of planar optical waveguides (blue slabs), separated by gray isolating layers, with the continuous transverse coordinate, x, and the discrete one, n.As shown by the arrow, light is coupled into the array along the z direction.The figure is borrowed from Ref.[81].

8 FIG. 24 .
FIG.24.Red and gray colors designate, respectively, stability and instability areas for PT -symmetric fundamental discrete solitons, produced by Eqs.(75) and (74), in the plane of the gain-loss coefficient, γ, and total power P , which is defined as per Eq.(73).The soliton solutions do not exist in the white area.The figure is borrowed from Ref.[90].