Stable VortexBright Soliton Structures in TwoComponent Bose Einstein Condensates
Abstract
We report the numerical realization of robust 2component structures in 2d and 3d BoseEinstein Condensates with nontrivial topological charge in one component. We identify a stable symbiotic state in which a higherdimensional bright soliton exists even in a homogeneous setting with defocusing interactions, due to the effective potential created by a stable vortex in the other component. The resulting vortexbright solitary waves, generalizations of the recently experimentally observed darkbright solitons, are found to be very robust in both in the homogeneous medium and in the presence of parabolic and periodic external confinement.
Introduction. Vortices in nonlinear field theory have a timehonored history (1). They are among the most striking features of superfluids, play a role in critical current densities and resistances of typeII superconductors through their transport properties, and are associated with quantum turbulence in superfluid helium (2). The advent of BoseEinstein condensates (BECs) 15 years ago (3); (4) has produced an ideal setting for exploring relevant phenomena. Since the experimental observation of matterwave vortices (5), by using a phaseimprinting method between two hyperfine spin states of a Rb BEC (6), the road opened for an extensive examination of vortex formation, dynamics and interactions. Stirring the BECs (7) above a certain critical angular speed (8); (9); (10); (11) led to the production of few vortices (11) and even of very robust vortex lattices (12). These structures have been produced by other experimental techniques, such as dragging obstacles through the BEC (13) or the nonlinear interference of condensate fragments (14). Later, not only unitcharged, but also highercharged structures were produced (15) and their dynamical (in)stability was examined. This field also has strong similarities and overlap with the emergence of vortices and even vortex lattices in nonlinear optical settings; see e.g. (16); (17).
Another remarkable possibility in both BEC (18); (19); (20) and in nonlinear optics, e.g. (21), is that of multicomponent settings. Matter waves exhibit rich phase separation dynamics driven by the nonlinear interatomic interactions between different species or states that make up the BECs. Longitudinal spin waves (22), transitions between triangular and interlaced square vortex lattices (23), striated magnetic domains (24); (25), and robust target patterns (26) have all been observed, as well as tunable interspecies interactions (27) and transitions between miscible and immiscible dynamics (28).
We interweave these two settings, motivated by (29) in which
darkbright solitons have been created in a
quasi1d
twocomponent
BECs. These structures were predicted (30) and extended to more
complex settings such as spinor condensates (31) (as darkdarkbright,
or darkbrightbright solutions), but were only realized experimentally in
2008. These are often termed “symbiotic solitons”, as the bright component
would be impossible to sustain under repulsive interatomic interactions
(i.e., defocusing nonlinearities, as considered here), unless the
darkcomponent creates an “effective potential”, of which the bright soliton
is a bound state.
The coupled bright solitary waves (32) and the
gap ones of (33) constitute additional examples of symbiotic
structures.
We consider higherdimensional realizations
(30)
i.e. vortexbright
solitons of various
topological charge in 2d, as well as in
3d
Physical Setup. The nondimensional Hamiltonian for a twocomponent condensate in the meanfield approximation reads (36):
(1) 
where is the pseudospinor order parameter, , is the diagonal matrix of chemical potentials associated with the conservation of the number of atoms and ; a related useful diagnostic is . is a matrix accounting for the effectively nonlinear interatomic interactions. For the and components of Rb we can use (26) , and . These determine, through the negative sign of det, the immiscible nature of the interactions leading to phase separation (19); (26). The confining potential is
(2) 
where is the parabolic component (often created magnetically) and the periodic (optical) lattice component. The time and length scales are and , where is the atomic mass and is an arbitrary frequency in Hz. For Rb with scattering length of nm, Hz, and choosing , the ratio between the actual and nondimensional number of atoms is , where is the dimensional interaction parameter. For a 2d reduction, the interaction parameter is (e.g. (38)) and taking , the amplification factor is . The equations of motion , where and interchanges rows with , for this infinitedimensional Hamiltonian system are
(3) 
The stability of stationary solutions is determined by the eigenvalues of the Hessian of the Hamiltonian, , and of . Negative eigenvalues of indicate energetic instability, since “dissipative” perturbations (e.g. from exchanges of atoms with the thermal cloud if the temperature deviates from zero) in the system can render them dynamically unstable, as can collisions with other eigendirections even in the pure Hamiltonian (zerotemperature) system. The linear stability of the latter system is examined through the eigenvalues of ; instability arises when since, due to the Hamiltonian structure, the eigenvalues are symmetric over both the real and imaginary axes. From prior experience which is confirmed again here, linear stability indicates evolutionary nonlinear stability in the meanfield model, at least for time scales on the order of tens of seconds (this is not generically the case for nonlinear static solutions of Hamiltonian systems).
The variation can be posed in the or the basis. The former is useful when the potential is axisymmetric, since then small excitations to a stationary solution of the form will have definite angular momentum . If we set then , , and , so a single index will indicate the angular momentum of the excitation with given eigenvalue . Hence, the spectrum of eigenvalues can be decomposed as the union of the spectra pertaining to angular momentum . We will also assume , so that and . It has been shown numerically (34) and analytically (37) that instability windows arise in a single component with topological charge only for wave numbers with . The null eigenvalues corresponding to gauge invariance appear in the spectrum of . For a single component in a parabolic trap, an anomalous mode for converges to zero as , accounting for translational invariance and leading to the energetic stability of the vortex without external potential. For each () an anomalous mode leads to windows of instability (34); (37). We show that these can be significantly suppressed, although the spectrum occasionally leads to small instability windows for a small fraction brightsoliton component, , for large with parabolic trap (no windows were observed without the trap).
Numerical Methods. Our methods extend those in, e.g., (39); (40). The spatial discretization in employs Chebyshev polynomials to represent dependence (41). The Fourier modes representing and make the Laplacian operators diagonal in these directions. To identify stationary states of (3), we first obtain an initial estimate via imaginarytime (i.e., replacing ) integration using a firstorder implicit/explicit Euler scheme with . We then refine the solution using Newton’s method. The linear system arising at each Newton step is solved using the matrixfree IDR(s) algorithm (42); (43), which requires only the action of the Hessian. To accelerate inversion, we precondition the system with the inverse Laplacian, using its block diagonal structure. Hence, we solve the system and update for . Fewer than 5 Newton iterations usually achieve an accuracy of .
For each stationary solution , we use the matrixfree Implicitly Restarted Arnoldi algorithm to iteratively compute the eigenpairs of the linearization to a specified tolerance (47). In order to find the desired eigenvalues we use inverse iteration, with the IDR(s) method and inverse Laplacian preconditioning to solve the linear systems, as above. Here, the preconditioner is taken to be , so that each iteration solves .
We used a resolution in of to represent nonaxisymmetric solutions and eigenvectors. For axisymmetric solutions, quantitative accuracy requires only 30 radial modes for , but up to 200 modes for larger . For eigenvectors, we use only modes in (see introduction) and identify quantitatively all expected invariant and negative directions and windows of instability from (34).
Results. In 2d () for , we first demonstrate in Fig. 1 the existence of an energetically stable (and hence also dynamically stable) vortexbright soliton state. This is so for all of the values that we have sampled. Notice the symbiotic nature of the state, as a bright soliton would be impossible to support under the repulsive/selfdefocusing interactions considered herein. Indeed, a similar state exists for vortices of higher topological charge , as shown in Fig. 1 for . The logarithmic scale shows that the soliton is more localized for larger . In this case, the negative energy modes (34) in the spectra associated to may lead to dynamical instability from complex quartets of eigenvalues as a result of HamiltonianHopf bifurcations which can occur upon collision with positive modes (44).
In a parabolic trap, the vortexbright structure remains dynamically stable. However, breaking of translational invariance produces a negative energy mode in the spectrum. Thus, an additional control parameter (through the atom number of the bright component) may lead to collisions of this mode with positive energy modes and, hence, rare isolated windows of instability arising for large and large ; see Fig. 2 (left) for such a window for as a function of when . A similar feature has recently been shown in the darkbright 1d analog of the vortexbright states (45).
For highercharge vortexbright structures, the same situation holds for the spectrum, while for windows of instability arise from the negative modes in the spectra of . As decreases, however, these windows of instability are generically suppressed by the increasing presence of the second bright component. Fig. 2 (right) depicts the growth rate of the spectrum for over parameter space. An example of an unstable solution perturbed in the growing excitation direction is depicted in Fig. 3. The vortices first split from the center, and begin to part and precess, but once they are far enough and the bright component is bimodal, they approach again and the bright component resumes unimodality. The sequence repeats, similarly to singlecomponent vortices (48).
When we impose an additional sinusoidal lattice potential, , the onecomponent () vortex may become unstable (due to resonant eigenvalue collisions and ensuing oscillatory instabilities), at least for sufficiently large (46). The same holds for a large mass ratio . However, below a critical , once again the bright component has a stabilizing influence.
The vortexbright structure is stable in 3d without the trap and with periodic boundary conditions in . Indeed this is immediately clear upon Fourier transforming in , since the spectrum of the Hessian decouples into an infinite family of subspectra equal to the 2d spectra shifted by , and hence it remains nonnegative. It is stable in the trapped case as well for and , and indeed for . When , the solution has another rotational invariance, and additional negative energy modes emerge for . For there are at least two additional negative energy modes, although this may not lead to dynamical instability. Upon addition of the lattice, ,the results are expected to be similar to 2d (up to considerations of the aspect ratio of the harmonic trapping). See Fig. 4 for an example with .
Discussion. We have generalized the darkbright, quasi1d soliton that has been predicted theoretically and observed experimentally in BECs to that of a vortexbright robust dynamical entity that emerges as a stable structure in both 2d and 3d condensates (although similar concepts could be directly applicable to the nonlinear optics of defocusing optical media). We also examined relevant structures in the presence of parabolic (magnetic) and periodic (optical) trapping and found that they remain stable. While instabilities may arise (e.g. for higher topological charge, , or as a result of the lattice), these are usually alleviated/suppressed by the presence of the second component in 2d and 3d.
It would be interesting to determine the robust existence of such waveforms, which are well within the reach of recent experiments, e.g. (28); (26). Our study suggests that highercharge vortices in a single component may be stabilized by an external bluedetuned laserbeam potential acting as the brightsoliton here. Hence, the stability of such vortices should be systematically examined in the presence of external potentials. Other themes such as multivortexbright soliton interactions and lattices would also be natural extensions of the present work.
Acknowledgments. PGK gratefully acknowledges support from NSFDMS0349023, 0806762 and the Alexander von Humboldt Foundation.
Footnotes
 We consider only repulsive interactions, where such structures are symbiotic. With an attractive interaction, where such states can be selftrapped with a vanishing tail, they were originally proposed in Z.H. Musslimani et al., Phys. Rev. Lett. 84, 1164 (2000).
