Vortex Solutions in a Binary Immiscible Bose-Einstein Condensate

We consider the mean-field vortex solutions and their stability within a two-component Bose Einstein condensate in the immiscible limit. A variational approach is employed to study a system consisting of a majority component which contains a single quantised vortex and a minority component which fills the vortex core. We show that a super-Gaussian function is a good approximation to the two-component vortex solution for a range of atom numbers of the in-filling component, by comparing the variational solutions to the full numerical solutions of the coupled Gross-Pitaevskii equations. We subsequently examine the stability of the vortex solutions by perturbing the in-filling component away from the centre of the vortex core, thereby demonstrating their stability to small perturbations.


I. INTRODUCTION
Quantized vortices have been realised in a variety of superfluid systems, including superfluid 4 He, phases of superfluid 3 He, and atomic Bose Einstein condensates (BECs) [1]; the latter was first achieved by using a phase imprinting method on a two-component BEC [2].Since this initial observation in BECs, vortices have been realised in a number of scenarios including rotating single [3,4] and multiple [5] component BECs, nucleation from a repulsive potential [6][7][8], and driving the condensate out of equilibrium in harmonic [9,10] and uniform [11,12] potentials.Of particular interest are the latter experiments [9][10][11][12], where the turbulent system manifests itself as many vortices distributed in a complex and disordered tangle.
In superfluid helium, where the typical size of a vortex core is of the order of 0.1nm [13], a range of methods have been theoretically and experimentally evaluated as ways to visualise the flow and track vortices.Many of these techniques involve the use of tracer particles, micron-sized solid particles which are suspended in the fluid and follow the local fluid velocity [13][14][15][16][17][18].A significant drawback of these techniques is that the particles are much larger than the vortices which they track, prohibiting detailed probing of quantum turbulence [19].By comparison, the most common way to image quantized vortices in a BEC is to allow the condensate to expand until the size of the vortex cores exceed the optical resolution limit and perform column-integrated imaging of the cloud [4,20].Advances in this technique now allow for imaging in multiple dimensions [21], and in real time [22,23].However, it remains an ongoing challenge to image a complex 3D distribution of many tangled vortex lines.
In this paper, we concentrate on a two-component BEC which is in the immiscible regime and where one component contains the majority of the total atoms.A vortex in the majority component is then known to become in-filled by the second component [2,47].Such a regime is a potential candidate for three-dimensional vortex detection in systems comprising complex vortex tangles, as the in-filling component can be tracked without destructively imaging the majority component.It is possible, however, that the presence of the in-filling component will modify the vortex states and their dynamics, and that this modification will primarily depend on the number of in-filling atoms and the inter-species interactions.It is important for the purpose of vortex detection to understand the regimes in which the in-filling component has no significant effect on the vortex dynamics, that is, is an effective passive tracer.It is also interesting to consider how, in more extreme cases, the in-filling might drive new regimes of vortex states and behaviour which have no analogy in conventional singlecomponent superfluids.In this paper we thus investigate the in-filled vortex solutions in an immiscible binary BEC, establishing their profiles through a variational approach and full numerical approaches, as well as their stability under perturbation.
The remainder of this paper is structured as follows: in Section II we introduce the equations of motion of the system, the coupled Gross-Pitaevskii Equation.In Section III we derive a variational approach to the in-filled vortex solutions, which allow the solutions to be established analytically.The variational solutions are compared to numerically-obtained solutions of the coupled Gross-Pitaevskii equation, showing excellent agreement.Section IV establishes the stability of the infilled vortices to small perturbations, before we present concluding remarks in Section V.

