Synthetic superfluid chemistry with vortex-trapped quantum impurities

We explore the effect of using two-dimensional matter-wave vortices to confine an ensemble of bosonic quantum impurities. This is modelled theoretically using a mass-imbalanced homogeneous two component Gross-Pitaevskii equation. By changing the mass imbalance of our system we find the shape of the vortices are deformed even at modest imbalances, leading to `barrel' shaped vortices; which we quantify using a multi-component variational approach. The energy of impurity carrying vortex pairs are computed, revealing a mass-dependent energy splitting. We then compute the excited states of the impurity, which we in turn use to construct `covalent bonds' for vortex pairs, leading to the prediction of impurity mediated bound states. Our work opens a new route to simulating synthetic chemical reactions with superfluid systems.


I. INTRODUCTION
Quantized vortices represent the fundamental excitations of superfluids -atomic gases formed from interacting particles that exhibit non-viscous transport phenomena. Early experimental work revealed the nature of superfluid vortices [1][2][3] in these systems, which has also spurred theoretical interest, since the Gross-Pitaevskii formalism facilitates the accurate modelling of the phenomenology of superfluid alkali gases [4,5].
Atomic Bose-Einstein condensates represent exceptionally pure physical systems, since it is possible to experimentally realize ground states with almost no noncondensate atoms present. This gives a unique opportunity to study the role of impurities in these systems, with recent pioneering experiments realizing both bosonic [6,7] and fermionic [8] polarons. Complementary to this, experiments with binary condensates have now achieved the trapping of one matter-wave inside another, here a degenerate Fermi gas formed of 6 Li atoms was confined inside a Bose-Einstein condensate of 133 Cs atoms [9]. Attractive atomic interactions provide another opportunity to understand the interaction of impurities with solitary waves in the nonlinear regime [10][11][12].
Homogeneous systems represent an important testing ground for theoretical ideas, and have been under active experimental pursuit. Experiments with cold atomic gases typically use magnetic or optical traps to provide spatial confinement. However, recently it has become possible to confine atomic gases inside potentials that realize an effective box trap, leading to homogeneous bosonic ground states in one [13,14], two [15], and threedimensional systems [16,17]. Degeneracy has also been achieved with Fermi gases confined in homogeneous potentials [18]. Homogeneous systems offer a unique opportunity to investigate superfluidity in a clean and precise way [19] and understand the non-equilibrium physics of uniform systems [20].
Quantum mechanical gases manifest superfluidity by the nucleation of quantized vortices when the gas un-dergoes rotation. Typically this leads to the formation of the Abrikosov lattice at equilibrium, however recent work has revealed that homogeneous [21], multicomponent [22,23] and density-dependent [24] gauge theories all exhibit novel vortex configurations. Since individual vortices are topologically protected, they represent an intriguing candidate for hosting impurities within the matter-wave [25]. In the context of the strongly interacting Helium fluids, charged impurities in the form of electrons were originally proposed as an extension to the Gross-Pitaevskii formalism [26], which has recently found renewed interest due to the ability to accurately numerically simulate the Gross-Clark model [27].
The weakly correlated regime represented by atomic Bose-Einstein condensates provides an important testing ground for understanding the dynamics of superfluid vortices [28], which in turn provides new theoretical [29,30] and experimental [31] insight into quantum turbulence and non-equilibrium phenomenology [32]. Understanding the complex dynamics of superfluids is assisted by knowledge of the few-body dynamics of vortices, in particular the realization of vortex dipoles [33] continues to provide important information [34,35], as well as the observation of solitonic vortices [36] in elongated superfluids. Multiple vortex states are generated by rotating the atomic cloud [37][38][39], leading at large rotation to the formation of the Abrikosov vortex lattice [40,41].
Quantum mechanical gases with several coupled interacting internal degrees of freedom present an ongoing opportunity to gain a deeper insight into exotic forms of superfluidity in microscopic systems [42]. Condensates with coherent couplings have also been an ongoing focus, with the dynamics of vortices in these models showing physics analogous to high energy physics [43][44][45][46][47][48][49][50][51]. Condensates with spin degrees of freedom have been shown to possess rich vortex physics, including stable multiply charged vortices [52,53] in the ferromagnetic phase of these systems, attributed to the existence of both mass and spin current, leading to two individual forms of circulation. It has also been possible to consider "half" quantized vortices in the antiferromagnetic phase, where the   circulation is still quantized but the unit of quantization becomes a fraction leading again to unique lattice phases [54,55], as well as exotic vortex structures with nontrivial core structures which have been studied in spinor condensates [56][57][58][59][60][61]. In multi-component systems, vortices are usually core-less [62] -the cores of vortices in one component are filled by the other component(s). This leads non-trivial vortex-vortex interactions [63,64] and consequently rich phase structures of vortex lattices other than the Abrikosov vortex lattice [47,62,[65][66][67].
Multi-component systems can also host more elaborate topological structures such as monopoles [68,69] consisting of a radially varying spin structure, which have been realized in the ferromagnetic spin-1 phase of a 87 Rb condensate [70]. It is also possible to consider even more elaborate and exotic excitations, such as skyrmions in 2D [71] and 3D [72][73][74][75][76] which arise as a spatially varying spin deformation of the ground state of spinor systems as well as knots that are characterized in terms of the Hopf invariant, a topological charge unique to these excitations [77][78][79]. Multi-component systems also host exotic composite topological excitations such as D-brane solitonsvortices ending on a domain wall [80][81][82]. Thus, multicomponent superfluids provide an ongoing resource for topology, facilitating the creation of these exotic excitations in a stable and controlled environment.
We consider a two-component system in which only the first component is condensed, while the second feels the superfluid component as an effective potential. We study how impurity atoms can be trapped and manipulated by individual and pairs of vortices in homogeneous superfluids. By changing the number of atoms in the individual components, we reveal the unusual phenomenology of this system; including distorted vortex profiles and mass-dependent splitting of the impurities energy. We then explore the properties of the excited states of the impurity, including the possibility of synthetic chemical bonds.
The paper is organized as follows. In Sec. II we introduce the model that describes our binary system, before examining the physics of individual vortices coupled to differing numbers of impurities in Sec. III, including a variational description of the system to explain the effect of the impurities on the vortex. We then compute the energy of pairs of vortices carrying impurities with different combinations of charges. Section. IV examines the exited states of the vortex-impurity system, revealing the possibility of vortex chemistry with superfluids. In Sec. V we summarize our findings.

II. THEORETICAL MODEL
We consider a system of N 1 atoms forming a Bose-Einstein condensate described by the state ψ 1 , selfconsistently collisionally coupled to a second component of N 2 impurity atoms, described respectively by the state ψ 2 . Then, after projecting the three-dimensional binary system into two-dimensions by integrating out the axial degree of freedom one obtains the coupled Gross-Pitaeveskii equations where g 11 = 2 a s /ma 2 ρ encapsulates the scaled threedimensional s-wave scattering length a s of the atoms with mass m, while g 12 defines the strength of the collisional coupling between the superfluid and the impurity. The normalization of each component obeys dr|ψ j (r)| 2 = N j , while the mass imbalance is N 2 /N 1 . Meanwhile, the total (conserved) energy associated with the model of Eqs. (1) is defined as (2) A single isolated two-dimensional vortex constitutes a hole on an otherwise homogeneous background, so it is important to separate the energy associated with the background and that of the vortex itself. As such we can write the total density of the superfluid as [83] the value of the superfluid far from the vortex core) while the term in brackets accounts for the vortex core region. Then, the atom number for the superfluid component can be expressed as where A 2D defines the area of the first component. The atom number N 1 is comprised of a constant background atom number A 2D n (0) 1 and the atom number of the vortex core. To obtain the energy of a single vortex in the two-component system, the energy associated with the homogeneous (vortex free) state and by writing E 0 /n (0) 1 higher order terms have been dropped, assuming that the size of the system L x,y satisfies L x,y ξ, where ξ = / √ mµ 1 is the size (healing length) of the vortex core. Then, the energy associated with the vortex E v+i = E − E 0 can be found by our expression for E 0 using Eq. (2), giving where the individual contributions can be written as Equations (3a)-(3c) will be used in Sec. III to calculate the energy of different vortex configurations.

