Vortex Synchronization in Bose-Einstein Condensates: A Time-Dependent Gross-Pitaevskii Equation Approach

In this work we consider vortex lattices in rotating Bose-Einstein Condensates composed of two species of bosons having different masses. Previously [1] it was claimed that the vortices of the two species form bound pairs and the two vortex lattices lock. Remarkably, the two condensates and the external drive all rotate at different speeds due to the disparity of the masses of the constituent bosons. In this paper we study the system by solving the full two-component Gross-Pitaevskii equations numerically. Using this approach we verify the stability of the putative locked state which is found to exist within a disk centered on the axis of rotation and which depends on the mass ratio of the two bosons. We also show that an analytic estimate of this locking radius based on a two-body force calculation agrees well with the numerical results.


I. INTRODUCTION
One of the most striking manifestations of the quantum-mechanical nature of superfluids under rotation is the formation of vortices [2,3]. Perhaps the most natural arena to controllably study the physics of vortices are Bose-Einstein Condensates (BECs) of alkali atoms [4,5,6]. For the simplest case where the condensate is composed of a single type of atom without spin degrees of freedom, a triangular lattice is formed [7,8]. On the other hand, for multicomponent systems (composed of mixtures of atoms or spinor condensates, for instance), the order parameter has additional degrees of freedom resulting in more complex vortex lattice structures (see, for instance, [9,10,11,12,13]).
For the classic problem of an ideal fluid in a container rotating at rate Ω, the steady-state local velocity is v = Ω × r where r is the distance from axis of rotation. Since this velocity has everywhere a nonvanishing curl (|∇ × v| = 2Ω), it is not a permissible flow for a superfluid which is supposed to be inherently irrotational (v = ∇θ with θ the phase of the SF order parameter). However, superfluids are well-known to mimic the classical rigid-body rotation on average by forming a vortex lattice where the density of these vortices is given by the Feynman relation [3,8] ρ v = mΩ πh (1) where m is the mass of the constituent bosons and Ω is the rate the superfluid is rotating which is equal to the rotational rate of the walls of the container. In a previous work [1], it was considered how the situation described above generalizes to the problem of twocomponent BECs composed of atoms having different masses. More specifically, Eq. (1) naturally generalizes for two-component systems to where m 1 and m 2 are the masses of the bosons in the two constituent condensates and Ω 1 and Ω 2 are angular rates at which the two superfluids are rotating. For the case where there is a negative interspecies scattering length, the attraction between species will lead to an attractive interaction between vortices of the two species. When this interaction is sufficiently large one has the situation where the vortices form bound pairs, forcing the densities of the two vortex lattices to be the same: For this case, Eqns. (2) imply that the two superfluids (taking without loss of generality m 1 > m 2 ) will rotate at different speeds Ω 1 < Ω 2 . This counterintuitive state results from the quantum mechanical nature of the superfluid and has no analog in the classical fluid case.
In [1] the existence of this state was argued by making an ansatz for the short-ranged interspecies interaction and performing a two-body force calculation using it. This gave a quantitative prediction for the distance from the center of the condensate at which the vortex pairs become unbound, resulting from the growth of the Magnus force, which is referred to as the locking radius. The goal of the current paper is to test these arguments by numerical integration of the full two-component Gross-Pitaevskii equation. We will verify that such locked states are stable for a range of parameters. Furthermore, we will see that the analytic prediction for the locking radius agrees well with the numerical results. This paper is organized as follows. First in Sec. II A we provide definitions and set the notation for the treatment of the Gross Pitaevskii equation. In Sec. II B we summarize the derivation of the locking radius previously given in [1] which is based on a two-body force calculation. Then in Sec. III we describe the split-operator technique utilized to propagate the Gross Pitaevskii equations in imaginary time. The main results of the paper are presented in Sec. IV. Here we provide the vortex lattice structures determined numerically, and compare them with the estimate for the locking radius. Finally, in Sec. V we provide a discussion of potential experiments to realize this effect and then conclude.

A. Two-component Gross-Pitaevskii Equations
Our analysis starts with the two-component Gross Pitaevskii energy functional in the frame of reference rotating at angular rate Ω d . This is given by E = and E 12 = g 12 d 2 r n 1 (r)n 2 (r).
In these equations, V 1 and V 2 are the confining potential of the BECs which we will take to be harmonic. The intraspecies and interspecies scattering strengths are defined as g 1,2 and g 12 respectively. The angular momentum operator, as usual, is defined as Varying this energy with respect to ψ 1 and ψ 2 and introducing the chemical potentials µ 1 and µ 2 to enforce particle number conservation gives the two-component Gross-Pitaevskii equations In Sec. III, we will describe how these coupled equations are solved numerically to find minima of the energy E.