II. THE COUPLED GROSS-PITAEVSKII EQUATION
A system which consists of two BECs which are dilute and weakly-interacting in the zero-temperature limit is well described by the coupled Gross-Pitaevskii Equation (CGPE) where In what follows, we will consider species 1 to be the majority component, and refer to species 2 as the in-filling component.For simplicity, we set the external potentials to zero, such that the ground-state of the majority component is a state of uniform density.Then, it is useful to work in natural units of the majority component: density is in terms of the uniform density n 0,1 , length is in terms of the healing length, ξ 1 = / √ n 0,1 m 1 g 11 , energy is given by the chemical potential, µ 1 = n 0,1 g 11 , and time is given by τ 1 = /µ 1 .The CGPE, Eqns.(1a) and (1b), may then be cast in dimensionless form where the dimensionless parameters are the ratio of the atomic masses, m = m 1 /m 2 , the ratio of the inter-and intra-species interaction parameters, g 12 = g 12 /g 11 , and the ratio of the inter-species interaction parameters, g 22 = g 22 /g 11 .In the following, we will consider a two-component system which is in the immiscible limit, corresponding to g 12 2 > g 22 in dimensionless variables.For convenience, we will drop the primes in the remainder of this paper.

A. Overview
We consider an in-filled vortex which is aligned along the z-axis of the system, and aim to determine its cross-sectional profile.Since the vortex is assumed to be uniform along z, this reduces to a two-dimensional problem.The core of a vortex is characterised by a region of depleted density containing a point of zero density, around which the phase winds by an integer multiple of 2π [68].In the immiscible regime, it is energetically favourable for the minority component to be located at regions where the density of the majority component is lower [41].Therefore, if the majority component contains a vortex, we expect that the minority component "infills" the vortex core.
Throughout the remainder of the paper we will set g 22 = 1.0 and g 12 = 1.1 unless otherwise stated, meaning that the two species are just in the immiscible regime.We also take the ratio of masses to be equal, corresponding to two components which are made up of identical atomic species in different hyperfine states.
First we outline our approach to obtaining the full numerical solutions before establishing the semi-analytic variational approach.

B. Numerical solutions
We obtain the full vortex solutions by numerically solving the CGPE; these solutions will allow us to later test the success of the variational approach.
We solve the CGPE in two-dimensions using an adaptive RK45 method, with an error tolerance of 10 −10 , implemented using XMDS2 [69].We do so in a square computational grid with 2 numerical grid points per healing length, up to 1000 time-steps.In order to obtain the vortex solutions we impose an azimuthal phase on to the majority component, and evolve the CGPE, Eqns.(2a) and (2b), under a Wick rotation t → iτ.Under this so-called imagingary time propagation [70], the family of GPEs are well established to evolve towards the lowest energy state of the system.To avoid issues with the phase at the edge of the (periodic) domain, we impose a circular hard-wall potential with radius R and amplitude 100µ in each component to force its density to zero; we have checked that the radius R is sufficiently large that this imposed boundary has no effect on the in-filled vortex solutions.
Example numerical solutions from the CGPE are shown in Fig. 1.For a small minority component (panels (a, b)), the vortex density profile in the majority component resembles that of a vortex in a single-component BEC [70]; the density is zero at the centre of the vortex and relaxes to the background density value over a lengthscale characterised by the healing length ξ.The minority component is localised within the vortex core as a narrow wavepacket whose width is consistent with the vortex core, i.e. the healing length ξ.For a large minority component (panels (c,d)), however, the vortex profile in the majority component is much broader and flat-bottomed, with the minority component forming a broad flat-topped profile of similar width.
Physically, it is evident that, due to the immisiciblity of the two components, it is energetically favourable for the minority component to sit at the vortex core so as to minimise overlap of the two components.As more atoms are added to the minority component, these cause the vortex core to broaden, again so as to avoid overlap of the two components.