III. NUMERICAL RESULTS
To understand the physical behaviour of the coupled superfluid-impurity system, we perform numerical simulations of the two-component Gross-Pitaevskii equation, Eq. (1), subject to the constraint that the phase distri- bution ϑ 1 (x, y) of the first component satisfies which constitutes N v vortices of charge q j placed individually at (x j , y j ). In our work we will take q j = ±1.
Although it is also possible to consider vortices of higher winding number, it is well-known that such configurations are energetically unstable for the model of Eq. (1). We work in the healing units, where the healing length is defined as ξ = / √ mn 1 g 11 , the units of energy are  (1) subject to the phase constraint of Eq. (4) for a single vortex with q 1 = +1 (viz. Eq. (4)). Throughout we have taken N 1 = 4×10 3 . Panels (a) and (b) present cross-sections of the density of the vortex-carrying superfluid |ψ 1 | 2 and the impurity |ψ 2 | 2 as a function of the number of impurity atoms N 2 respectively. The effect of changing the number of impurity atoms is quite dramatic -the shape of the vortex core is observed to change from the characteristic funnel shape presented in panel (i) showing the density and phase for N 2 = 1 to a barrel-shaped hole in the fluid, as presented in panel (ii) for N 2 = 391. Then, individual density cross-sections are presented in panel (c) for the superfluid and impurity components. One can see the qualitative difference between the vortices' core shape (red dash-dotted) for N 1 = 4×10 3 and N 2 = 391 (blue dashed). The corresponding impurity wave functions are plotted in the corresponding colors and styles; while the individual impurity densities are displayed in panels (d-f) for N 2 = 1, 41, 391. The final panel Fig. 1 (g) shows a calculation of the width (standard deviation) of the vortex core (black solid) and the impurity wave function (grey dashed).

