Two-Dimensional Solitons and Vortices in Linear and Nonlinear Lattice Potentials Two-Dimensional Solitons and Vortices in Linear and Nonlinear Lattice Potentials

Periodic potentials, such as photonic crystals and optical lattices, have shown great abil - ity to manipulate the dynamics of photonic and atomic waves. The interplay of the peri odic potentials and material nonlinearity (self-focusing or defocusing) can create and stabilize several types of solitons, including ordinary and gap solitons, which populate, respectively, in the semi-infinite gap and finite bandgaps of the corresponding linear spectrum. Besides, lattice defects have also been used to construct solitons. This review reports the generation of two-dimensional (2D) solitons in lattice potentials with local defects, under the self-focusing nonlinearity. The numerical analysis demonstrates a novel kind of embedded solitons (or intraband solitons), which are continuous families of 2D localized modes (different from isolated solutions reported before in usual embed -ded-soliton models) embedded into the first and second Bloch bands of the underlying linear spectrum, and pinned to the defect, which determines the spatial position of the modes. We call these modes embedded defect solitons . Further, double, triple, and quadru ple lattice defects can also support stable dipole-mode solitons and vortices. In addition to that, combined linear and nonlinear lattice potentials are also used to construct stable fundamental solitons, gap ones, and solitary vortices. them pseudopotentials in that context), in the form of nonlinear lattices (NLs). The NLs, as the induced pseudopotentials with periodically modulated local strength sign of the nonlinearity. In optics, NLs may be built as material structures with spatially periodic modulations of the local Kerr nonlinearity or other nonlinear coefficients. In terms of BEC, NLs may be engineered via the Feshbach resonance mechanism by using spatially periodic external fields to induce the local nonlinearity [17]. The the soliton dynam ics in NLs, and in NL-LL combinations, have been comprehensively in recent review


