Introduction

Multidimensional localized modes in nonlinear dispersive media, which are usually considered as solitons, is a topic of great interest in many areas of physics1,2,3,4,5,6,7,8,9,10, especially in nonlinear photonics and in the studies of matter waves in Bose-Einstein condensates (BECs). In addition to their significance to fundamental studies, two- and three-dimensional (2D and 3D) solitons offer potential applications, such as the use of spatiotemporal optical solitons as bit carriers in data-processing schemes11 and the design of matter-wave interferometers using 3D matter-wave solitons12,13,14,15,16.

Unlike 1D solitons, which are usually stable objects17,18, stability is a major issue for their 2D and 3D counterparts. The ubiquitous cubic self-focusing nonlinearity, which readily produces formal 2D and 3D soliton solutions, gives rise to the wave collapse (critical collapse in 2D and supercritical collapse in 3D)19,20,21,22,23, which entirely destabilizes the respective soliton families. For this reason, the first example of solitons which was predicted in nonlinear optics, viz., the Townes solitons24, i.e., 2D self-trapped modes supported by a cubic self-focusing nonlinearity, are subject to the instability caused by the critical collapse, which is why they have not been created in experiment. Multidimensional solitons with intrinsic vorticity (vortex rings), which also formally exist in cubic media (in particular, vortex counterparts of the Townes solitons25), are subject to an even stronger splitting instability initiated by azimuthal perturbations1.

Thus, the stabilization of multidimensional fundamental and vortex solitons is an issue of great importance in various physical settings. In a 2D geometry with cubic self-focusing nonlinearity, the stabilizing factor must break the specific scale invariance, which underlies the critical collapse. The most straightforward possibility for that is provided by the addition of a spatially periodic (lattice) potential, which fixes a particular spatial scale (the lattice period). As a result, it was predicted that such potentials provide for the stabilization of both fundamental and vortex solitons26,27,28,29,30,31,32,33,34,35. In most cases, stable vortices are created in the form of a ring-shaped chain of local density peaks, with the vorticity carried by the superimposed phase profile which features a phase circulation of 2πS, with integer S being the respective topological charge.

In many cases, the natural interactions in the underlying physical media are repulsive, hence they cannot create regular solitons, which play the role of the ground state in the physical system, although it is possible to create gap solitons as a result of the interplay of the lattice potential and a repulsive nonlinearity36,37,38,39,40. In particular, inter-atomic forces in bosonic gases are normally repulsive41, as well as the exciton-exciton interactions in polariton condensates in semiconductor microcavities42,43,44,45.

The overall repulsive interaction may be switched into attraction by means of a Feshbach resonance (FR), in bosonic gases46,47,48,49 and polariton condensates50 alike. The FR applies to the inter-component interactions in binary (spinor or pseudo-spinor) condensates51,52,53. This suggests a new approach to the creation of solitons in spinor condensates: while each component features self-repulsion, the FR-induced attraction between them makes it possible to create symbiotic solitons, supported solely by the attraction between the two components, which may even overcome the intrinsic repulsion in each of them54,55. In particular, for microcavity polaritons this condition may be satisfied in a narrow spectral range, as explained in more detail below. The formation of localized symbiotic modes in the presence of a lattice potential was studied too, but only in 1D settings56,57.

The objective of the present work is to develop the concept of symbiotic solitons for 2D spinor condensates with an inter-component cubic attraction and intra-species repulsion. To secure the stability of the 2D solitons against the critical collapse, the system must include a lattice potential, which can be easily implemented in experiments for both atomic BEC (as an optical lattice)41 and polariton condensates58 (in the latter work, the lattice was used for the experimental creation of 2D exciton-polariton solitons).

The results presented below are obtained by means of a combination of an analytical variational approximation (VA) and a systematic numerical investigation of the existence, stability, and dynamics of the solitons. We find that the VA provides very good accuracy in predicting fundamental and vortex solitons, both symmetric and asymmetric ones, with respect to the two components. We also consider self-trapped modes with hidden vorticity (HV), i.e., a soliton with opposite signs of the vorticity in its two components59. The solution families include both gap solitons and regular ones, as concerns their placing relative to the bandgap spectrum of the linearized system. It is found that, in the presence of a lattice, all the fundamental solitons are completely stable solutions, while the vortices have a stability region when the cross-component attraction is stronger than the intra-species repulsion. The HV modes may be stable in the opposite case.

It is further worth noting that if the filling factor of the lattice is too small or the lattice potential is too weak, 2D solitons supported by the cubic nonlinearity with either attractive or repulsive sign suffer a delocalization transition, i.e., they cannot exist as stationary states trapped by the lattice60. In this work, we do not aim to exactly identify delocalization boundaries.