B. Variational calculation
To build intuition and explore the physics of the vortex impurity model, we introduce a variational approach to understand the static properties of the mass-imbalanced system. As such we introduce three variational parameters α, β and σ which can be physically interpreted as the vortex curvature, inverse vortex radius and impurity length scale, respectively. Then a general purpose variational ansatz (centered at the origin) describing the coupling between the two components ψ var 1,2 for a vortex of charge q with varying mass imbalance N 2 /N 1 will be taken as The ansatz given by Eqs. (5a) and (5b) will be used in what follows to compute the variational energy of the system, as well as to obtain an approximate analytical expression for the angular momentum of the vortex. In order to compare the values of the energy obtained numerically from the Gross-Pitaevskii equations (1) and the renormalized energy of Eqs. (3a)-(3c), we use the super- 1 − |ψ var 1 | 2 which consequentially allows us to perform a semi-analytic computation of the energy. Proceeding, we can substitute Eqs. (5a)-(5b) into Eqs. (3a)-(3c) to yield the individual contributions to the total variational energy E var v+i , yielding the semi-analytic expressions To evaluate the various integrals that contribute to the total energy (Eqs. (3a)-(3c)), we can exploit the cylindrical symmetry of the problem and evaluate the energy in polar coordinates. The integrals required to obtain these results are listed in appendix A. Note that the limits of integration (arising from the vortex phase) in Eq. (6a) re-quire a little care. The lower limit is the (scaled) length scale of the vortex, the healing length ξ v . This quantity can be approximated from the chemical potential of the homogeneous system, µ var 1 = g 11 n (0)  Eqs. (5a) and (5b). Then ξ v = / 2mµ var 1 . The upper limit (∝ R) of this integral is simply the 'radius' of this system, which can be obtained directly for a given simulation using R = A ⊥ /π where A ⊥ is the area of the two-dimensional numerical box. We wish to compare the variational solutions to the energy, Eqs. (6a)-(6d) to the exact values computed using Eqs. (3a)-(3c) and the numerical solutions to the Gross-Pitaevskii model. We use a least squares method to procure the variational parameters α, β and σ for a given number of atoms (N 1 , N 2 ) and interaction strengths g 11 , g 12 . As such we wish to calculate (1). Equation (7) can then be used to obtain a set of variational parameters for particular atom numbers and scattering lengths. Figure 2 presents the solutions for α/ξ, β and σ/ξ obtained from Eq. (7). The first two of these are the vortex variational parameters, plotted as α/ξ (blue crosses) and 1/β (green triangles). One can see that the vortex 'radius' 1/β is less for small concentrations of impurity atoms, and increases approximately linearly for larger numbers of impurities. The vortex curvature α on the other hand increases for small concentrations of impurities, before leveling off at larger impurity numbers. There appears to be some correlation between the vortex curvature α/ξ and radius 1/β, which seems to occur for N 2 100, i.e. at the point where the curvature α/ξ levels off and the radius 1/β starts to grow (grey dashed line). A possible explanation for this is that at this concentration of impurity atoms the shape of the vortex starts to significantly change, developing the flatbottomed profile (Fig. 1(c)), leading in turn to competition between the two terms appearing in the first fraction in square brackets in Eq. (5a). The impurity length scale σ/ξ is also shown (open red circles). An accompanying fit is displayed (black dashed) to σ fit /ξ ∝ √ N 2 , which can be inferred from Eq. (6b), showing good agreement.
The expressions for the variational energy, Eqs. (6a)-(6d) can be compared to the exact values found from the numerical solutions of the two-component Gross-Pitaevskii equation, Eq. (1) and the renormalized energy of Eqs. (3a)-(3c). This is presented in figure 3. Panel (a) shows a comparison of the total variational energy (green squares) and the total energy calculated from the exact numerical solutions (red crosses) for ψ 1,2 as a function of the number of impurity atoms, showing good agreement. Then, panels (b-e) show the individual contributions to the total energy, while the two kinetic terms for the vortex and impurity, E j,kin are shown in (b) and (c), respectively, and the two interaction energies E 11,vdW and E 12,vdW are shown in (d) and (e), respectively. In general there seems to be good agreement between the two energies amongst the various contributions to the total energy, with the exception of the kinetic energy contribution, panel (b). This small difference likely originates from approximating the length scale ξ v of the vortex using the chemical potential of the homogeneous system.
The variational approach can also be used to calculate an expression for the angular momentum of the superfluid system. The expectation value of the z-projection of the angular momentum is defined as L z = dr ψ * varLz ψ var . Using the representation ψ var = √ n var exp(iϑ var ) to separate the density n var and phase ϑ var variables, the angular momentum can then be computed asymptotically using Eq. (5a) and the assumption R/α 1 as The angular momentum of ψ 1 is calculated and compared in panel (f) of Fig. 3 with three different approaches. The first two data sets show L z computed using the Gross-Piteavskii solutions (blue plusses) and the variational angular momentum L z var , Eq. (8) (green circles). We can see that all three approaches are almost independent of N 2 , as we would expect, and are only a few percent away from the theoretical value L z Ana = N 1 (red crosses).