C. Variational solutions using a super-Gaussian ansatz
It is our aim here to establish a semi-analytic approach to the in-filled vortex solutions, for both components, using a variational method.Variational methods have been employed to find stationary solutions of a variety of systems, including vortex cores in a single-component BEC [68,71], bright solitons in BECs of attractive [72][73][74][75][76][77] and dipolar [78] atomic species, quantum droplets in vanilla [79,80] and dipolar BECs [81], and bosonic quantum impurities [82].The main advantage of these variational methods is the relative ease with which stationary solutions may be obtained, by comparison with the computational requirements of finding solutions to the full CGPE, and the results can often give useful physical insight into the properties of the solutions [70].
It can be shown that the energy functional corresponding to Eqns.(2a) and (2b) may be written as where In this section we proceed by substituting an ansatz solution into the energy functional, Eqn.(4); this results in an energy function which depends on parameters controlling the density profile of the two components.We minimise this function to find the variational solution.
The choice of ansatz for each component must satisfy the following limits.Firstly |ψ 1 (r)| 2 → n 0,1 and |ψ 2 (r)| 2 → 0 as r → ∞, which is to say that, well away from the vortex core, the density of the majority component relaxes to the density of the uniform background and the density of the minority component reduces to zero.Secondly, |ψ 1 (r)| 2 → 0 and 2,max as r → 0, where ψ 2 2,max is the peak density of the minority component.Moreover, the ansatz should be able to capture the range of profiles illustrated in Fig. 1 -from narrow, high-curvature profiles to broad, flattened profiles.Here we choose to base our ansatz on the super-Gaussian functions, which is sufficiently versatile to satisfy these criteria.
We write the ansatz solution for the in-filled vortex in polar coordinates (r, ϕ) as and where q is the charge of the vortex, λ is a parameter which characterizes the width of the vortex core and in-filling component.The exponential terms are the super-Gaussian function, and therein the exponent α controls the shape of the function.The pre-factors and normalise the ansatz wavefunctions to N 1 and N 2 respectively, where Γ (s, x) is the incomplete Gamma function [83], defined in Eqn.(A15).According to the super-Gaussian function, for α < 1 we obtain a cusp profile, while for α = 1 we recover a vanilla Gaussian curve, and for α > 1 the curve is a flattopped.Example fits of this ansatz to numerical solutions of the CGPE can be found in Fig. 1; clearly we see that the ansatz suitably captures both the narrow, high-curvature and the broad, flattened profiles presented.
We proceed analytically by substituting Eqns.( 6) and ( 7) into the energy functional, Eqn.(4), and integrating out the spatial dependence.This is a non-trivial calculation and further details can be found in the Appendix.To prevent the integrals involving the majority component from diverging, we consider the energy of the atoms within a finite distance R of the vortex core, and without loss of generality we place the vortex core at the origin [68].The resulting equation for the energy functional is where S = (R/λ) 2α , and we have introduced the Spence function, Sp, defined in Eqn.(A7), the Exponential integral, Ei, defined in Eqn.(A9), and the Euler-Mascheroni constant, γ EM ≈ 0.5772.
While approximations may exist, such that we can minimise Eqn.(10) analytically, we choose to numerically minimise the energy function as it appears in Eqn.(10).In order to do this, we apply the quasi-Newton method of Broyden, Fletcher, Goldfarb and Shanno [84], available in the scipy Python library.Thus, for given atom numbers N 1 and N 2 , the variational solution for the in-filled vortex is specified by the two parameters (α * , λ * ) which minimises the variational energy.These solutions agree well with the full numerical solutions, as evident in Fig. 1.
The variation of the variational solution parameters (α * , λ * ) with the number of atoms in the in-filling component, N 2 is shown in Fig. 2. For low N 2 , the values of λ * and α are small, giving rise to narrow, high-curvature profiles such as in Fig. 1(a,b).As N 2 is increased, the values of λ * and α * grow, indicating the broadening and flattening of the profiles in both components, such as the profiles in Fig. 1(c.d).

D. Accuracy of the variational solution
We now perform a quantitative assessment of the accuracy of the varational solution for the in-filled vortex compared to the full numerical solution of the CGPE.First we consider the energy of the solution.We compute the normalised error in the energy, where E(λ * , α * ) is the energy of the variational solution (the energy functional Eqn.(10) evaluated at (λ * , α * )) and E CGPE is the energy of the numerical solution.We also compute the normalised maximum deviation in the density of the variational solution from the numerical solution, where ψ k,Ansatz is the k-th component of the variational wavefunction, and ψ k,CGPE is the k-the component of the numerical wavefunction.In order that this statistic is not affected by the implementation of the circular hard-wall potential in the numerical solution, we compute this in the region r < R/2, well away from the hard walls.To judge the goodness of the shape of the variational wavefunction, we also calculate the normalised error in the Full Width at Half Maximum (FWHM) of the density profile, which is given by (13) The results of these metrics is plotted in Fig. 3.We see that the energy of the variational solution is accurate to within 5% of the numerical solutions throughout the full range of N 2 considered.In panel (c) we note that the normalised maximum error in the density of the in-filling component is largest for small N 2 .We suggest that this is due to the fact that the maximum value of the in-filling component is relatively small here [see panel (e)], which causes the normalised error to grow quickly.We observe that, while the variational solution for the in-filling component under-estimates the peak density, this error is mainly symptomatic of small N 2 , and for larger N 2 the maximum value of the in-filling density is in good agreement.
Of particular note is the close agreement between the variational solution and numerical solution for the majority component.This can be seen both in the deviation of the density profiles [panel (b)], and in the deviation of the FWHM [panel (d), blue pluses].The main motivation of this work is to establish how a second component might affect the ground state of a majority component which contains a vortex.With this in mind, we might regard the excellent agreement of the variational solution and the numerical solution in the majority component as being more important than the good agreement of the in-filling component.