Results

The model

The dynamics of both two-component atomic BEC and coherent microcavity polariton gases (with losses appropriately compensated by gain), are modeled by the symmetric pair of coupled nonlinear Schrödinger equations41,42,

Here ϕ(x, y, t) and ψ(x, y, t) are the coherent fields of the two components, with ∇2 ≡ ∂2/∂x2 + ∂2/∂y2 being the 2D Laplacian and V(x, y) an effective lattice potential. It is assumed, as noted above, that the self-interaction of each component is repulsive, while the cross-interaction is attractive, accounted for by g > 0; another mechanism that may give rise to effective attraction in the polariton fluid was proposed61. The equations are written in the scaled form, so that the coefficients in front of the Laplacian and self-repulsion terms are normalized to unity. Equal coefficients in front of the Laplacians in both equations imply that the spinor describes a pair of different hyperfine states of the same atom in the bosonic gas. The same symmetry refers to the two spin states of excitons in the polariton gas62.

The lattice potential is taken in the usual form, V = −[cos(2πx/P) + cos(2πy/P)], with spatial period P. Below, generic numerical results are adequately presented for P = 10. The depth of the potential, Vmax − Vmin = 4, is fixed by means of the remaining scaling invariance. At the origin, x = y = 0, where the center of the soliton will be placed, the lattice potential has a local minimum. It is easy to check that the shape of the lattice different from quadratic, adopted here (e.g., hexagonal, radial, or anisotropic) will not essentially affect the results30.

In our analysis, we numerically sought for stationary localized solutions of Eqs (1) and (2), in the form of

with corresponding chemical potentials λ and μ and real-valued functions u and v obeying the stationary equations,

Numerical solutions were constructed on a finite-size grid that was sufficiently large to avoid influence of boundaries on the localized solutions. Stationary solutions were found starting from a localized input, using the imaginary-time propagation method63,64. The stability of these stationary solutions is then tested by simulations of perturbed evolution of Eqs (1) and (2) in real time.

Symmetric solitons and vortices

First, we look for symmetric modes with identical stationary wave functions u and v of the two components, with equal norms

In this symmetric case, the overall nonlinearity present in Eqs (1) and (2) can be characterized by a single effective coefficient g − 1. Thus, g > 1 makes the effective nonlinearity self-attractive. For g < 1 the strength of the repulsive inter-component interaction is larger than the attractive cross-component interaction and thus the net effect in this case is an overall defocusing nonlinearity. Note that the net nonlinearity vanishes in the symmetric case for g = 1.

Numerical results for the symmetric solutions are summarized in Fig. 1. We have found that the symmetric fundamental solitons [see a typical example in Fig. 1(a)] are stable in their entire existence region, as shown by Fig. 1(b). This conclusion complies with the fact that, for g > 1 and g < 1, the N(μ) dependences obey, severally, the Vakhitov-Kolokolov (VK) and anti-VK criteria, i.e., dN/dμ < 0 and dN/dμ > 0, which are necessary, although not sufficient, stability conditions for solitons supported, respectively, by attractive22,65 and repulsive66 nonlinearities.

Figure 1
figure 1

Symmetric solitons and vortices.

(a) The amplitude distribution in a stable symmetric fundamental soliton, for g = 1.1 and M = 20. (b) Dependence of the norm (M ≡ N) on the chemical potential (λ ≡ μ) for the family of fundamental symmetric solitons at different values of the cross-attraction coefficient, g. Solid lines depict stable numerically found solutions, while triangles represent the corresponding results produced by the VA as per Eqs (10)–(13). (c,d) Comparison of numerically found (solid lines) and VA-predicted (dashed lines) profiles of the fundamental symmetric solitons at g = 0.8, M = 40 (c) and g = 1.2, M = 40 (d). Panels (e) and (f) display, respectively the amplitude and phase distributions in a numerically found stable symmetric vortex solitons for g = 1.1 and M = 20. (g) Dependence of the norm (M ≡ N) on the chemical potential (λ ≡ μ) for the symmetric vortex solitons for different values of g. Solid and dashed lines represent stable and unstable vortices, respectively, while triangles depict the corresponding VA results, as per Eqs (16)–(19). At g > 1, the vortices are stable in a finite interval of values of the norm, see Eq. (6). (h) Dependence of the stability boundaries of vortices, Nmax and Nmin, on g. The gray regions in (b,g) represent the (narrow) first and second bands of the lattice-potential spectrum, respectively.

Full size image