C. Vortex pairs
As well as individual vortices, it is also useful to consider pairs of vortices; either with the same sign (q 1,2 = ±1) or with differing signs (q 1 = 1, q 2 = −1) (see Eq. (4)). Then the resulting phase distribution for the vortex pair is where the vortices are separated by a distance x = d. By changing the separation d and mass imbalance ratio N 2 /N 1 , we can explore the effect that impurity atoms have. Figure 4 shows the computed impurity energies E 2 = drψ * 2 [−( 2 /2m)∇ 2 + g 12 |ψ 1 | 2 ]ψ 2 /N 2 for different vortex charges and mass imbalances as a function of vortex separation. Note that the asymptotic energy E 2 (d ξ) (i.e. the energy of two well separated vortices) has been subtracted for a given atom number N 1 to aid comparison. At large mass imbalances N 2 /N 1 = 1/4×10 3 there is a small splitting between the like-sign and vortex dipole energies (blue data, (a)) while for larger imbalances the splitting is more pronounced N 2 /N 1 = 1/500 (green data, (a)) and N 2 /N 1 = 1/250 (red data). We attribute this difference to the short-range nature of the vortices densities, which for a vortex dipole manifest as a depression when the vortices approach, resulting in a lower impurity energy than for a pair of like-sign vortices. Then, the three panels (i)-(iii) show example stationary states obtained from the solution of Eqs. (1). Here (i) and (ii) show the density |ψ 1 | 2 and phase ϑ 1 of the first component, respectively; while (iii) shows the delocalized impurity density |ψ 2 | 2 . These calculations reveal a simple relationship between the mass imbalance N 2 /N 1 and the splitting of the impurities energy, with smaller mass imbalances leading to larger differences in the impurities energy at short separations. Since we know that a vortex pairs energy is logarithmic in their separation d such that E 12 ∝ q 1 q 2 ln(R/d), larger vortices will comparatively feel an increased repulsion or attraction at the same distance d than those with smaller core sizes, consistent with the known behaviour of vortices. Figure 4 demonstrates that the presence of impurities will contribute an attractive interaction between vortices.

IV. COVALENT BONDS BETWEEN VORTICES
In this section we consider the possibility of bonding between a pair of like-sign vortices. Motivated by our findings in the previous section where the interaction between two vortices with both like and different signs were considered, we consider how the excited states of such a system can contribute an additional attractive force, such that the impurity mediates a bound-state between the vortex pair.

