Stable two-dimensional soliton complexes in Bose-Einstein condensates with helicoidal spin-orbit coupling

We show that attractive two-dimensional spinor Bose-Einstein condensates with helicoidal spatially periodic spin-orbit coupling (SOC) support a rich variety of stable fundamental solitons and bound soliton complexes. Such states exist with chemical potentials belonging to the semi-infinite gap in the band spectrum created by the periodically modulated SOC. All these states exist above a certain threshold value of the norm. The chemical potential of fundamental solitons attains the bottom of the lowest band, whose locus is a ring in the space of Bloch momenta, and the radius of the ring is a non-monotonous function of the SOC strength. The chemical potential of soliton complexes does not attain the band edge. The complexes are bound states of several out-of-phase fundamental solitons whose centers are placed at local maxima of the SOC-modulation phase. In this sense, the impact of the helicoidal SOC landscape on the solitons is similar to that of a periodic two-dimensional potential. In particular, it can compensate repulsive forces between out-of-phase solitons, making their bound states stable. Extended stability domains are found for complexes built of two and four solitons (dipoles and quadrupoles, respectively). They are typically stable below a critical value of the chemical potential.

can be seen as a result of a counterflow produced by anomalous spin-dependent velocity determined by SOC [30]. These SVs and MMs are stable, severally, in the system with XPM stronger than SPM, and vice versa. Both SV and MM species are stable in the case of equal strengths of the SPM and XPM terms (the Manakov's nonlinearity [31]). When the SV or MM is stable, it represents the system's ground state (in the absence of SOC, the ground state is missing in the unstable BEC), provided that the total norm of the pseudo-spinor wave function does not exceed a critical value (in fact, it is the norm of the Townes soliton [32], which is an unstable solution of GPE in 2D in the absence of SOC [19,20]). Above the critical norm, where the centrifugal counterflow, induced by the anomalous velocity [30], cannot overcome the attraction-induced density flux, the collapse still takes place, albeit being modified by the SOC. In [18] it was shown that SOC stabilizes solitons, gap-solitons and half-vortices of different symmetries in Zeeman lattices. Later, similar results for stable 2D solitons of the SV and MM types were reported for the SOC of a more general type, combining the Rashba and Dresselhaus couplings [33], as well as for the binary BEC of "heavy atoms", for which the kinetic energy (second derivatives in GPE) may be neglected in comparison with the SOC energy, represented by first-order cross-derivatives in the GPE system [34]. Also demonstrated [35] were the stabilizations of MM states by SOC in a system including additional nonlinearity which represents the beyond-mean-field corrections of the Lee-Huang-Yang type [36,37,38,39,40,41,42]. Eventually, the creation of metastable solitons of both SV and MM types was predicted in the 3D system with SOC [43]. Furthermore, it was demonstrated that the SOC with reduced dimensionality, viz., oneor two-dimensional Rashba coupling acting in the 2D or 3D space, is sufficient for the stabilization of the 2D and 3D solitons, respectively [44]. The SOC acting in a spatially confined area is also sufficient for maintaining stable 1D [14] and 2D [45] solitons.
More generally, the stabilizing effect of gauge fields is known even since earlier times. In [46] it was shown that stable solitons in one-component 2D BECs can be sustained by the rotating trap emulating a gauge field. Confinement of a spinor BEC in the quasi-relativistic limit by a gauge field was reported in [47]. Recently, in [48] it was discovered that field components, taken in a pure or non-pure gauge form, affect the soliton dynamics in very different ways. Localization characteristics of nonlinear states are determined by the curvature originating from the non-pure gauge field, while the pure gauge field strongly affects the stability of the modes. The respective solutions represent envelopes formed by a non-pure gauge field, which modulates stationary carrier states created by a pure gauge field. In coupled 1D nonlinear Schrödinger equations with SU(2)-invariant nonlinearity, a non-uniform pure gauge field preserves the integrability of the system [16]. The total effect of a SOC gauge field can lead not only to the enhanced stability of solitons in media with attractive interactions, but also to the existence of novel stable localized modes -quasisolitons -in purely repulsive condensates, both in the free space and in the presence of traps [48].
In addition to the BEC realm, SOC may be emulated in terms of nonlinear optics. In that connection, stabilization of 2D spatiotemporal optical solitons in a dual-core planar waveguide with the help of the counterpart of the BEC effect, which is represented by temporal dispersion of the inter-core coupling, was predicted in Ref. [49].
Although SOC appears to be a universal mechanism for the stabilization of 2D solitons, its efficiency is limited. It cannot stabilize excited states, obtained by adding equal vorticities to both components of SV and MM solitons in attractive media [29], although stable quasisolitons with embedded vorticity in both components in media with repulsive cubic interactions have been found [48]. Possibilities of maintaining bound states of separated solitons by SOC are not yet known either. For this reason, a relevant possibility, which we pursue in this work, is to consider SOC in the 2D geometry, subject to spatially periodic modulation (solitons and their interactions in a 1D BEC with a helicoidal SOC were considered in [15]). In the experiment, this setting can be implemented by means of an appropriate OL, which determines the spatially periodic SOC structure. As we demonstrate below, this system, with a helicoidal orientation of the local SOC [see equation (1)], gives rise to a vast stability region for fundamental 2D pseudo-spinor solitons. Further, it also produces stable dipole and quadrupole bound states of fundamental solitons with opposite signs, which cannot be found in the absence of the spatially-periodic modulation. We also demonstrate that stable solitons form sets of states, generated from a seed one by symmetry transformations.
The rest of the paper is organized as follows. The pseudo-spin system of coupled GPEs is formulated in Section 2, where we also produce its linear spectrum, in the numerical form in the general case, and analytically in the long-wave limit. Symmetries of the system are considered in Section 3. Systematic numerical results for fundamental solitons, dipoles, and quadrupoles, including the analysis of their stability and symmetry, are presented in Section 4. The paper is concluded by Section 5.

The model and linear spectrum of the system
We describe the evolution of the 2D spinor BEC in the presence of SOC by coupled GPEs for a spinor wave function Ψ = (Ψ 1 , Ψ 2 ) T : Here σ = (σ 1 , σ 2 , σ 3 ) is the vector of the Pauli matrices, σ 0 is the identity matrix, α is the strength of SOC and the position-dependent unit vectors n and m determine the SOC direction. We consider a helicoidal SOC with n = (cos θ(x, y), 0, sin θ(x, y)) and m = (− sin θ(x, y), 0, cos θ(x, y)). The above choice of vectors m and n yields the components of a non-Abelian spin-dependent gauge field with a position-independent commutator, In what follows, we concentrate on a spin-orbit coupling realized with a lattice symmetry. In this choice the coordinate-dependent phase θ(x, y) is given by with periods 2π/λ x,y , as shown in figure 1.  (3), with λ x = λ y = π. Middle and right panels: plots of cos θ(x, y) and sin θ(x, y), which determine directions of local vectors n and m.
The system includes the attractive nonlinearity with equal strengths of the selfand cross-interactions in equation (1). We stress that no Zeeman splitting, which may stabilize solitons in other SOC models [15,48], or external potential is present in equation (1). Although the stabilization of fundamental solitons in the present setting may be provided by the uniform SOC [29], it is demonstrated below that the spatially periodic modulation of local SOC is necessary to create stable multisoliton complexes (dipoles and quadrupoles), which do not exist in the uniform space.
The model (1), (3) has three length parameters. Two of them, 2π/λ x,y , are the lattice constants. The third one, 1/α, usually referred to as the "spin-flip length" and characterizing the coordinate dependence of a spinor describing spin 1/2 particle [3] can, in general, be scaled out, i.e., the dynamics is governed by λ x,y /α ratios. Nevertheless, for the analysis of different limiting cases it is convenient to keep α and λ x,y as separate parameters.
Before analyzing nonlinear states supported by equation (1), it is instructive to consider the spectrum of linear modes determined by Hamiltonian H, which represents the kinetic energy H kin and spin-orbit coupling H soc in equation (1) in the form: Here has the spin-diagonal form. The other contribution, H soc , can be presented as the sum of two terms: defined as: and Hereafter we omitted the explicit spatial dependence where it is not required directly. As coefficients in equation (1) are periodic functions of x and y, we seek for eigenfunctions in the form of Bloch waves where r = (x, y), C ↑,↓ (x, y) = C ↑,↓ (x + 2π/λ x , y) = C ↑,↓ (x, y + 2π/λ y ) are components of a time-independent spatially-periodic spinor wave function, k x,y are the Bloch momenta in the first Brillouin zone, µ lin (k) is the chemical potential of the linear system, and we omitted the Bloch band index since we are interested only in the lowest bands. The substitution of the Bloch ansatz (9) in the linearized version of equation (1) leads to an eigenvalue problem, which determines the bandgap spectrum of the system, µ lin (k). We begin with a particle in the free space, where we can set λ x = λ y = 0 and θ = π, the coefficients C ↑,↓ (x, y) are position-independent, and the values of k x,y are not restricted. According to equation (10) the chemical potential in this case becomes µ lin (k) = α 2 +k 2 /2±α k, where k ≡ (k 2 x +k 2 y ) 1/2 and the spinor components, for different values of µ lin (k), satisfy conditions where upper/lower sign corresponds to the sign in front of αk in µ lin (k).
In what follows we concentrate on a square lattice with parameters λ x = λ y ≡ λ and begin with displaying two numerically calculated lowest bands of the linear spectrum in figures 2(a,b). The solitons of our interest, formed by the attractive nonlinearity, belong to the semi-infinite gap with the values of µ below the bottom of the lowest band. As one can see from two cuts of the linear dispersion relation, shown in figures 2(a,b), the minimum of µ lin (k) is achieved on the ring of radius k min in the (k x , k y ) plane. The radius of this ring, being a non-monotonous function of α, shown in figure 2(c) by the black curve, may shrink to zero at certain values of α, while it monotonously increases at a sufficiently strong SOC.
For small Bloch momenta, |k| λ, and small SOC strength, α λ, one can find the spectrum in an approximate analytical form by means of averaging with respect to relatively rapid oscillations of cos θ(x, y) and sin θ(x, y), while keeping the slowly oscillating function exp (ikr) in ansatz (9). A nontrivial result can be obtained in the asymptotic regime defined by the following relation between the two small parameters: k ∼ α. Under this condition, in the zero-order approximation the spinor wave function in equation (9) does not depend on the coordinates, C (0) ↑,↓ = const, an estimate for the lowest-order corrections to it being δ (C ↑,↓ ) ∼ αk/λ 2 . Then, the kinetic energy term in Hamiltonian (5) yields an obvious lowest-order contribution to the eigenvalue, that is k 2 /2. In the same approximation, a contribution of terms ∼ iα∂ x,y in (6) can be found as eigenvalues of the respective matrix with entries, cos θ(x, y) and sin θ(x, y), replaced by their spatially average values, cos θ(x, y) ≡ C and sin θ(x, y) ≡ S, and −i∂ x,y replaced by k x,y , the result being ±α √ C 2 + S 2 k. As concerns corrections to the eigenvalues produced by corrections to the eigenfunctions given by this δ (C ↑,↓ ), they scale as δµ ∼ (αk) 2 , being negligible in comparison with the terms written above. In fact, the current approximation is similar to the lowest-order perturbation theory in quantum mechanics, which produces the first contribution to the eigenvalue of energy, induced by a perturbation Hamiltonian, omitting a correction to the stationary wave function.
Thus, collecting the above terms, the lowest-order approximation for the spectrum at small k and α is obtained as This result is written with respect to µ lin (0) because the calculation of µ lin (0) is more cumbersome, as the contribution to it from the last two terms in Hamiltonian (4) requires a calculation of the first correction to the spinor eigenfunction, the result being For the SOC modulation in equation (3) with λ = π, a straightforward calculation using the well-known decomposition exp(iz cos ζ) = J 0 (z) + ∞ l=1 i k J l (z) cos(lζ) [50] with Bessel functions J l (z) yields C = J 2 0 (π/2) ≈ 0.223 and S = 0. Therefore, the analytically predicted spectrum (12) takes the following form: In the same approximation, the structure of the constant spinor eigenstate in equation (9) is determined by the ratio of its components (cf. equation (11)): The branch of dispersion relation (13) with the bottom sign has a minimum, µ min = α 2 /2 at a finite Bloch momentum, The derivative dk min /dα = J 2 0 (π/2) exactly coincides with its numerically found counterpart at α → 0, as shown in figure 2(c). The same figure suggests that the analytical approximation remains accurate enough up to α 0.5, which is corroborated by a detailed comparison with numerical data (not shown here).
Finally, the exact numerical calculation of µ at the bottom of the Bloch band, demonstrates that it is equal α 2 /2 for all α. Therefore, the relevant for soliton existence band edge µ be = α 2 /2. The black curve is the numerical result, while the short red line shows the analytical prediction for α λ, as given by equation (15). Note that, at large α, the radius k min increases as α, as expected for the uniform SOC.

Symmetry properties of solitons
Now we turn to exact soliton solutions of equation (1) in the form of Ψ = e −iµt ψ(x, y), with t-independent complex spinor ψ = (ψ 1 (x, y), ψ 2 (x, y)) T . Thus, we consider the nonlinear eigenvalue problem: where the linear part H of the model's Hamiltonian (1) is given by equation (4). The consideration below is limited to the case of a square lattice cell with equal x− and y−periods: λ x = λ y = λ, and hence the general symmetry properties pertain to this choice. To describe families of nonlinear solutions (which in our case do not bifurcate from the linear limit), we first address the symmetries of H. To this end we introduce time reversal operator T acting on ψ as the complex conjugation, T ψ = ψ * , x− and y− spatial-inversion operators P x and P y , acting as as well as P = P x P y mapping, (x, y) → (−x, −y), the exchange operator, and two operators of translations by half-period along the x and y directions: Now we can identify the symmetries of H.
Using properties P x JP x = P y JP y = PJ = JP, one can verify that β is a generator of a cyclic group of order 8 (β 8 = 1), which is isomorphic to Z 8 [51,52]. Both PT and β symmetries do not contain shifts, hence they are referred to as local symmetries. Meantime we observe that there exists also the discrete translational symmetry over the primitive lattice vectors, i.e., u 1 = (2, 0) and u 2 = (0, 2) in our case, as well as the symmetries related to the shifts over half-period and defined as σ 3 P x τ x τ y and σ 1 P y τ x τ y . Since here we are interested only in localized solutions, these lattice translations will not be considered. Taking into account that β 4 = −1, one identifies as the group of local symmetries of H. We emphasize that although the group G loc is obtained for the particular choice of θ(x, y), the following analysis is not restricted to this choice and depends only on the group G loc , which can be realized for a broad variety of the H−Hamiltonians. In G loc we use notation T F = β 2 PT = iσ 2 T for the spin-1/2 time reversal, and refer to elements with "+" ("−") sign as symmetries (anti-symmetries). Then, one can distinguish three cases as follows. If, for a given solution ψ of equation (16), there exists an element g of G loc such that ψ g = gψ is a solution of equation (16) with the same µ, we define ψ g as a g-transformed solution. If, additionally, ψ g = ψ (or ψ = gψ), such a solution is termed a g-symmetric one. If a given solution features ψ = gψ for all g ∈ G loc , then one says that symmetry breaking occurs. Although we cannot exclude a possibility of the symmetry breaking in our case, the extensive numerical studies of localized states of equation (16) have not revealed such a possibility. Therefore, we focus below on the existence and properties of "symmetric" solutions, i.e., ones generated by all symmetry transformations from G loc . Generally speaking, a symmetry of a linear Hamiltonian does not yet imply the symmetry of the complete nonlinear Hamiltonian. For a nonlinear solution to be generated by a symmetry g, the nonlinearity must be compatible with that symmetry. There exist two possibilities for that to occur [53]. Either the nonlinearity is gsymmetric, i.e., in our case, [g, ψ † ψ] = 0 for any ψ, or it is weakly g-symmetric, i.e., ψ † ψ = (gψ) † gψ is verified for a g-transformed solution gψ (although not necessarily for any ψ). It is straightforward to verify that the Manakov's (SU(2)-symmetric) nonlinearity in equation (16) is weakly symmetric with respect to local elements of G loc . Therefore, below we look for nonlinear modes generated by the symmetry group G loc .
Not all symmetries from G loc result in physically distinct solutions. For any symmetry involving antilinear time reversal T , the anti-symmetry amounts to a trivial phase shift. For instance, if solution ψ PT is a PT -symmetric one, then iψ PT is anti-PT -symmetric. Therefore, in the analysis of symmetric solutions of equation (16) presented below, we exclude all anti-symmetries. This leaves us with only four transformations: PT , T F , βPT , and β 3 PT , that may generate physically different solutions. Among these elements only the PT -symmetry allows for symmetric solutions: ψ PT = PT ψ PT . Now, let us consider ψ PT and its β 2 -transformation: ψ T F = β 2 ψ PT = T F ψ PT = iσ 2 Pψ PT . One readily concludes that the absolute value of the first (second) component of ψ T F is the P-inverse absolute value of the second (first) component of ψ PT . In this sense, these solutions are topologically similar. The topological similarity is also verified for the pair of βPT -and β 3 PT -symmetric solutions related by the β 2 transformation. On the other hand, PT -and βPT -symmetric solutions have essentially different distributions of the fields in both components.
The above analysis is applicable for any family of localized nonlinear modes. Although the detailed analysis of numerically obtained families (fundamental, dipole-, and quadrupole-like) will be done in the next Section, here we can illustrate their main symmetry-related properties. These solitons, found by the Newton relaxation method, that permits obtaining entire families parameterized by the chemical potential µ, are characterized by their norm U (µ) and size w(µ), defined as: and the peak amplitudes of spinor components a 1,2 (µ) = max |ψ 1,2 | . In figure 3 we show PT -symmetric solutions (panels (a) and (c)) and the respective β-transformed solutions (panels (b) and (d)) for the families of fundamental solitons (panels (a) and (b)) and dipoles (panels (c) and (d)). It is straightforward to verify that both U (µ) and w(µ) are equal for a given PT -symmetric soliton and its β ntransformation. It is also clear that the linear stability of a PT -symmetric state coincides with the stability of its β n transformation: if ψ PT is stable (unstable), the same is valid for β n ψ PT . Therefore, the nonlinear modes shown in figures 3 (a) and (b) (and in figures 3 (c) and (d)) are either both linearly stable or both unstable, in the former case being stabilized by the nonlinearity.

Soliton families and their stability: numerical results
The simplest family of self-trapped solutions supported by the 2D helicoidal SOC landscape is represented by fundamental solitons. To be stable, such solitons have to be centered at a maximum of the periodic SOC phase, θ(x, y) = π, with the chemical potential taking values below the edge of the first band in the linear spectrum, i.e., at µ ≤ µ be [see the dashed line in figure 4(a)]. The dominant component (either ψ 1 or ψ 2 , see Sec. 3) of such a soliton typically has a bell-shaped profile with weak modulations on the top of it, while a weaker component may feature a very complex profile with curved nodal lines, see, e.g., figure 4(c). Both components of the wave function have nontrivial phase distributions, as shown in the insets of figure 4(c). When the chemical potential approaches the bottom of the first linear band, the soliton expands and its width diverges, see a representative w(µ) dependence in figure 4(a). In this regime, spatial modulations of the soliton shape become much more pronounced for both components [ figure 4(d)]. Peak amplitudes a 1,2 (µ) of both components decrease with the increase of µ and vanish exactly at the edge of the allowed band [ figure 4(a)]. In contrast, the soliton's norm, U (µ), shows a non-monotonous behavior, clearly indicating that solitons exist only above a certain threshold value of the norm. This threshold decreases with the increase of the SOC strength α. At α → 0, the dependence U (µ) approaches a constant value, U T ≈ 5.85, which is the above-mentioned Townes' soliton norm [19,20,32].
The growth of deformation and non-monotonous character of the U (µ) dependence appearing with the increase of α suggests that the helicoidal SOC may play a stabilizing role for fundamental solitons. We have checked the stability of such states by simulating equation (1) with small input complex noise added to stationary profiles, up to very large times. As a result, it was found that large parts of U (µ) families, with negative slope ∂U (µ)/∂µ, represent stable fundamental solitons. These findings are summarized in the existence-and-stability diagram displayed in figure 4(b), where the solitons exist at µ ≤ µ be , and are stable for all values of µ ≤ µ cr , the border of the stability domain exactly corresponding to the ∂U (µ)/∂µ| µ=µcr = 0 boundary, in agreement with the well-known Vakhitov-Kolokolov criterion [19,20,54].
The existence of stable fundamental solitons residing at the maxima of the periodic θ(x, y) function suggests that equation (1) allows to construct complexes of the solitons, which cannot be done with the help of spatially-uniform SOC [29,33,34]. We have found that it is indeed possible to create stable complexes of spatially separated solitons with opposite signs of ψ, i.e., dipoles and quadrupoles. Examples of the simplest dipole soliton complex, that may be considered as a pair of two solitons with opposite signs, and properties of dipole families are presented in figure 5. The dipole is built of two soliton-like density maxima placed at adjacent maxima of θ(x, y) and, accordingly, separated by distance ≈ 2π/λ = 2. The solitons interact through their tails, the interaction getting stronger with the increase of chemical potential µ, due to weakening localization of each soliton [cf. figures 5(c) and (d)]. Even though phase distributions in the dipole modes are quite complex, due to spatially dependent SOC, one can see in the insets that the two density maxima in both components have, "on average", phase shifts of π, which corresponds to the opposite signs of the two solitons. Because dipole solitons are excited nonlinear states, their chemical potential does not reach the bottom of the band. Instead, the dipoles exist below a certain cutoff value, µ = µ co , whose α−dependence is shown in figure 5(b) by the blue line. At the cutoff position, amplitudes of both components of the dipole's wave function, ψ 1,2 , do not vanish, and the width remains finite, but dependence U (µ) acquires a vertical tangential direction [for illustrative purposes, in figure 5(a) we plot half of total norm, as it is comparable with that of fundamental soliton]. Because the capability of SOC landscape to compensate repulsive forces acting between solitons with opposite signs reduces with the decrease of α, the cutoff value µ co rapidly decreases at α → 0 [ figure 5(b)].
An important finding is that the dipole solitons enabled by helicoidal SOC may be stable. The stabilization is possible below a critical value of the chemical potential, µ ≤ µ cr [see the stable (black) segment of the U (µ) branch in figure 5(a)], while, close to the cutoff, the dipoles are unstable [see the red segment of the U (µ) branch]. The entire existence and stability domains for the dipole solitons is presented in figure 5(b): they exist at µ ≤ µ co , and are stable at µ ≤ µ cr . In the instability domain, the most typical destabilization mechanism is spontaneous breakup of symmetry between two density maxima, accompanied by oscillations with a growing amplitude, leading to the collapse or formation of a fundamental soliton. Finally, we have found quadrupole soliton complexes, composed of four solitons with alternating signs of ψ, residing at neighboring maxima of θ(x, y). Examples of quadrupoles and properties of their families are presented in figure 6. Note that in such states diagonal pairs of solitons have identical signs. In the course of the evolution, such pairs of mutually attractive solitons may fuse, destroying the whole structure. For this reason, such states are more fragile than their dipole counterparts, and the stability region for them is narrower, with a more complex shape. Quadrupoles also exist below the respective cutoff value of the chemical potential, µ ≤ µ co , see figure 6(a,b).
Due to the complexity of the stability domain, we were able to identify it only for a limited set of values of the SOC strengths, α. For instance, in figure 6(b) the quadrupoles are stable between two upper red dots, and also below the lower one. Stable segments of the U (µ) dependence in figure 6(a) are black ones. This stability was found solely in a limited range, 3 α 4, but not for larger SOC strengths.
Finally, in addition to the dipoles and quadrupoles, it is possible to construct multipole states arranged into lines, which are more stable than the quadrupoles. Figure 6. (a) The norm U (µ), size w(µ), and component amplitudes, a 1,2 (µ), of quadrupole solitons versus chemical potential µ at α = 3.5. Black and red segments in the U (µ) dependence correspond to stable and unstable states, respectively. (b) Existence and stability domains of the quadrupole solitons in the (α, µ) plane. They exist at µ ≤ µ co < µ be . Since the shape of the stability domain is very complex, in panel (b) we only show borders of the stability domain at α = 3.5. Panels (c) and (d) display examples of profiles of the absolute value and phase (insets) of quadrupoles with µ = 3 (c, stable) and µ = 5.6 (d, unstable) for α = 3.5. For convenience for quadrupole solitons the SOC profile was shifted by half of the period π/λ along both x-and y-axes.

Conclusions
In this work, we have introduced a 2D system which includes the cubic self-attraction and SOC (spin-orbit coupling) subject to the spatially-periodic modulation, as a physically relevant dynamical model of BEC in the binary ultracold bosonic gas. The system gives rise to two-component (pseudo-spinor) fundamental solitons in a vast region of the corresponding parameter space. It also essentially expands the variety of stable self-trapped modes in the 2D setting by the creation of families of robust dipole and quadrupole states built of two or four fundamental-like solitons with alternating signs of the wavefunction. Stability of these solutions and properties of their families, such as the norm, widths, and amplitudes, have been analyzed, parameterizing the families by the chemical potential. We have also shown that, for each species of the localized state (fundamental, dipole, or quadrupole), additional solutions can be generated from a given one by symmetry operations forming an Abelian eight-element group.
These results may be extended by considering axially symmetric (radial) patterns of the spatial SOC modulation, instead of the periodic lattice patterns. In the axisymmetric system, it may be interesting, in particular, to consider azimuthal mobility of solitons and interactions between them. A challenging possibility is to develop a threedimensional generalization of the system with the spatially periodic SOC modulation.