References
 L.M. Pismen, Vortices in Nonlinear Fields, Oxford Science Publications (Oxford, 1999).
 R.J. Donnelly, Quantized Vortices in Helium II, Cambridge University Press (New York, 1991); D.R. Tilley and J. Tilley, Superfluidity and Superconductivity, IOP Publishing (Philadelphia, 1990).
 M.H. Anderson et al., Science 269, 198 (1995).
 K.B. Davis, et al., Phys. Rev. Lett., 75, 3969 (1995).
 M. R. Matthews et al., Phys. Rev. Lett. 83, 2498 (1999).
 J. E. Williams and M. J. Holland, Nature 401, 568 (1999).
 K. W. Madison et al., Phys. Rev. Lett. 84, 806 (2000).
 A. Recati, F. Zambelli, and S. Stringari, Phys. Rev. Lett. 86, 377 (2001).
 S. Sinha and Y. Castin, Phys. Rev. Lett. 87, 190402 (2001).
 I. Corro et al. PRA 80, 033609 (2009).
 K. W. Madison et al., Phys. Rev. Lett. 86, 4443 (2001).
 C. Raman et al., Phys. Rev. Lett. 87, 210402 (2001).
 R. Onofrio et al., Phys. Rev. Lett. 85, 2228 (2000).
 D. R. Scherer et al., Phys. Rev. Lett. 98, 110402 (2007).
 A.E. Leanhardt et al., Phys. Rev. Lett. 89, 190403 (2002); Y. Shin et al., Phys. Rev. Lett. 93, 160406 (2004).
 Yu. S. Kivshar et al., Opt. Commun. 152, 198 (1998).
 A. Dreischuh et al., J. Opt. Soc. Am. B 19, 550 (2002).
 C.J. Myatt et al., Phys. Rev. Lett. 78, 586 (1997).
 D. S. Hall et al., Phys. Rev. Lett. 81, 1539 (1998).
 D.M. StamperKurn et al., Phys. Rev. Lett. 80, 2027 (1998).
 D. Rand et al., Phys. Rev. Lett. 98, 053902 (2007).
 H.J. Lewandowski et al., Phys. Rev. Lett. 88, 070403 (2002).
 V. Schweikhard et al., Phys. Rev. Lett. 93, 210403 (2004).
 H.J. Miesner et al., Phys. Rev. Lett. 82, 2228 (1999).
 J. Stenger et al., Nature 396, 345 (1998).
 K.M. Mertes et al., Phys. Rev. Lett. 99, 190402 (2007).
 G. Thalhammer et al., Phys. Rev. Lett. 100, 210402 (2008)
 S.B. Papp, J.M. Pino and C.E. Wieman, Phys. Rev. Lett. 101, 040402 (2008).
 C. Becker et al., Nature Phys. 4, 496 (2008).
 Th. Busch and J.R. Anglin, Phys. Rev. Lett. 87, 010401 (2001).
 H.E. Nistazakis et al., Phys. Rev. A 77, 033612 (2008).
 V.M. PérezGarcía, J. Belmonte Beitia, Phys. Rev. A 72, 033620 (2005).
 S.K. Adhikari and B.A. Malomed, Phys. Rev. A 77, 023607 (2008).
 H. Pu et al., Phys. Rev. A 59, 1533 (1999).
 D. L. Feder, et al., Phys Rev, Lett. 86, 564 (2001).
 C.J. Pethick and H. Smith, BoseEinstein condensation in dilute gases, Cambridge University Press (Cambridge, 2002).

R. Kollár and R. L. Pego,
http://www.math.cmu.edu/
CNA/Publications/publications2009/017abs/017abs.html (preprint).  W. Bao, D. Jaksch, P. A. Markowich. J. Comp. Phys. 187 (2003), 318342.
 C. Huepe et al., Phys. Rev. A 68, 023609 (2003).
 T. Kapitula, K.J.H. Law, P.G. Kevrekidis, SIAM J. Appl. Dyn. Syst. 9, 34 (2010).
 L. N. Trefethen, Spectral Methods in MATLAB, SIAM (Philadelphia, 2000).
 Peter Sonneveld and Martin B. van Gijzen, SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 10351062 (2008).
 http://ta.twi.tudelft.nl/NW/users/gijzen/IDR.html
 J. C. van der Meer, Nonlinearity 3, 1041 (1990).
 S. Middelkamp et al., arXiv:1005.3789.
 K. J. H. Law, et al. Phys. Rev. A 77, 053612 (2008)
 http://www.caam.rice.edu/software/ARPACK/
 H.M. Nilsen and E. Lundh, Phys. Rev. A 77, 013604 (2008).