A. Padé approximation for a single vortex
Let us begin by writing the time-independent form of Eq. (1) using the transformation ψ j (r, t) = exp(−iN j E j t/ )φ j (r), giving and φ 2 → φ 2 /ξ. We then perform one additional scaling, φ j → N j φ j which allows us to understand the parameter regimes in terms of the coupling between the superfluid and the impurity. This yields the coupling parameters (corresponding to the dimensionless coefficients of |φ 1 | 2 and |φ 2 | 2 in Eq. (11) and (10), respectively) as C 1 = 2mg 12 ξ 2 n 1 N j / 2 and C 2 = 2mg 12 N 2 / 2 , while the units of energy are scaled as 2mN j ξ 2 E j / 2 for the energy of component j. The ratio of coupling constants yields C 2 /C 1 = (N 2 /N 1 )(1/n 1 ξ 2 ), which in the limit N 1 N 2 allows us to decouple Eqs. (10)- (11). Then, by separating the phase and amplitude of the vortex as φ 1 = N1E1 g11 e iqθ f (ρ) we obtain a simplified form of Eq. (10) To build the Padé approximation, we can solve Eq. (12) in two limits subject to the boundary conditions f (ρ → 0) = 0 and df (ρ)/dρ| ρ→∞ = 0. This gives the pair of asymptotic solutions here the constant a appearing in Eq. (13) is the first fitting constant used to build the Padé approximation [84]. We note that the authors of [25] used a piece-wise function to perform a variational analysis of the vortex's allowed bound states. Such an approach is useful for analytical treatments, however since we wish to obtain the excited states numerically, the Padé approximation provides an accurate approximation to the numerical solution of the Gross-Pitaevskii model, and hence a closer agreement to the true numerical spectrum of the impurity. Then, an appropriate approximation for a single vortex is written as which satisfies both asymptotic boundary conditions given by Eqs. (13) and (14). Figure 5 shows a comparison between the numerical solutions to Eq. (12) (blue solid line) and the Padé approximation, Eq. (15) (open orange circles). The two fitting parameters are found to be a ∼ 0.81 and b ∼ 0.88. Accompanying this, the (scaled) solution φ GPE 1 obtained from Eqs. (10)-(11) with N 1 = 500, N 2 = 1 and C 2 = 2mg 12 N 2 / 2 = 1 is shown for comparison (left vertical axis). The dashed grey data shows the difference = (|f num | − |f P |) 2 /|f num | 2 (right vertical axis) between the numerical and Padé approximations, which being small, supports this choice.

B. Impurity potential model
Having the analytic Padé approximation of Eq. (15) for f P (ρ), we can proceed to solve Eq. (11). Since we assume (N 2 /N 1 )(1/ξ 2 n 1 ) is small, we can treat φ 1 in Eq. (11) as an external potential by omitting any back reaction between the vortex and the impurity. We are left with the two-dimensional quantum mechanical problem of calculating the bound states of φ 2 of an impurity by the axially symmetric Schrödinger potential g 12 f P (ρ) 2 . Then we have here we have changed m → m i , the impurity mass to distinguish it from the impurity angular quantum number defined below, and the potential U(ρ) = g 12 |φ P 1 | 2 − N 1 E 1 /g 11 satisfies U(ρ → ∞) = 0 while the energy is E = N 2 E 2 − N 1 E 1 g 12 /g 11 . Then, the Schrödinger potential U(ρ) obeys U(0) = −N 2 E 2 g 12 /g 11 , U(ρ → ∞) = 0, hence the potential exists on the interval −N 2 E 2 g 12 /g 11 ≤ U(ρ) < 0. Since the potential U(ρ) is axisymmetric, we can separate ρ and θ as here m ∈ Z defines the quantum number of angular momentum for the impurity, with the associated operator L z = −i∂ θ . Normalization for ϕ(ρ) is then imposed through 2π Eq. (16), we are left with the following one-dimensional Schrödinger problem where we have introduced the operator then we haveQ †Q = −( 2 /2m)(d 2 /dρ 2 + 1/4ρ 2 ), while the term m 2 2 /2m i ρ 2 in Eq. (18) defines the centrifugal potential felt by the impurity. To understand the basic properties of the impurity, let us first consider m = 0. Note that the operatorQ †Q is positive semi-definite. Therefore, when U(ρ) vanishes (g 12 → 0), the eigenvalue E satisfies E g12→0 ≥ 0. The lowest eigenstate (n = 1) has E = 0 and should satisfyQϕ = 0. This can be easily solved by ϕ(ρ) ∝ √ ρ. This in turn corresponds to constant φ 2 , which is clearly non-normalizable. Now we consider small but finite g 12 by slightly increasing g 12 from zero. This has the effect of lowering the Schrödinger potential U(ρ) everywhere since U(ρ) < 0. Then, because we can continuously change g 12 from 0 to g 12 > 0, there is a corresponding continuous shift of the eigenvalues. Namely, the ground state energy should satisfy E 1,0 < 0 for g 12 = 0. Therefore, the Schrödinger equation at ρ → ∞ reads d 2 ϕ 1,0 /dρ 2 = 2m i |E 1,0 |ϕ/ 2 , the solution for which is a bound state whose asymptotic wave function is To summarize, we have shown that there exists at least one bound state regardless of the form of the potential U(ρ). On the other hand, the existence of (bound) excited states ϕ n,0 (ρ) for n > 0 sensitively depends on g 12 . In general the greater g 12 is, the more excited states exist and hence there are more possibilities to form bound states with a particular potential. On the other hand, for m = 0, the centrifugal potential always lifts the Schrödinger potential since it contributes a positive energy to the total potential felt by the impurity. Therefore, a sufficiently large g 12 is required to accommodate bound states for ϕ n,0 (n > 1).