Introduction
In their many realizations in various areas of mathematics and physics, solitons are localized "waves of translation" resulting from the balance between linear (dispersive) and nonlinear effects in the medium. For the first time, what is nowadays called as soliton was reported by John Scott Russell as a solitary water wave as early as in 1834. In his experiments, Russel observed that when propagating along a narrow channel, a rounded, smooth, and welldefined heap of water kept its shape over a long distance. This remarkable phenomenon is now called as Russell's solitary wave in fluid dynamics. The soliton phenomenology had started to develop after the advent of modern computers in 1960s. Solitons are, by definition, solitary waves that behave like a "particles," keeping their amplitude, shape, and velocity in the course of propagation and even after colliding with other solitons. Nowadays, soliton phenomena have been identified in various physical systems. For instance, in astrophysics, a robust long soliton behavior appears as a vortex in the rotating atmosphere, which can effectively explain the Great Red Spot in the Jupiter's atmosphere. Rigorously speaking, solitons refer to integrable partial differential equations [1]. In reality, most governing equations concerning soliton phenomena are not integrable [2]. In the latter case, stable pulses are named "solitary waves," rather than "solitons," in mathematical literature, while in physics, the name of solitons is applied to robust pulses in nonintegrable models equally well. In fact, they can readily demonstrate dynamics similar to that in integrable systems over experimentally relevant timescales [3]. In this brief review, we adopt the term soliton to refer to localized waves supported by the balance of linear dispersion and nonlinearity.
Theoretical and experimental studies of the solitons dynamics have drawn a great deal of interest in the course of several past decades. It is now known that solitons are universal collective excitations in very diverse areas, ranging from mathematics to physics and chemistry, and even to biology. In physics, the studies of solitons chiefly refer to self-trapped light pulses and beams in nonlinear optics, strain waves in elastic media, matter waves in superfluids (especially, solitons in Bose-Einstein condensates-BECs), various species of magnetic solitons, and other nonlinear waves in magnetic media, etc. The studies in these areas are focused on the existence, stability, mobility, and interactions of solitons. They are most often found in one-dimensional (1D) physical systems, although they exist in multidimensional settings too. In the two-dimensional (2D) case, the ubiquitous cubic self-focusing readily gives rise to formal solutions for solitons. However, the same setting brings about the critical collapse [4]; therefore, fundamental solitons in the free 2D space (Townes solitons [5]) are unstable, as a separatrix between collapsing and decaying modes. Vortex solitons exist too in the 2D free space [6], while being vulnerable to a stronger splitting instability against azimuthal perturbations [7,8]. Vortex solitons are a special type of solitons containing screw phase dislocations. To stabilize solitons in such physical systems, researchers have turned to adding different types of periodic potentials, alias linear lattices (LLs) to the systems under the consideration [9].
Suggested by the ability of crystals to mold the flow of electrons in solids, different periodic potentials, such as ones induced by photonic crystals and lattices in optics [10], gratings built into plasmonic waveguides [11], optical lattices in BECs [12][13][14] or degenerate Fermi gases [15], etc., were used to gain control of the dynamics of photonic, plasmonic, and atomic waves, respectively. Periodic potentials, when combined with the self-focusing or defocusing material nonlinearity, have been shown as excellent candidates for stabilizing various types of solitons [8,9,13,14,16], including ordinary ones, lying in the semi-infinite gap (SIG) of the corresponding linear spectrum, and gap solitons (GSs) populating finite bandgaps.
Periodic potentials have recently been extended to make them nonlinear (sometimes called pseudopotentials in that context), in the form of nonlinear lattices (NLs). The NLs, as the names suggest, induced pseudopotentials with periodically modulated local strength and/ or sign of the nonlinearity. In optics, NLs may be built as material structures with spatially periodic modulations of the local Kerr nonlinearity or other nonlinear coefficients. In terms of BEC, NLs may be engineered via the Feshbach resonance mechanism by using spatially periodic external fields to induce the local nonlinearity [17]. The studies of the soliton dynamics in NLs, and in NL-LL combinations, have been comprehensively summarized in recent review [18].
This review aims to provide a brief survey of our results reported in the context of theoretical studies of 2D solitons and vortices in linear and nonlinear lattice potentials [19,20]. In Section 2, we introduce a dynamical model of lattice potentials with defects and self-focusing nonlinearity and report the existence and stability of a new type of 2D embedded (alias intraband) solitons, which we call as embedded defect solitons (EDSs), as they are embedded, as localized modes, into the first and second Bloch bands of the underlying linear spectrum, being pinned to local lattice defects. Stable dipole-mode solitons and vortices supported by multiple defects are also reported in this section. The variational approximation (VA) based on the Gaussian ansatz, along with numerical methods, is presented to study 2D ordinary solitons in the setting combining the LL and NL with incommensurate periodicities of the two lattices in Section 3. GSs and vortex solitons in the same setting are studied too by means of numerical simulations. Finally, we formulate conclusion in Section 4.