B. Inter-species vortex attraction and locking
In this section, for completeness, we provide an estimate of the locking radius of the vortex-bound state based on a two-vortex calculation. This calculation is presented in more detail in [1]. We first consider the state for the case where the interspecies interaction is large and all of the vortices of one species are bound with that of the other due to the strong short-ranged attractive force. A bound pair of vortices is depicted in Fig. 1. The vortex binding causes the two superfluids to rotate at different rates (because of the different masses and the O X FIG. 1: A bound pair of vortices occurring in the locked state composed of a mixture of BECs with m1 > m2. The vortex in the heavier species is denoted with an 'x' while that in the lighter species is denoted with an 'o'. The center of the condensate is taken to be to the left of this bound pair. The attractive short-ranged interspecies force F 1,2 rstr serves to bind the vortex pairs together. This is counterbalanced by the Magnus force F 1,2 mag which increases from the center of the condensate. Note that the healing length of the superfluid is larger than the sphere representing vortices in this figure Feynman relation), creating a Magnus force which tries to rip the bound pair apart. The Magnus force is balanced by the short-ranged interspecies vortex interaction resulting from the overlap of the vortex cores. However, since the Magnus force grows linearly with the distance from the center of the condensate, it will eventually overcome the short-ranged interspecies attraction. The point at which this occurs we refer to as the locking radius.
We will now put the previous arguments on more quantitative footing. A vortex sitting at rest in a superfluid flowing at velocity v will experience a force perpendicular to the flow the so-called Magnus force [8], whereκ is a unit vector centered on the vortex pointing out the plane. Assuming the system is composed of entirely locked vortex lattices, we have for the vortex densities ρ v . This, via the Feynman relations, gives where Ω 1 and Ω 2 are the angular rotational rates of the two superfluids. Since m 1 > m 2 the superfluids will rotate at different rates which will lead to the Magnus forces pulling the bound pair apart. The short-ranged interspecies interaction arising from E 12 defined in Eq. (5) will counteract the Magnus force. We refer to this as the "restoring force". In order to obtain an analytic expression for this interaction, we take a Gaussian ansatz for the density profile about a vortex. Specifically, for a vortex in species α centered at r 0 , we take where λ α is on the order of the superfluid coherence length [1]. We take for the value of this length parameter. This is a slight modification of the analysis in [1] where the simpler case of λ α = ξ α was taken. We choose the the nonuniversal constant in Eq. (11) so that the ansatz in Eq. (10) provides a better fit to the density surrounding a vortex. For a more detailed discussion of this, see the Appendix . Inserting this ansatz for two vortices separated by distance d into Eq. (5) we find where we have dropped terms which do not depend on the vortex separation. The interspecies force immediately follows from the derivative of this interaction energy and is Balancing the forces on each vortex in the frame of reference rotating at the drive frequency, we have F 1 mag = F 1 rstr and F 2 mag = F 2 rstr , as is illustrated in Fig. 1. Since the restoring force acting on either species has the same magnitude we have that F 1 mag = F 2 mag . The Magnus force on species α in the frame rotating with the vortex lattice at frequency Ω v is given by which grows linearly with the distance from the center of the condensate r. Also, note that the restoring force will not depend on the position in the condensate. A bound vortex pair will become unstable when the Magnus force is equal to the maximum possible value of the restoring force. This can be worked out to be r c = 1 2e A more detailed derivation of this expression can be found in Ref. [1].