A. Overview
Until now we have concentrated on calculating the stationary state of the in-filled vortex.We now turn to considering how stable the vortex solution is to perturbation.Specifically we will consider the response to perturbing the in-filling component and how localised the minority component remains within the vortex core; this is particularly relevant when considering the possibility to use the minority component as a tracer of vortices.

B. Perturbing the in-filling component
We prepare the perturbed state by forming the in-filled vortex solution as previously but then instantaneous translate the in-filling component by a distance x 0 along the x-axis, relative to the vortex core.We then evolve this system using Eqns.(1a) and (1b) in real time.We consider in-filling components with two different atom numbers, N 2 = 10 and N 2 = 100.The number of atoms in the majority component, N 1 is again chosen so that, far away from the trapping potential or the vortex core, the density is unitary.For the systems which we consider, the number of atoms in the majority component is approximately 5 × 10 3 and 5 × 10 2 larger than the number of We track the location of the in-filling component through two approaches.Firstly, we consider the centre of mass of the component; we define this as Secondly, we track the peak of the in-filling condensate's density, r 2 = max |ψ 2 | 2 .Example trajectories are plotted in Fig. 4 for an initial perturbation of 4ξ, and for three different parameter sets -(a) N 2 = 10 and g 12 = 1.1,(b) N 2 = 100 and g 12 = 1.1, and (c) N 2 = 10 and g 12 = 2.2.In all three of the cases considered, we observe that the in-filling component is stable to small perturbations away from the vortex core.One would expect this since, as the two components are immiscible, it is energetically favourable for the peak density of the minority component to be attracted to the minimum density in the majority component, i.e. the vortex core.We see that this is true over a wide range of in-filling atom numbers, as well as a range of inter-species interaction strengths.In all cases the in-filling component undergoes an irregular trajectory in the x − y plane; this is due to the non-trivial potential it experiences from the majority component.We see that, in the weakly immiscible system (g 12 = 1.1) shown in rows (a) and (b), the centre of mass of the system with a larger number of in-filling atoms is subjected to a more tightly confined trajectory.This is because the overlap interaction term (g 12 |ψ 1 | 2 |ψ 2 | 2 ) in the energy functional, Eqn. ( 4), grows faster with larger N 2 .Similarly, we find that for systems with a comparable number of atoms [N 2 = 10 in rows (a) and (c)], the centre of mass of the system with the stronger interspecies interaction strength undergoes a more constrained orbit.In rows (b) and (c) we observe that the trajectory of the peak density fluctuates more from the vortex core than the trajectory of the centre of mass.We suggest that, by comparison with row (a), the effective trap which the in-filling component experience (from the interaction potential with the majority component) has a larger radius and shallower gradient in the centre; thus there is more "sloshing" of the in-filling component, leading to greater variance in the position of the density peak.We observed that the vortex undergoes a negligible translation from the origin after perturbing the in-filling component.The fact that this translation is very small is due to the large imbalance between the number of atoms in each component.
We trace the coarse-grained density of the in-filling component in Fig. 5.The effect of instantaneously perturbing the in-filling component away from the vortex core is to generate sound waves, which propagate from the edge of the vortex core.Since the wave-front is an area where the density of the majority component is depleted (although it is non-zero, unlike the vortex core) it carries a small amount of the infilling component away from the centre of the vortex core.This wave-front collides with the hard-wall trapping potential, and is reflected back into the centre of the trap, interfering with itself.A result of this is that the density waves away from the vortex core have velocities which are radially both inward and outward.Over time, this leads to a redistribution of the in-filling component: while a large amount of the in-filling component remains within the vortex core throughout the simulation, the small amount which is displaced approaches a radial distribution which is approximately uniform [see Fig. 5 (g)-(i)].It is clear that, at late times, the majority of the in-filling component remains strongly localised within the vortex under perturbation, and that it's possible suitably traces the position of the vortex core.[85] Despite the fact that the in-filling component is only perturbed in the x direction, we also observe motion in the y direction.This is due to the fact that, in the majority component, the vortex imposes a velocity field about the origin.We may consider this superfluid velocity field as the velocity at which a particle would be advected [68], and hence, due to the small overlap between the two components, the in-filling component is subjected to a velocity field in the x and y component, as well as the oscillations which are due to perturbing the component away from the vortex core.This effect, combined with the non-trivial shape of the interaction potential experienced by the in-filling component from the density depletion in the majority component, lead to the in-filling component tracing out an irregular trajectory in the x − y plane.It is clear that the coupled solution is stable against small perturbations of the infilling component, and that the in-filling component remains localised within the vortex core.