In fact, the fundamental solitons found at g < 1, which correspond to the effective repulsive nonlinearity, are gap solitons36. It can be readily checked that values of λ ≡ μ in Fig. 1(b), corresponding to g < 1, fall into the first finite bandgap of the spectrum of the linearized symmetric Eq. (4), while those corresponding to g > 1 belong to the semi-infinite gap. These gaps are separated by the (first) narrow band, which is displayed in Fig. 1(b) by a gray-shaded stripe. Unlike the fundamental solitons, the vortices corresponding to g > 1 and g < 1 fall, respectively, into the first and second finite bandgaps, which are separated by the narrow second band, as shown in Fig. 1(g). It is usually assumed that gap solitons have a loosely bound shape, with many inner oscillations36. Nevertheless, they may also feature a tightly bound shape, similar to that of regular solitons, and in that case they can be effectively fitted by means of the VA39,67, similar to what occurs in the present system.

Also included in Fig. 1(b,c) are the analytical results produced by the VA for the fundamental solitons, as detailed in the Methods Section below. The VA-predicted results agree extremely well with their numerical counterparts (very small deviations are observed for stronger nonlinearity and/or larger norms).

While, as mentioned above, lattice potentials often create vortex states with a complex structure, composed of four density peaks for the sates with vorticity S = ±126, simple crater-shaped stable vortices may be found too68. Symmetric crater-shaped vortices have been found in the present case too, with a typical example displayed in Fig. 1(e,f). As Fig. 1(g) shows, the version of the VA developed for the vortex solitons (see Methods Section for details) also provides a very good accuracy, in comparison with the numerical findings.

As said above, all vortices (here, we consider only those with S = ±1) have been found in the first and second finite bandgaps (but not in the semi-infinite gap). The vortices are unstable for g < 1, when they belong to the second bandgap. On the other hand, Fig. 1(g) demonstrates that families of the symmetric vortex solitons are stable at g > 1 in the first gap, in a finite interval of

where Nmax decreases with the increase of g, see Fig. 1(h). Direct simulations demonstrate that unstable vortices spontaneously evolve into robust randomly vibrating single- or multi-peak patterns, which are asymmetric with respect to the two components.

Asymmetric solitons and vortices

Asymmetric solitons, with unequal u and v components, form an especially interesting class of modes in the present model. We start their consideration with the case of g = 1, in which symmetric solitons cannot exist. Basic results for this case are displayed in Fig. 2. Note that the width of the profile with the larger norm [Fig. 2(a)] is larger than that of its counterpart [Fig. 2(b)] with smaller norm. For a quantitative analysis we define the asymmetry ratio, R ≡ N/M. Then, with g = 1 and the asymmetry ratio 0 < R < 1, i.e., roughly speaking, |ψ|2 > |ϕ|2, Eqs (1) and (2) suggest that the defocusing and focusing nonlinearities dominate in the former and latter equations, respectively. This conclusion is also confirmed by the behavior of the chemical potentials λ and μ, cf. their behavior in Fig. 2 with that observed in Fig. 1(b) for g < 1 and g > 1, respectively. Thus, asymmetric solitons with g = 1 may be considered as bound states of gap solitons and regular ones, i.e., complexes of the semi-gap type56, the respective chemical potentials, λ and μ, falling, respectively, into the first finite gap and semi-infinite one.

Figure 2
figure 2

Asymmetric solitons.

Amplitude distributions of the u- (a) and v- (b) components of the stable asymmetric fundamental soliton with g = 1 and N = 10, M = 20 (i.e., the asymmetry ratio is R = 0.5). Panels (c) and (d) display dependences between the larger norm (M) and the two chemical potentials, λ and μ, respectively [see Eqs (5), (12) and (13)], for different values of the asymmetric ratio, R, at g = 1. Solid lines represent stable numerical solutions, and triangles depict the VA-predicted analytical results from Eqs (10)–(13).

Full size image

The analytical results the VA produces for these asymmetric modes are also found to be in excellent agreement with the numerical findings, as can be seen from the comparison of the analytical and numerical results for varying chemical potentials λ and μ in Fig. 2(c,d). A detailed comparison (not shown here in detail) demonstrates that, as the u component becomes broader with decrease of the asymmetry ratio R, this component is more strongly affected by the underlying lattice potential, which slightly worsens the VA accuracy for stronger nonlinearities in the limit of R → 0.

Similar to their symmetric counterparts, it has been found that the asymmetric fundamental solitons are stable in their entire existence region. Intuitively, this conclusion agrees with the fact that the λ(M) and μ(M) dependencies in Fig. 2(c,d) satisfy the anti-VK and VK criteria, respectively.