III. NUMERICAL METHODS
Since there are only rare occasions when the timedependent Gross-Pitaevskii equation permits analytic solutions, numerical simulation is often the method of choice for theoretically studying Bose-Einstein condensates (for a recent account numerical solution of the GPE, see Ref. [14]). In this section, for simplicity, we will only consider the single-component case, noting that the generalization to the two-component case is straightforward. To this end, the equation we wish to solve is which describes the evolution of ψ in imaginary time, τ = it. Under long enough evolution ψ will relax to the ground state of the Gross-Pitaevskii energy functional. In the above equation H is given by We use a split-operator method to evolve the order parameter ψ as in Eq. (16). The idea behind the splitoperator method is to approximate the evolution operator through imaginary time interval ∆τ , U (∆τ ) = e −H∆τ , by a product of terms which are easily diagonalizable. Neglecting for the moment the rotational term in Eq. (17), H can be written as the sum of two terms, These terms are easily diagonalized in momentum and position space respectively. The wave function ψ can then be advanced in time by ∆τ by which is accurate to second order in ∆τ . The order parameter can then be evolved by taking successive Fourier (and inverse Fourier) transforms of ψ and multiplying by the factors e − 1 2 T ∆τ , e −V ∆τ , and e −T 1 2 ∆τ respectively. Such Fourier transforms account for the bulk of the computational cost in this algorithm, thus using the efficient Fast Fourier Transform algorithm is crucial.
A complication in the above occurs due to the nonlinearity of the GPE. That is, V in our above prescription for time evolution depends on the density n = |ψ| 2 , and it is at first unclear for what time this quantity should be evaluated. It is shown in [15] that, provided we use the most updated version of the time-dependent density n = |ψ| 2 , Eq. (18) will retain its second order accuracy. The final complication occurs from the rotational term in H which is We have neglected this term thus far since it is diagonalized in neither position nor momentum space and therefore cannot be included in either T or V . However, we note that R commutes with both T and V so we can write ≈ e − 1 2 T ∆τ e −R∆τ e −V ∆τ e − 1 2 T ∆τ ψ(τ ).
Then we can perform a similar split-operator decomposition of the additional term as Evolution of ψ by this factor can then be performed by taking the partial Fourier transform of ψ, that is transforming over the x variables but leaving the y variables unchanged (or vice-versa). This completes the overview of the numerical method used to solve the GPE. As stated before, the generalization to the two-component case is straightforward.

IV. RESULTS
Next, we discuss the results of the numerical simulation. We first provide the parameters which were used for the computations. For simplicity, we restrict our attention to the simplest case where g 1 = g 2 , and we fix the interspecies scattering strength such that |g 12 |/g 1 = 2/3. Furthermore, we take the number of particles in each species to be the same: N 1 = N 2 . We set the dimensionless parameter defined asg ≡ m h 2 g 1 N 1 to beg = 2 × 10 4 which is in line with values from typical experiments [21]. We take the two trapping potentials to be harmonic and adjust their curvatures ω 1 , ω 2 so that the density profiles of the two species have the same Thomas-Fermi profiles. Finally we rotate the system at 0.9 times the critical rate at which the condensate becomes unstable due to centrifugal forces. We discretize the system on a 200 × 200 grid, and propagate the system in imaginary time intervals of ∆τ = 0.01 1 hωx . We first consider the simplest case where the masses of the two species are the same. For this state we take the initial wavefunction to be a perfect triangular lattice of vortices with density given by the Feynman relation, Eq. (1). This structure is then relaxed by evolving the wavefuntions in imaginary time using the methods described in Sec. III. As expected these vortex lattices remain fully locked. The relaxed structures show small deviations from the perfect triangular initial structure due to the effects of the trap [16]. The density profiles of these are shown Fig. 2 in panels (1a) and (1b). The positions of the vortices are determined by analyzing the phases of the relaxed wave functions. Using this relaxed structure as the initial state, we change the mass ratios and propagate the wavefunctions in imaginary time until convergence. Specifically, we consider the ratios of m 1 /m 2 = 1.2, 1.4. and 1.6 as shown in Fig. 2.
To compare these numerical results to our estimate for the locking radius described in Sec. II B, we need to tailor Eq. (15) to the case of a harmonic trap. We take the density profiles used in Eq. (15) to have the Thomas-Fermi form: where n 0 is the density at the center of the trap and R T F is the Thomas-Fermi radius (note that we are only considering the case when the two condensates have the same radius). Inserting this profile into Eq. (15) (taking the correct dependence of the coherence lengths on the density) one finds where r 0 c is Eq. (15) evaluated for parameters at the center of the trap. This equation can then be solved for r c to obtain the renormalized value of the locking radius This indicates that near the center of the condensate we will have r c ≈ r 0 c as expected. Also, the locking radius will never exceed the radius of the condensate as expected. The reduction of the bare value of the locking radius can be qualitatively understood as follows. The Magnus force is proportional to the superfluid density while the restoring force is proportional to this density squared. Therefore near the edge of the condensate where the density is considerably smaller than its value at the center, the Magnus force will be favored thereby suppressing the locking radius.
Shown in the third column of Fig. 2 are the positions of the vortices of the two species, labeled with x's and o's. The width of these labels are roughly the size of the coherence length of the condensates. Superimposed on this is the locking radius predicted by Eq. (24) (using Eq. (24) for the bar locking radius) shown as a dotted line. This shows that the analytic results provide an excellent estimate of the locking radius. Note that due to the strong interactions between the two condensates, an unbound vortex in one species will create a local minimum in the other. Such features can be seen in columns c and d of Fig. 2. These local depletions should not be mistaken for vortices which are defined by the phase behavior of the wavefunctions.