V. CONCLUSIONS
We have considered a two-component Bose Einstein Condensate which is in the immiscible regieme, where one component (the majority component) contains a vortex and the other component (the minority component resides in the vortex core.For low in-filling atom numbers, the vortex profile is not significantly different from that of a single-component vortex, while a larger number of in-filling atoms leads to a broadening and flattening of the vortex core.We have pre-sented an ansatz for the wavefunction of each component in a uniform system based on a super-Gaussian function.Following a variational approach using this ansatz, we were able to shown that the parameters which minimise the GPE energy functional lead to wavefunctions which are in excellent agreement with the numerical solutions obtained by evolving the full coupled GPE equations for a range of atom numbers.This approach may be extended in the future to include trapping potentials or vortex pair solutions. We then proceeded to consider the response of the coupled vortex solution to perturbation.We were able to ascertain that the solution is stable against perturbations of the in-filling component away from the vortex core, for a range of atom numbers and inter-species interaction strengths.This work was partly motivated by the prospect of using the minority, in-filling component as a passive tracer of vortex lines in atomic BECs.Our work confirms two essential criteria for such a prospect -firstly, that the in-filling component remains localised in the vortex core (even under perturbation) and secondly, for suitably small atom numbers, has no significant affect on the vortex profile or back-action on the majority component.Further work is needed to establish how the in-filling component behaves in more complex vortex configurations, such as three-dimensional vortex tangles.This work was also partly motivated by whether the infilling component can alter the vortex properties and potentially open up new physical regimes of vortex dynamics.Indeed, the significant change to the vortex core profile for large numbers of in-filling atoms suggests a significant affect on the vortex-vortex interaction.This interaction underpins many macroscopic vortex phenomena such as quantum turbulence, Abrikosov vortex lattices, Onsager vortex states and the Berezinskii-Kosterlitz-Thouless transition.Studying how the in-filling component modifies the vortex-vortex interaction is an avenue for further work. where We give brief details on computing these in the following subsections.

Kinetic Terms
In order to compute the kinetic terms, Eqns.(A2a)-(A2c), we must find the gradients of ψ 1 and ψ 2 which are given by The second term is the kinetic energy of the azimuthally circulating fluid motion [68], which here is given by The Maclaurin series expansion of the integrand confirms that this integral does not diverge at the origin.Then it can be shown that where γ EM is the Euler-Mascheroni constant, γ EM ≈ 0.5772.Then the lower limit of the integral in Eqn.(A8) is , which tends to γ EM /2α as ε → 0, since the final term is a sum of positive powers of ε.Hence, the resulting form of the kinetic energy of the azimuthal motion is This is readily computed on making the substitution u(r) = (r/λ) 2α , resulting in

Interaction Terms
The final terms to compute are due to the interaction terms of the energy functional, Eqns.(A2d)-(A2f).It's possible, however, to save some work in noticing that these integrals contain terms of three forms:

FIG. 1 .
FIG. 1. Example density profiles of the in-filled vortex solutions for (a),(b) N 2 = 2 and (c),(d) N 2 = 1000.In panels (a) and (c) the blue solid curves show the numerically obtained density profile of the majority component, |ψ 1 | 2 , while the red dashed curves show the density profile given by the variational solution.In panels (b) and (d) the blue solid curve shows the numerically obtained density profile of the minority component, |ψ 2 | 2 , while the red dashed curve show the corresponding variational solution.In each case, g 12 = 1.1, g 22 = 1.0 and m = 1.0.

FIG. 2 .
FIG. 2. Parameters which minimise the energy of the ansatz, Eqn.(10), as the number of atoms in the in-filling component, N 2 , varies.Panel (a) contains the energy-minimizing width, λ * , while panel (b) contains the energy-minimizing exponent of the super-Gaussian function, α * .These solutions have scaled intra-species interaction g 22 = 1.0, inter-species interaction g 12 = 1.1, and mass ratio m = 1.0.

FIG. 3 .
FIG.3.Comparison of the numerically obtained ground state, and the ansatz wavefunction predicted by minimising the energy function, Eqn.(10).Panel (a), we plot the normalised error between the energy of the ansatz wavefunctions and the energy of the numerically obtained ground state, Eqn.(11).We plot the normalised maximum deviation of the ansatz wavefunction from the ground state wavefunction, Eqn.(12), for the density of the majority [panel (b)] and the in-fill [panel (c)] components.In panel (d), we plot the normalised error in the FWHM, Eqn.(13), for the majority, k = 1, component (blue pluses), and the in-fill, k = 2, component (red circles).In panel (e), we plot the maximum values of the in-filling density for both the ansatz (orange squares) and the numerical groundstate (purple triangles).

FIG. 4 .
FIG.4.Example trajectories of the in-filling component which is perturbed from the centre of a vortex core in the majority component.Column (i) shows the trajectory of the centre of mass of the infilling component, while column (ii) shows the trajectory of the peak of the in-filling component's density.In row (a), we consider a system where g 12 = 1.1 and N 2 = 10; in row (b) we consider a system where g 12 = 1.1 and N 2 = 100; in row (c) we consider a system where g 12 = 2.2 and N 2 = 10.In all cases, g 22 = 1.0 and m = 1.0.The location and size of the vortex in the majority component is indicated by the contour lines; the lines are separated by one-tenth of the background density of the majority component.In each case, the centre of the in-filling component is perturbed by 4ξ.

FIG. 5 .
FIG. 5. Normalized histogram showing the radial distribution of the infilling component density, |ψ 2 | 2 , after an initial perturbation of 4ξ from the vortex core, for R = 64ξ.Panels represent different snapshots in time: (a) t = 0, (b) t = 10, (c) t = 60, and (d) t = 120.Green crosses in panel (d) show the time-averaged distribution, averaged over the latter half of the simulation, 500 ≤ t ≤ 1000.Insets show the density profile of ψ 1 (blue), and ψ 2 (red), at corresponding times.

1 −α 2 2 [ 1
, Eqn. (A2a), is then given byE Kin,1a = πA 2 α 2 e −u du,(A5)where the second line is obtained by making the substitution u(r) = (r/λ) 2α .The solution isE Kin,1a = πA 2 + u(R)] exp [−u(R)] + Sp e −u(R) − 1 ,(A6) where we have introduced the Spence function[83] to consider the analytic continuation of Ei(x) along the negative real axis[83], given by E 1 (x) = −Ei(−x) which leads to the series expansion Ei(−x) = γ EM + log x + A11)The final kinetic term is due to the infill component, and is given byE Kin,2 = πmB 2 α 2 t) is the mean-field wavefunction of the kth component, k = 1, 2, and m k is the mass of the atomic species in the k-th component.Each of the components are independently subject to an external trapping potential V k (r).The intra-species interactions of the first and second components are parameterised by g 11 and g 22 respectively, while the inter-species interactions are given by g 12 .The density profile of each component is n 1 (r, t) = |ψ 1 (r, t)| 2 and n 2 (r, t) = |ψ 2 (r, t)| 2 , and each component is normalised to N 1 and N 2 atoms.