For asymmetric vortices, the spatial profiles of the u and v components feature small differences, as seen in Fig. 3(a,b), while the phase distributions are virtually identical in both components, provided that they carry the same vorticity, S = +1 or −1, in both components (not shown here). Similar to what is shown above for the fundamental solitons in the case of g = 1, when the solutions are asymmetric (M ≠ N), the repulsive and attractive interactions dominate in the different components, resulting in different dependence between the chemical potentials and the norm, as shown in Fig. 3(c).

Figure 3
figure 3

Asymmetric vortices.

Amplitude distributions of the u- (a) and v- (b) components of the asymmetric (R = 0.5) vortex soliton at g = 1 and M = 20. (c) Dependence of the norm (M) on the chemical potentials (μ and λ) of the two components of the vortex for different asymmetry ratios at g = 1. Solid and dashed lines denote stable and instable solutions, respectively, while triangles represent the VA predictions, as per Eqs (16)–(19). Regions of stability and instability for the vortices are shown for g = 1 in (d), g = 1.1 in (e), and g = 1.2 in (f).

Full size image

The comparison of the VA for asymmetric vortices and numerical results is shown in Fig. 3(c). In the asymmetric case, the vortex component (u) with a larger norm is broader than its counterpart with a smaller norm (v). Also in this case, due to the influence of the underlying lattice, the agreement of the exact numerical solution with the simplest vortex ansatz, given by Eq. (14), gets worse with increasing vortex size. Therefore, the agreement of the VA and numerical results for the v component is better than for the u component, see Fig. 3(c). Generally, the size of the vortices is larger than that of the fundamental solitons, therefore the overall agreement of the variational and numerical results for asymmetric vortices in Fig. 3(c) is somewhat poorer than for the asymmetric fundamental solitons in Fig. 2(c,d).

Figure 3(d) shows stability and instability regions for the vortices, depending on the asymmetry ratio, for g = 1. At R → 1, the vortices become unstable, as the nonlinearity effectively disappears when the solutions approach the symmetric case with g = 1. At R → 0, the model reduces to a single-component system with a defocusing nonlinearity. The smallest norm for stable vortices is found around R = 0.7.

Figure 3(e) shows an example of the stability region for vortices for g > 1, viz., with g = 1.1. Starting from the symmetric case with R = 1, the stability range initially slightly widens with the increase of the asymmetry (decrease of R). Then, the stability region gradually decreases with further decrease of R, until all stable solutions disappear at R < 0.6. For g = 1.2, as shown in Fig. 3(f), the stability region decreases monotonously with increasing asymmetry (decrease of R), again shrinking to nil at R ≈ 0.6. At g > 1, smaller asymmetry ratio R implies domination of the attractive nonlinearity in Eq. (2), eventually leading to an instability of the v component. It is worth noting that the fact that vortices are unstable for larger norms in Fig. 3(e,f) is in contrast to the results for g = 1, shown in Fig. 3(d). Roughly speaking, the stability regions for g > 1, shown in Fig. 3(e,f) are rotated by 90° in comparison with their counterpart shown in Fig. 3(d) for g = 1. At g < 1, all asymmetric vortices are unstable, similar to what is stated above for the symmetric ones.

Hidden-vorticity (HV) modes

For the HV vortex states, in which the two components have opposite vorticities, results are summarized in Fig. 4. The vorticities of the two components are S = +1 and −1, corresponding to the phase maps shown in Figs 4(b) and 1(f), respectively. In the symmetric case (M = N), both components have the same spatial profiles, as shown in Fig. 4(a). In this case, the mode carries zero angular momentum, therefore it is called the HV state59.

Figure 4
figure 4

Hidden vorticity modes.

(a) Amplitude and (b) phase distributions of symmetric (R = 1) stationary solutions for the u-component of an vortex-antivortex pair at g = 0.8 and M = 35. Dependences of the norm (M) on the chemical potential (λ) of vortex-antivortex pairs for (c) different nonlinearity at R = 1 and (d) different asymmetric ratio at g = 0.8. Solid lines are stable numerical solutions, while dashed lines are unstable solutions. Triangles depict analytical results from Eqs (16)–(19).

Full size image

We have found that the vortex-antivortex pairs are unstable for the overall-focusing nonlinear system, with g > 1, but a remarkable fact is that they may be stable for the defocusing system, with g < 1, as shown in Fig. 4(c). Recall that all states with explicit vorticity, i.e., S = 1 in both components, are entirely unstable at g < 1, as shown above. We have also found finite stability ranges for the vortex-antivortex pairs at g < 1 in the weakly asymmetric case, with asymmetry ratios R slightly smaller than 1, as shown in Fig. 4(d). If the norm difference between the two components becomes too large, the vortex-antivortex pairs become unstable.