V. EXPERIMENTAL CONSIDERATIONS AND CONCLUDING REMARKS
The main requirement for realizing the locked state is a BEC composed of a binary mixture of atoms having different masses and a negative scattering length. Such a transition could be tuned with an interspecies Feshbach resonance which have been found in Li-Na [17] and Rb-K [18] mixtures. Such mixtures have respective mass ratios of 3.3 and 2.2. Another promising experimental system are mixtures of two isotopes of a particular atom. For instance, the interspecies scattering lengths of different species of Yb have been analyzed in [19] and are often found to be negative. Since the mass ratios for different isotopes are closer to unity, having a strong attractive interaction (often requiring a Feshbach resonance) is unnecessary to reach the vortex locked state for this case.
We also note that these results are closely related to the experiment described in [20]. Here a single-component BEC is stirred by a rotating optical lattice which acts as vortex pinning sites. When the optical lattice is rotated at the speed for which the density of the pinning sites matches the density of the vortex lattice predicted by Eq. (1), a completely locked state is observed. Away from this resonance, a similar analysis to the above will predict a disk of bound vortices.
In conclusion, we report the confirmation of the putative vortex locked state proposed in [1]. For this state, the two superfluids and the stirring potential all rotate at different rates, exhibiting an unusual effect due to the quantum mechanical nature of superfluids. In this paper, we showed that such a state exists within a disk centered on the axis of rotation and whose size agrees well with an analytic estimate. Note that our numerical analysis did not assume anything about the vortexvortex attraction (unlike our theoretical analysis, which assumes Eq. (12), and evolves the Gross-Pitaevskii equa-tions directly). Our results (both analytical and numerical) rely on approaching this state from the fully locked state. Experimentally, this is probably most easily realized by controllably adjusting an interspecies Feshbach resonance. Alternatively, one can use an optical lattice to control the effective masses of the atoms by varying the lattice depth. PHY-0456720 and PHY-0803371, and The Research Corporation Cottrell Scholars program (GR).

APPENDIX: DENSITY PROFILE FOR A SINGLE VORTEX
In order to find the short-ranged interspecies vortex interaction, we need to know the behavior of the density of the condensate about a vortex. To this end we consider the Gross-Pitaevskii equation for a single component BEC having a vortex at the origin. That is, we write ψ = f e iθ and take θ = ϕ where ϕ is the azimuthal angle from polar coordinates. Substituting this into the GPE leads to the following equation dictating the density profile The density n = f 2 resulting from the numerical solution of this equation is shown in Fig. 3.
The numerical solution shows that the density behavior close to the vortex core (r ≪ ξ) is n(r) ≈ 0.340 n 0 r ξ 2 . On the other hand, the far distance behavior is found to be n(r) = n 0 1 − r ξ 2 . To make our work amenable to analytic treatment, we take the following ansatz for the vortex profile: where λ is a parameter on the order of the coherence length. Note that while this ansatz has the correct form close to the vortex core, the long distance behavior differs considerably. Fortunately our problem of vortex locking is dominated by the short-distance behavior, and we choose λ so that the two densities (numerical and ansatz) agree at r = ξ which requires λ = 1.781ξ, as shown in Fig. 3.
Our positive results, confirming the vortex-locked state, also confirm our intuition that the origin of the phenomena is in the short-ranged attraction between vortices. As explained in Ref. 1 (but without proof), the algebraic decay of the superfluid order parameter of a single votex does not imply that vortices of one species, when in a lattice, exhibit a power-law decaying force on the vortices on the other species. Unlike the single-species vortex-vortex force, which is the result of the inductive (kinetic) energy term in the Gross-Pitaevskii equations, the interspecies force is a result of a density-density interaction. The density suppression due to a single vortex occurs since the superflow of the vortex effectively increases the mass terms V 1 , V 2 in Eq. (4). But in a lattice of vortices, the combined superflow vector is nearly zero (i.e., negligible compared toh/m α ξ α ), and, therefore, so is the respective density suppression.