Introduction
As mentioned above, the model with the integration of periodic potentials and material nonlinearity (self-focusing or defocusing) can give rise to various types of solitons, such as ordinary solitons and GSs, existing, respectively, inside the SIG and the finite bandgaps of the corresponding linear spectrum.
Furthermore, it was also demonstrated that, counter-intuitively, specific species of solitary waves, known as embedded solitons, may exist inside spectral Bloch bands populated by radiation modes [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. Generally, embedded solitons cannot exist in continuous families because of the resonance with radiation waves. However, isolated embedded modes are possible, as the rate of the decay into radiation may vanish at particular spectral points. Under some specific conditions, continuous families of embedded solitons were constructed too [30][31][32][33]. Recently, the formation of such solitons was demonstrated in a model combining the quadratic (χ (2) ) nonlinearity and a complex lattice potential, whose real and imaginary parts are subject to the condition of parity-time (PT) symmetry, representing symmetrically placed and mutually balanced local gain and loss [34]. While embedded solitons have been widely investigated in general 1D models, it remains an open issue to discover them inside 2D Bloch bands. Previously, stable intraband solitons have been found in a 2D potential which consists of a periodic lattice in one direction and a harmonic-oscillator trap applied in the orthogonal direction [35].
Besides the perfect lattices, lattice defects have also been introduced to help constructing solitons. In particular, solitons supported by local defects in the form of the localized χ (2) nonlinearity have been studied in [36][37][38]. Recently, solitons pinned to local PT-symmetric defects, inserted into the 1D medium with the uniform χ (2) nonlinearity, were considered too [39]. Along these lines, the current section aims to report the numerical analysis of 2D solitons in lattice potentials with local defects, under the action of self-focusing nonlinearity. Numerical results demonstrate a novel kind of embedded solitons (or intraband solitons), which form continuous families of 2D localized modes (different from isolated solutions reported before in general embedded-soliton models), existing inside the first and second Bloch bands of the underlying linear spectrum, in terms of the propagation constant, and pinned to the defect. We call them as embedded defect solitons (EDSs) accordingly. It is necessary to point out that, for localized modes pinned to local lattice defects, one can define the location of their propagation constant relative to the underlying bandgap spectrum of the infinite uniform LL, as the latter is defined as the spectrum of radiation waves into which the localized mode may decay.
Since the EDSs exist with the norm exceeding a finite minimum value (and the propagation constant of linear defect modes cannot be located inside a Bloch band), they do not emerge from a continuation of linear localized defect modes. Systematic direct simulations verify that the EDSs are totally stable in the first Bloch band, and partially stable in the second. Along with the regular defect solitons indwelling the SIG, and gap defect solitons (GDSs) localized inside the first finite bandgap (the spectral interval between the first and second bands) and pinned to the same defect, the intraband EDSs constitute a continuous "superfamily" throughout the overall bandgap structure. It is important to identify the stability of the GDSs under the self-attractive nonlinearity, since previous results for stable gap solitons supported by perfect 2D lattices were limited to the case of self-repulsion, while the case of self-attraction was assumed to be always unstable [8,40,41]. Although GSs bifurcating under the self-attractive nonlinearity from edges of adjacent Bloch bands into the first finite bandgap have been constructed in [42], their stability has not been considered there. In addition, 2D fundamental GSs under the self-attractive nonlinearity were found to be completely unstable in the first finite bandgap, while a family of dipole-mode gap solitons is stable in a part of the first bandgap as long as the depth of the lattice potential takes values above a certain threshold [43].

Model
The model is the 2D mean-field Gross-Pitaevskii (GP) equation for the BEC wave function (or the nonlinear Schrödinger equation for the amplitude of the electromagnetic wave), ψ(x, y, z), (in BEC, propagation distance z is replaced by time t), with the self-focusing cubic nonlinearity, and ∇ 2 ≡ ∂ x 2 + ∂ y 2 . The perfect lattice potential with depth 2ε and period π is taken in the general form, V OL = ε [ cos ( 2x ) + cos ( 2y ) ] . The defect is defined by δ ( x, y ) = 1 at x 2 + y 2 < r 0 2 and δ ( x, y ) = 0 at x 2 + y 2 ≥ r 0 2 in Eq. (1), with a "hole" of radius r 0 = 1.2 in the potential. Numerous simulations demonstrate that, for this shape of the defect with fixed ε = 2, fully stable solitons can only be found in the interval of 1.1 ≤ r 0 ≤ 1.3 , for different shapes of the "hole." The contour plots of the LL potential with the single or compound defects are shown in Figure 1.
Stationary solutions to Eq. (1) with propagation constant -μ (or chemical potential μ, ∈ terms of BEC ) are sought for as ψ = ϕ ( x, y ) exp ( − iμz ) , with amplitude ϕ ( x, y ) obeying the stationary equation, It should be noted that the fundamental and dipole-like solitons are represented by real solutions for ϕ ( x, y ) , while the vortical counterparts (vortices that can be supported by quadruple defects, as shown below) correspond to complex solution.
The linear spectrum of Eq. (1), obtained by a numerical solution of the linearized version of Eq. (2), is displayed in Figure 2. It can be seen that, for the given parameters, a growing number of defect-induced isolated eigenvalues emerge in the second finite bandgap with the increase in the number of the defects (single → double → triple), although the first bandgap remains unchanged. In the following Sections 2.3 and 2.4, rather than discussing apparent  quasi-linear defect modes generated from these isolated eigenvalues, we will analyze localized states that have no linear limit but are supported by the defects, in the first finite bandgap and in the two Bloch bands adjacent to it. Numerically, such states were constructed as solutions of Eq. (2), using the Newton's method in the domain of size 30 × 30, with a spatial grid of 192 × 192 points. The initial wave was the usual Gaussian, with amplitude A, width W, and total power (or number of atoms, in the BEC) Stability of the solutions was then checked in direct simulations of the perturbed evolution in the framework of Eq. (1).

Numerical results for the single defect
The existence of continuous families of localized modes (fundamental EDSs), pinned to the single defect (the one shown in Figure 1a), is illustrated in Figure 3. The numerical solution proves that the stability region of EDSs lies inside the first and second Bloch bands and is connected to the family of GDSs that are pinned to the same defect, but with the propagation It can be seen from Figure 3 that the combined family of the pinned solitons exists above a finite minimum value of the total power, N min = 1.56, which exactly corresponds to the junction point between the first Bloch band and the first finite bandgap. This implies, as mentioned above, that the family of localized modes (the regular solitons, EDSs and GDSs) supported by the single defect does not appear as a nonlinear evolution of any linear defect mode because the latter would lead to N→0.
Numerical calculations prove that the EDS is totally stable within the first Bloch band, while being only partially stable in the second band. A noteworthy feature is that, in the first Bloch band, the stability of the intraband solitons pinned to the defect satisfies the Vakhitov-Kolokolov (VK) criterion, which is relevant for all kinds of localized modes supported by a self-focusing nonlinearity [44], ∂ μ / ∂ N < 0 . On the other hand, the EDS branch features ∂ μ / ∂ N > 0 in the second Bloch band, where the stability region is very narrow (actually, it may be a region of weak instability, which is, however, equivalent to stability in a possible experiment). Of course, there is no rigorous proof of the applicability of the VK criterion to the current model. Shapes of stable EDSs inside the first and second Bloch bands (at points marked A and B in Figure 3), and an unstable EDS at point C, are depicted in Figure 4. In the whole first band, the form of the EDSs is similar to that shown in Figure 4a. Furthermore, in the second band (as shown in Figure 4b and), the localized modes feature a more complex structure, with pronounced side peaks, in line with the principle that the shape of the Bloch modes is more complex too in the same band. As a matter of fact, the expanding tail of the soliton displayed in Figure 4c initiates the onset of the instability of this soliton. Direct simulations demonstrated that the unstable modes decay into radiation (not shown here in detail).
The current model can support a new type of nonlinear localized states found inside Bloch bands of the spectrum induced by the 2D lattice (2D solitons of the embedded type were found in Ref. [35], but the lattice was one-dimensional in that case). Actually, the EDS family traversing the first Bloch band serves as a missing link between the two famous kinds of solitons created by the lattice potential, which are regular solitons nested in the SIG, and gap solitons populating the first finite bandgap (in the current case, both of them are pinned to the single defect). The GDS family was also found in the first finite bandgap, being linked to the intraband EDS solutions between the two nearby Bloch bands (as in Figure 3). Our simulations have confirmed that these gap solitons, also pinned to the defect, are stable under the self-focusing nonlinearity [defined as in Eq. (1)], in contrast to their counterparts in ideal lattices, where gap solitons are generally considered to be stable only under the self-defocusing nonlinearity [8,40].
Furthermore, the GDSs were also explored under the self-defocusing nonlinearity, and it was concluded that they are always unstable (not displayed here in detail). The reason is that the effective mass is negative for the gap soliton under the self-defocusing nonlinearity, therefore inverting the feature of the soliton-defect interaction, and rendering unstable the localized state of the gap soliton in the case of the attractive defect. Additional analysis strengthens this argument: if the local defect was replaced by one with the opposite sign, the pinned states of gap solitons were demonstrated to be stable under the self-defocusing nonlinearity and unstable in the case of self-focusing. The negative mass of gap solitons explains a number of other counterintuitive dynamical effects featured by these modes [45][46][47][48]. Figure 5 presents a typical example of the stable pinned GDS, corresponding to the point D in Figure 3, found in the first finite bandgap and supported by the self-focusing nonlinearity.
Our simulations further suggest that stable GDSs pinned to the defect cannot be created in the relatively narrow second finite bandgap (see Figure 2), in agreement with the general trend of gap solitons to be unstable in the second bandgap [13,14,49].

Numerical results for multiple defects: stable dipoles and vortices
The above analysis was also extended to the cases of double, triple, and quadruple lattice defects, such as those displayed in   The triangular defects (see Figure 1c) are able to support stable triangularly shaped dipole modes too, see examples displayed in Figure 7. Typical vortices supported by lattice potentials can be constructed as orthogonal forms of four intensity peaks [8,9] and an almost hollow site at the pivot (vortices of this kind are usually called onsite-centered ones), featuring the phase differences of π/2 between adjacent peaks, which amounts to the total phase circulation of 2π (i.e., topological charge 1). The quadruple defects (see Figure 1d) can create such stable solitary vortices, characteristic examples of which are shown in Figure 8, while the entire "superfamily" is presented in Figure 9. As indicated in Figure 9, direct simulations demonstrate that the pinned vortices are unstable in the second Bloch band, being stable elsewhere.

Introduction
Besides the LL potentials, NLs have also been proposed to support localized modes in recent years. In optics, the photonic-crystal structures can induce an effective combined LL-NL potential, formed by congruent periodic spatial modulations of the linear refractive index and local nonlinear Kerr coefficient. In terms of BEC, the NLs can be introduced using external magnetic fields that affect the nonlinearity coefficient via the Feshbach resonance. Results produced by studies of localized modes (solitons and vortical ones) in NLs and LL-NL combinations have been reviewed in [18]. Experimentally, NL-supported optical solitons have been created at an interface of two lattices [50].
A natural consideration of the combined LL-NL system is the one with different or incommensurate periodicities between the two lattices. In the 1D model, both numerical simulations and analytical approximations have demonstrated that such a combination of competing linear and nonlinear potentials can support various localized modes, both ordinary solitons and GSs [51].
In particular, intriguing results were obtained for existence borders of the solitons as functions of the LL-NL incommensurability and for the empirical "anti-Vakhitov-Kolokolov" (anti-VK) stability criterion for GSs, which is expressed as the dependence of the chemical potential, μ, on the norm, N, of the soliton: dμ/dN > 0 (ordinary solitons supported by the self-focusing nonlinearity in the SIG obey the VK criterion, dμ/dN < 0 [44]). Here, we focus on the case of 2D solitons supported by the incommensurate LL-NL combinations. A VA approach and numerical methods are applied to ordinary solitons, while GSs in the first finite bandgap are studied only numerically, as the analytical consideration would be too cumbersome. Vortex solitons, in the semi-infinite and finite gaps alike, are constructed in a numerical form too.
We study the physical settings including both full 2D and quasi-1D (Q1D) lattice potentials, the latter depending on the single coordinate (the Q1D LL potential is sufficient to stabilize 2D ordinary solitons in diverse 2D self-focusing media [52][53][54]). Actually, the combination of the incommensurate LL and NLs (with different periods) makes the medium an effective quasi-crystal for nonlinear excitations. Theoretical work on fundamental solitons and solitary vortices has been done for linear quasi-periodic potentials [55], and experiments with 2D photonic quasi-crystals counterpart have been reported recently [56,57]. Another relevant work [58] dealt with 2D solitons supported by combined crossed Q1D linear and nonlinear periodic potentials.

The model
The 2D theoretical model combing the periodic LL potential and NL pseudopotential is based on the scaled Gross-Pitaevskii (or nonlinear Schrödinger) equation for the BEC wave function, or the local amplitude of the electromagnetic wave in optics, ψ(x, y, z): where z is propagation distance (or time t in BEC), ∇ 2 = ∂ x 2 + ∂ y 2 , the LL period is set to be 1, the period of the NL is 2/q (q is the incommensurability index), and the NL strength is normalized to g = ±1 . The center of the soliton will be fixed at point x = y = 0 , thus g = +1 and −1 pertain, severally, to the self-attraction and self-repulsion nonlinearity, which support ordinary solitons in the SIG or GSs in finite bandgap(s), respectively.
The remaining parameter in Eq. (3) is the LL strength, ε. The relevant band structure, for a fixed strength, ε = 7.4, in the reduced Brillouin zone [1], produced by linearizing Eq. (3), is shown in Figure 10.
Eq. (3) corresponds to the general model with the LL potential and spatially constant nonlinearity at q = 0. The model is commensurate, with respect to the LL and NL, at q = 2, and subharmonically commensurate once the ratio between the LL and NL periods is 1:2 at q = 1. The model features the full incommensurability (the system is totally quasi-periodic) when q becomes an irrational number, while, actually, the incommensurability is emulated by setting q = 1.5, which yields the period ratio 3:4. The Q1D model corresponds to discarding the linear potentials cos(2πy) and/or nonlinear terms cos(qπy).
Stationary solutions of Eq. (3) with chemical potential μ (or, in terms of optics, propagation constant −μ) can be looked for as ψ = ϕ ( x, y ) exp ( − iμz ) , with wave function ϕ ( x, y ) obeying the stationary equation, The VA is based on the Lagrangian of Eq. (4), which is

Fundamental solitons
The ordinary solitons, which are expected to form in the SIG under the action of the selffocusing nonlinearity ( g = +1 ) in Eq. (3), are first sought by means of the VA based on the Gaussian ansatz, ϕ ( x, y ) = Aexp [ − ( x 2 + y 2 ) / ( 2 W 2 ) ] , with the corresponding norm N = ∬ | ϕ ( x, y ) | 2 dxdy π A 2 W 2 . Substituting this ansatz into the Lagrangian produces the following result, expressed in terms of N and width W: and the respective variational equations, ∂ L eff / ∂ W = ∂ L eff / ∂ N = 0 :  Typical shapes of stable solitons created inside the SIG by the full 2D model, as well as by its version with the Q1D linear potential, are shown in Figures 12 and 13. Obviously, these shapes correspond, respectively, to quasi-isotropic and strongly elongated localized modes, resembling those found previously in other 2D models stabilized by LL potentials [8,9,18,58].
It can be seen from Figure 11 that the VA predicts the ordinary solitons with a reasonable accuracy, except near the edge of the SIG, where the Gaussian ansatz is irrelevant, due the complex shape of the soliton. Furthermore, direct simulations of the dynamics of disturbed solitons prove that the solitons' stability correctly follows the VK criterion, dμ/dN < 0 (the dependence μ(N) is obtained numerically). These properties for the family of ordinary solitons resemble those reported previously in the 1D combined LL-NL model [51]. Figure 14 summarizes the stability of the 2D ordinary solitons in the SIG for all four realizations of the model (2D or Q1D linear and/or nonlinear potentials). The unstable solitons suffer decay into radiation waves (rather than transforming into stable solitons). A notable feature is that replacing full 2D NL by its Q1D version results in expansion of the stability areas, which can be explained by the fact that, in the Q1D NL, ordinary solitons attempt to "elude" the destabilizing locally self-repulsive nonlinearity only in one direction, rather than in two.
In addition, we tried to test the mobility of the solitons by studying their evolution initiated by a sudden "kick," that is, multiplication of the wave function of a stable quiescent soliton by the phase-tilt factor, exp [ i ( k x x + k y y ) ] , with the vectorial kick parameter k. The solitons were found to be immobile under the action of the Q1D LL or NL potential, except for an obvious case when both the LL and NL have a collinear Q1D structure, while the kick is applied in the unconfined direction. A sufficiently strong kick would destroy the soliton, instead of setting it in motion. This type of solitary vortices is generally stable, owning to the weak interaction between the peaks. More densely packed vortex patterns can be constructed too, but they were found to Figure 13. The same as in Figure 12, but in the case of the Q1D linear potential (combined with the 2D nonlinear potential), for g = +1, q = 1, μ = −2.91, and N = 1.68. The shape of the solitons in the model with the quasi-1D nonlinear potential is similar to the present one. be unstable. Studies of other models have also proved that the vortices with inner "voids" are more likely to be stable [8,18].
Numerical simulations have verified that stable vortices may only be created for values of incommensurability index close to q = 0, 1, and 2, namely, within intervals of the half-width Δq ≈ 0.1 around these values. The latter observation makes sense, as both the linear and nonlinear potentials have minima at or close to sites where the power (density) peaks are located, at such values of q. The stable vortices of the type shown in Figure 15 can only be observed (in the SIG) under the fully 2D LL potential, and the NL, however, can be taken in 2D or Q1D form. Figure 16 shows μ(N) curves for the families of the vortex modes. The analysis proves that the stability of such solutions totally obeys the VK criterion, that is, the families with dμ/dN < 0 are stable.    Figure 17 shows an example of stable GS generated sufficiently far from edges of the finite bandgap (the shapes of the GSs are quite similar in the model combining the 2D LL and Q1D NL). It is seen that the numerically found fundamental GSs feature, as usual, is in more complex shapes than that in the ordinary solitons. The μ(N) curves for the GS families are shown in Figure 18a. Unlike the ordinary solitons (cf. Figure 12), GSs always feature dμ/dN > 0, thus complying with the anti-VK criterion [51].

Fundamental solitons
Direct simulations demonstrate that the GSs become stable sufficiently deep inside the finite bandgap and unstable near its edges (unstable GSs suffer decay into radiation). The numerically found stability borders for the GSs in the full 2D NL and Q1D versions of the current model are displayed in Figure 18(b). It is observed that, ion the contrary to the ordinary solitons in the SIG, the stability region for GSs becomes extremely narrow for the Q1D NL, compared to the case of the full 2D NL. The latter feature seems obvious, since, in contrast to the ordinary solitons, and the NL may provide necessary support to the GSs.
The GS stability areas gradually shrink to zero with an increase in q, as can be seen in Figure  18b for the variant of the Q1D NL model (we envisage similar situation for the full 2D NL model, while numerical troubles stem from expanding the stability diagram at larger values of q). This tendency can be understood due to the fact that the fast oscillating NL field averages itself to zero at large q; thus, the broad (see Figure 17) GS cannot feel the action of the nonlinearity. The ordinary solitons within the SIG do not follow this trend (cf. Figure 14) because, as q increases, these solitons can shrink in a single cell of the structure, remaining confined to a self-attractive nonlinearity region. Lastly, numerical simulations demonstrate that, similar to what was found above for the ordinary solitons, GSs are immobile objects too (not shown here in detail).  Figure 19 displays an example of stable vortex solitons created in the first finite bandgap. Such vortex soliton features the form similar to its counterparts in the SIG, cf. Figure 15, namely, it is composed of four peaks which is separated by twice the LL period, while the phase distribution carries the vorticity. Opposite to the case in the SIG, the solitary vortices inside the finite bandgap may be stable only when both LL and NL are taken in the full 2D form (i.e., the vortices are unstable if the Q1D NL is used).

Solitary vortices
Like their counterparts in the SIG, stable solitary vortices in the finite bandgap are found as long as the incommensurability index is close to q = 0, 1, and 2. Families of these vortices are displayed in Figure 20 by means of the μ(N) curves. As in the case of the fundamental GSs, the stable vortices in the finite bandgap obey the anti-VK criterion, dμ/dN > 0.

Conclusion
We have studied localized states in the framework of the nonlinear Schrödinger equation by introducing single, double, triple, and quadruple defects into the 2D lattice potentials in the self-focusing medium. The model can be realized as lattices in nonlinear photonic media, and BECs trapped in the OL. A new type of 2D embedded (alias intraband) solitons, which we call embedded defect solitons (EDSs), has been constructed. These modes are stable in the first Bloch band and partially stable in the second band. Contrary to the fact that 2D fundamental gap solitons may be stable only under the self-repulsion in uniform lattices, stable gap solitons pinned to the defect are found under the self-attractive nonlinearity. Stable dipole-mode solitons and vortices maintained by multiple defects have been identified too.
Localized states supported by the model of 2D nonlinear photonic crystals or BEC under the effect of linear and nonlinear lattices with different (generally incommensurate) periods have also been investigated. The combined incommensurate linear and nonlinear lattices may be considered as a "nonlinear quasi-crystal." Both fully 2D periodic potentials and their Q1D reductions were considered, and the VA (variational approximation) was developed for ordinary (fundamental) solitons in the semi-infinite gap. The stability regions for the entire sets of ordinary solitons and gap ones as well as the vortex solitons (built as "hollow" four-peak complexes) of both types have been produced. New directions can be naturally explored in the future. Introducing defects in complex lattices, whose real and imaginary parts obey the PT symmetry, one may consider the stabilization of localized modes and, in particular, of 2D gap solitons that were found to be unstable in uniform PT-symmetric lattices [59]. Search for stable vortex solitons with higher topological charges and, possibly, "supervortices" (vortex rings built of compact localized vortices) [60,61] may be relevant too. On the other hand, it may be interesting to extend the analysis to higher-order bandgaps, aiming to construct solitons and solitary vortices in those gaps.