D. Covalent bonds between two vortices
In this subsection we compute the excited states of two neighboring vortices. We consider two vortices placed at (x, y) = (±d/2, 0). For the background configuration of φ 1 , we take an product Ansatz φ (2) 1 (x, y, d)= and we solve Eq. (16) with the replacement U(ρ) → g 12 (|φ The presence of one vortex af-fects the other vortex, and the distortion cannot be ignored when the separation d is small. Therefore, in what follows we will take d ξ for the above product Ansatz to be a valid approximation. As explained in Sec. IV B, we omit any back reaction between the field φ 1 and φ 2 by keeping (N 2 /N 1 )(1/ξ 2 n 1 ) 1. In the previous subsection IV C, it was explained that each background φ 1 vortex represents an attractive potential for the impurity atoms. Therefore, when we place two such vortices on the plane, the impurities are attracted by such a pair of vortices. However, when the displacement d ξ, the individual vortex wave functions no longer overlap, so the energy spectra are the same as those of a single vortex, On the other hand, when the vortex separation d is large but finite, the degeneracy is lifted and instead the spectra are slightly split. Since the problem is no longer axially symmetric, we consider instead Eq. (11) with φ 1 replaced by φ (2) 1 , from which we obtain the numerical solutions. Let us start with the simple case of C 1 = 1 for which each vortex can have only a single 1s orbital. We vary d/ξ from d = 3ξ to d = 12ξ with step size δd = ξ/2. For each given d, we numerically solve Eq. (11) and determine the energy eigenvalues. The results are summarized in Fig. 7. As expected we have two states, a state with even symmetry, while the other is an odd function for the parity x/ξ → −x/ξ. The former is the ground state whereas the latter is the excited state, since it has one linear node on the y axis. The two states are almost degenerate when d is sufficiently large (see d 10ξ of Fig. 7). These states are slightly below the ground state energy of the single vortex, E 1s , for a large but finite d, and asymptotically approaches E (1) 1s from below. Therefore, at asymptotically large d, having impurities for both the ground state and the ex-  1s . The blue disks correspond to states which are even (no nodes) for x/ξ → −x/ξ, while the orange square are for those which have odd parity. We show the two wave functions for d/ξ = 3.5. The spatial region of the density plots is x/ξ ∈ [−10, 10] and y/ξ ∈ [−10, 10].  cited state results in an attractive force between the two background vortices of φ 1 . By decreasing d ( 9ξ), however, the ground and the first excited states show different behaviors. The ground state energy is a monotonically decreasing function of d in the limit d → 0. The decrease of the ground state energy occurs because the impurity wave function can localize around both of the vortices and its length scale increases compared to the single vortex background. Namely, sharing the impurity between two vortices gives rises to an attractive force. This is nothing but a covalent bonding-like effect for the vortices. Since the (1, 0) wave function is axially symmetric, we may call this the σ s-s bond. Then, we find a numerical fit of the covalent energy of the ground state as a function of d ( 5ξ) as On the other hand, the energy eigenvalues of the first excited state increase as d is decreased. This is the so-called anti-bonding effect common for states of odd parity. The presence of a linear node prevents the wave function from becoming shallow, instead the wave function becomes steeper in this region. Next, we consider the case of 2mg 12 ξ 2 n 1 N 1 / 2 = 5. As is shown in Fig. 6, each background vortex can localize the impurity not only as the ground state (1s) but also several higher excited states (2p, 2s, 3d, 3p, 3s). In general, a wave function of an excited state has a larger length scale compared to the ground state wave function, so covalent bonding effects start to appear for vortices placed with a larger separa-tion compared to states with lower energy. Furthermore, when a bound state has non-zero angular momentum (m = 0), its wave function is not axially symmetric but has directivity. As we will see in subsection IV E, these properties of the higher excited states will give us additional contributions to the covalent bonding effect. Let us begin with the ground state and first excited state, namely the 1s-1s mixtures. The numerical results are shown in Fig. 8. Qualitatively, there is no difference from the case with 2mg 12 ξ 2 n 1 N 1 / 2 = 1. The two states are almost degenerate in energy when the separation d is sufficiently large. As d decreases, both states decrease in a similar fashion from the asymptotic energy E (1) 1,0 −2.85. Then, when d is small enough, the spectrum undergoes splitting. The ground state energy continuously lowers due to the σ s−s bonding, while the energy of the first excited state stops lowering, and instead starts increasing due to the anti-bonding effect. As discussed above, the wave functions for larger g 12 are smaller than those for smaller g 12 . Therefore, the spectral splitting occurs at relatively small d ( 4.5ξ) in comparison with the case of C 1 = 1.
Although we studied the impurity bound states and bonding of vortex-vortex states, the choice of the vortex charge q will in fact not affect the bound-states, since the impurity feels the condensate density through a potential that depends on the density of the first component, (see Eq. (16)). However, the dynamics of the coupled system would be quite different depending on the choice of q.

E. Bonding at finite angular momentum
Let us next look at higher excited states (2p) in which two independent states (n, m) = (1, ±1) can coexist. Fig. 9 shows the numerical results. When the vortex-vortex separation d ξ, there are four degenerate states. Each state has directivity like a dipole, so that a total of four states can be geometrically constructed, (see right-most four panels of Fig. 9). We denote the individual 2p wave functions by the states |2p, ← , |2p, → , |2p, ↑ , |2p, ↓ expressing the direction of individual dipole axes. We consider the following products (viz. Eq. (21)) of these four states When the two vortices are well separated, the four states given by Eq. (23a)-(23d) are almost degenerate. However, they are not exactly degenerate as long as the distance d is finite. Then, ordering these four states by their energy, we have the lowest (second excited) state |• , Eq. (23a) (here the arrows indicate the direction on the bond axis, a blue disk in Fig. 9), the second from below (the third excited) state is denoted by | (orange triangle, Eq. (23b)), while the third from below (fourth excited state) state is denoted by | where the arrows are orthogonal to the bond axis (yellow squares, Eq. (23c)). The fifth excited state is denoted by | where the arrows are orthogonal to the bond axis (green diamonds, Eq. (23d)).
The eigen energies of the four states simultaneously decrease as d decreases with respect to their energy, with their order unchanged. However, on the interval 6 < d/ξ < 6.5, level crossings occur. Namely, the state | becomes the third from fourth while the state | is lifted from the third to fourth state, (see inset of Fig. 9).
The latter state | is further lifted and becomes the fifth excited state as d decreases. The state |• does not participate in mixing, remaining at the lowest level. The |• states energy suddenly decreases, caused by the increasing overlap of the individual lobes of each vortex. We denote this σ p−p bonding. Following this, the remaining three states also split around d/ξ 5. The energy of the second lowest state | also continuously decreases since two lobes of one orbital can hybridize with two lobes of the other. We denote this π p−p bonding. On the contrary, overlaps of lobes cannot occur due their geometry for the remaining two states | and | . For these states smaller displacements d cause their respective energies to increase due to the anti-bonding effect.
Finally, we explain the higher excited states involving the 2s and 3d orbitals, which are shown in Fig. 10. The two states 2s are denoted by and the interchanged labels appearing in the first ket on the r.h.s of Eq. (25) indicates the odd parity version of |2s, + , a general notation we adopt for the remainder of this work. When the vortex displacement d is large, the mixtures of the two 2s orbitals have almost degenerate energies which asymptotically converge to 2mN 2 ξ 2 E (1) Since the 2s orbitals are axially symmetric, the spectral flow of these almost degenerate states are qualitatively similar as for the 1s-1s bound states explained above. However, a difference does arise from the length scales of the wave functions. Since the 2s orbital is larger than the 1s one, the splitting of the two degenerate states due to the σ s−s anti-bonding begins at relatively large d ( 9ξ) displacement for the 2s state. However, by carefully examining the spectral flow, we find that the σ s−s bonding occurs in two stages. The first is due to overlap of the wave functions outer shells (d 9ξ), while the second occurs instead due to the overlap of the wave functions inner shell (d 3.5ξ). On the other hand, two vortices with 3d orbitals have four almost degenerate states whose energies asymptotically converge to mN 2 ξ 2 E (1) 3d / 2 −0.18. Since the 3d orbital has four lobes, we symbolically express these states as |3d, + and |3d, × . Amongst these four 3d states, the lowest energy state is denoted where a single lobe of one orbital facing to a single lobe of another has common color (an orange square in Fig. 10). Then, the second lowest energy state is where two lobes of one orbital are facing towards the two lobes of another, and have a common colour (green diamond in Fig. 10). The energy eigenvalues of both of these states decreases when the vortex separation d goes down towards zero. On the other hand, the odd parity states correspond to brown circles, while correspond to open blue squares. The energy of the |•, − and | , − states are lifted by the anti-bonding effect in the region where d is small. In the asymptotic limit (d ξ), the four states are almost, but not exactly degenerate. The observed ordering is E + < E + < E − • < E − . However, as d decreases, level crossings occur between the | , + and | , + states. In the region d 5ξ, the lowest energy state is | , + . The lowering of energy is due to the σ p−p bonding for | , + while for the | , + this effect originates instead from π d-d bonding. Therefore, the level crossing implies that the latter bonding is more efficient than the former one. In the region where d is small enough, these states are relatively smaller than the anti-σ s−s bonded 2s-2s mixture with odd parity. We also find another type of covalent bond, δ p−p which is a bound state of 3p orbitals (upside-down purple triangle) appearing second from the top (d 4ξ) in Fig. 10. This involves bonding among four lobes, each from a 3p orbital. Additionally, there are further higher excited states of 3p-3p and 3s-3s which we do not show here.

V. CONCLUSIONS
In this work, we have modelled the interaction of a superfluid vortex with quantum impurity atoms. This has been accomplished using a two component (binary) Gross-Pitaevskii formalism, where the superfluidimpurity interaction has been handled via a simple density-density interaction. The superfluid-impurity interaction has been characterized in the mass-imbalanced regime, and it was found that the presence of the impurity atoms can have a dramatic effect on the superfluid vortex; leading to distorted vortex profiles even at modest mass imbalances. The properties of vortex pairs have been also studied in the presence of impurities, and it was shown that there is a splitting of the impurities energy that depends on the signs of the vortex pairs. The splitting of the impurities energy has been observed to be greater at smaller imbalances.
In the second part of this work we considered the excited states of the coupled system. Using a Padé approximation, a model for the impurity potential generated by the vortex was constructed that allowed us to diagonalize the impurity component in the limit of small mass imbalances, subsequently allowing the examination of the excited states of the system. In general, the number of accessible excited states depends on the details of the coupling between the two components as well as the particular choice of mass imbalance. For deeper potentials, more excited states become accessible. The impurity orbitals of this two-dimensional system possess a two-dimensional hydrogen-like character, with higher excited states carrying finite angular momentum, and being more diffuse in character as the zero (ionization) energy limit is approached. The excited impurity orbitals were then used to construct the eigenstates of a pair of likesign vortices. The spectral flows of the impurities lowest energy excited states revealed an attractive power law behaviour, raising the possibility of covalent bonds between vortices. The spectral flows of the excited (angular momentum carrying) states also revealed a rich behaviour. Chemical bonding of pairs of vortices possessing one or two units of angular momentum display physical effects analogous to their chemical counterparts, such as the appearance of bonding and anti-bonding states, as well as spectral level crossings.
There are several interesting routes for future investigations. Extending this analysis to the case of fermions is a novel and worthwhile avenue, since this system is a good candidate for producing a novel matter-wavetrapped degenerate Fermi gas, analogously to the experiment of DeSalvo et al., Ref. [9]. The impurity states of this system also represent good candidates for quantum information applications, such as long-lived qubits analogously to those proposed for dark solitons [85,86]. Our analysis of the impurity mediated chemical bonds could also be used as a basis to construct more elaborate forms of controllable synthetic matter, such as Toda lattices [87], vortex-impurity crystals and the impurities associated superfluid-insulator transitions. Further, collisions between multiple vortex carrying impurities would allow for the simulation of synthetic chemical reactions. Our work also represents a new tool in the growing field of atomtronics [88], where the ability to construct synthetic chemical bonds yields a useful new paradigm in this emerging discipline.