The VA results for the vortex-antivortex pairs are given by the same equation, Eqs (18) and (19), as for their vortex-vortex counterparts. The analytical results agree very well with the numerical finding when M = N, see Fig. 4(c). For M ≠ N, the difference between the analytical and numerical results becomes more apparent due to the larger size of the vortex components, as seen in Fig. 4(d).

Discussion

We have numerically and theoretically investigated the existence and stability of symbiotic solitons and vortices in the 2D two-component spinor system with repulsive intra-species and attractive cross-component interactions, in the presence of an underlying lattice potential. We have found that the fundamental solitons are always stable, in both symmetric and asymmetric cases (equal or different norms of the two components). Vortex solitons with the same sign of the vorticity in both components are stable when the attractive cross-component interaction is stronger than the intra-species repulsion. On the other hand, the modes with hidden vorticity, i.e., opposite signs of the vorticity in the two components, can be stable only when the repulsive interactions dominate. These novel types of self-trapped modes are of interest for the understanding of general principles of pattern formation in 2D nonlinear settings, and they can be realized in physical systems. In atomic gases, the inter-particle interaction can be finely tuned in the vicinity of a FR (Feshbach resonance)46, using the strongly dispersive nature of the scattering length near the two-particle resonance (a negative scattering length is associated with an attractive interaction).

This FR control may also be realized for microcavity polaritons, where the frequency dependence of the scattering matrix element has been analyzed in detail for excitations which are placed spectrally below the fundamental exciton resonance69,70. In a narrow spectral range close to the two-particle resonance associated with the formation of a bound biexciton state, the role of the attractive cross-component interaction can be finely tuned so that it may dominate over the repulsive inter-component interaction. We note that in this spectral range, in addition to the intrinsic loss rate for polaritons, the system will also suffer excitation-induced losses (a similar side effect of the FR is known in atomic BECs). However, these losses may be compensated in the presence of an exciton reservoir, or by an external pump laser, which makes it possible to maintain the quasi-stationary behavior. In this case, the stationary modes predicted in the present work should be observable in semiconductor microcavities.

There remain a number of points for extension of the work, such as hybrid bound states, with a fundamental self-trapped state in one component, and a vortex in the other. A challenging problem is the consideration of a three-dimensional version of the present system, which may be realized in an atomic BEC.

Methods

The variational approximation (VA)

The Lagrangian of the stationary equations (4) is

For fundamental solitons, the simplest Gaussian ansatz is adopted, with amplitudes A, B and radial widths a, b:

This ansatz produces the following norms, Eq. (5), in the two components: M = πA2a2, and N = πB2b2. In the asymmetric case, M ≠ N, we define the asymmetry ratio R = N/M as noted above, restricting it to 0 < R ≤ 1. The substitution of ansatz (8) into Lagrangian (7) and subsequent integration yields the following effective Lagrangian:

The first pair of the variational equation following from Eq. (9) makes it possible to express a2 and b2 in terms of M and N: ∂Leff/∂(a2) = ∂Leff/∂(b2) = 0, i.e.,

For given M and N, a2 and b2 can be found by means of a numerical solutions of the algebraic system (10) and (11). The chemical potentials, λ and μ, do not appear in Eqs (10) and (11). They are produced by the second pair of the variational equations: ∂Leff/∂M = ∂Leff/∂N = 0, i.e.,

The first objective of the VA is, for given g, to identify a region in the plane of (M, N) where the numerical solution of Eqs (10) and (11) produces physically relevant solutions, with a2, b2 > 0.

Similarly to the approach for the fundamental solitons outlined above, for vortices we adopt a natural generalization of the Gaussian ansatz:

whose norms are M = πA2a4 and N = πB2b4. Substituting ansatz (14) into Lagrangian (7), the following effective Lagrangian is obtained, cf. Eq. (9):

Then, the variational equations for a2 and b2 are produced as ∂Leff/∂(a2) = ∂Leff/∂(b2) = 0, i.e.,

Finally, the chemical potentials for the two components of the vortex are produced according by the remaining variational equations, ∂Leff/∂M = ∂Leff/∂N = 0, i.e.,

Additional Information

How to cite this article: Ma, X. et al. Two-dimensional symbiotic solitons and vortices in binary condensates with attractive cross-species interaction. Sci. Rep. 6, 34847; doi: 10.1038/srep34847 (2016).