Vortices with massive cores in a binary mixture of Bose-Einstein condensates

We analyze a notable class of states relevant to an immiscible bosonic binary mixture loaded in a rotating box-like circular trap, i.e. states where vortices in one species host the atoms of the other species, which thus play the role of massive cores. Within a fully-analytical framework, we calculate the equilibrium distance distinguishing the motion of precession of two corotating massive vortices, the angular momentum of each component, the vortices healing length and the characteristic size of the cores. We then compare these previsions with the measures extracted from the numerical solutions of the associated coupled Gross-Pitaevskii equations. Interestingly, making use of a suitable change of reference frame, we show that vortices drag the massive cores which they host thus conveying them their same motion of precession, but that there is no evidence of tangential entrainment between the two fluids, since the cores keep their orientation constant while orbiting.


I. INTRODUCTION
Vortices in quantum fluids are topological excitations characterized by quantized circulation [1] which are present in a number of nonlinear field theories and models [2], ranging from superfluid media [3,4] and quantum optics [5,6] to superconductivity theories [7,8] and Josephson-junction arrays [9,10] and play a key-role in fundamental effects such as superfluid turbulence [11], the Berezinskii-Kosterlitz-Thouless transition [12], fractional statistics [13], and in the development of a fully quantized field theory for topologically complex excitations [14][15][16]. Among the plethora of different physical systems where vortices can be experimentally investigated, ultracold quantum gases provide a particularly controllable and versatile platform [17,18] for the study and the observation of the rich phenomenology associated to their formation [19,20], dynamics [21,22], and interactions [23]. Vortices in Bose-Einstein condensates (BECs) were first obtained by means of a phaseimprinting method involving two hyperfine spin states of 87 Rb [24] but, at present, can be produced also by stirring the BEC above a certain critical velocity [25][26][27], dragging barriers through the BEC itself [28] or interfering multiple condensate fragments [29].
Solitons are a kind of localized excitations which, because of the competition between dispersion and nonlinearity, propagate keeping their shape unaltered, even if a two-soliton collision occurs [30]. Soon after the achievement of Bose-Einstein condensation, different types of solitons have been described and observed [18,31,32]. To our purposes, of particular importance are those systems where a bosonic binary mixture features dark-bright * Correspondence to: andrea.richaud@polito.it soliton configurations [33][34][35]. These structures, first predicted in Ref. [36], are frequently termed as "symbiotic solitons" [37] because the bright component, being endowed with repulsive intraspecies interaction, could not exist if the dark component did not play the role of an effective confining potential.
The same symbiotic relationship was shown to constitute the mechanism underlying the robustness of vortexbright soliton complexes. First observed in 2000 by the JILA group [38], they represent the topological extension of the dark-bright soliton configuration to the case where a component hosts one or more vortices [39]. The aforementioned study paved the way to a series of further investigations which highlighted, among various aspects, the spontaneous generation of vortex-bright soliton structures [40], the possibility, for the effective potential well corresponding to the vortex core, to support not only bound states [41], but also multi-ring excited radial state complexes [42], and a rich dynamical scenario for the bright-solitary component [43].
Within an analytical framework and by means of extensive numerical simulations, our work aims at analyzing the static and the dynamical properties of vortexbright soliton complexes, i.e. how the presence of massive solitons within the cores of two corotating vortices affect the equilibrium distance characterizing their motion of precession around the trap center, the role of the interspecies repulsion as an antagonist to the centrifugal force acting on the solitons, the functional dependence of the angular momenta carried by two species, of the vortex healing length, and of the characteristic radius of the massive cores on the mass of the latter.
If one considers repulsive intraspecies (g a , g b ) and interspecies (g ab ) interactions such that the immiscibility condition g ab > √ g a g b is fulfilled, the dynamical picture of the mixture [in which the order-parameter fields of arXiv:1908.06668v3 [cond-mat.quant-gas] 24 Nov 2019 the two species obey two coupled Gross-Pitaevskii equations (GPEs)], indeed reduces to the much simpler equations of two point-like vortices with nonzero-mass cores.
Noticeably, the latter are found to exhibit an evident Lorentz-like form since, in the presence of vortex cores occupied by a second species, the vortex-motion equations are equivalent to those of a pair of massive charges acted by a transverse magnetic field. With negligible fractions of the minority component, one recovers the Helmholtz-Kirchhoff equations for planar point-like vortices [44].
The outline of the manuscript is the following: in Sec. II, we present an analytical model for the dynamics of massive vortices in a confined system which incorporates the effect of the virtual vortices resulting from the boundary condition of vanishing normal velocity. In particular, we derive a formula giving the equilibrium distance distinguishing the motion of precession of two corotating massive vortices. Sec. III is devoted to the presentation of the two coupled stationary GPEs which provide a good description of the bosonic binary mixture in the meanfield approximation. In Sec. IV, we show how the presence of massive cores (i.e. species-b atoms trapped within species-a vortices) affects the equilibrium distance of the pair of corotating vortices. We also show that the interspecies repulsion tends to counterbalance the centrifugal force acting on species-b atoms. In Sec. V, we address the angular momenta of the two components and provide analytical formulas that well capture their functional dependence on the number of species-b atoms (which, in turn, is directly proportional to the mass of the cores). By means of a suitable change of reference frame, we show that the cores, although following the same motion of precession of the vortices, keep their orientation constant while orbiting. This circumstance witnesses the fact that there is no tangential entrainment between the two fluids. In Sec. VI, we present an heuristic but effective system of equations that well reproduces the functional dependence of the vortex healing lengths and of the cores' characteristic radius on the number of speciesb atoms. Eventually, Sec. VII is devoted to concluding remarks.

II. POINT-LIKE VORTICES IN A CIRCULAR BOX
In this section, we review some results concerning the dynamics of point-like vortices and we introduce a model for the dynamics of vortices whose cores host point-like masses (hence the name massive vortices, as opposed to the traditional massless vortices).

A. Massless vortices
The Hamiltonian of N point-like massless vortices in an ideal unbounded fluid is given by [44] where ρ * is the fluid planar density, z j = x j + iy j ∈ C is the position of the jth vortex in the ambient plane and k j = n j h/m f is its strength (n j ∈ Z is the vortex quantization and m f is the mass of the fluid particles).
In the following, we will specialize our discussion to the case of N = 2 vortices. When one considers bounded systems, Hamiltonian (1) modifies due to the presence of the confining potential. In the case of a box-like potential (this type of confinement is within the reach of current experimental trapping techniques, see, e.g., Refs. [45][46][47]), the presence of a boundary confining the fluid is accounted for by means of the virtual charge method, i.e. by introducing a suitable configuration of virtual vortices. With this premise in mind, the Hamiltonian of N = 2 point-like massless vortices in an ideal fluid confined in a circular box of radius R reads [48] In this framework, the coordinates of each vortex constitute a pair of canonically conjugate variables, and motion equations can be obtained by means of the Poisson Brackets involving, in turn, the canonical brackets {x i , y j } = δ i,j /(ρ * k j ) (see for example [49]).

B. Massive vortices
If one wants to introduce into the model the fact that the vortex cores host point masses, it is convenient to move to the Lagrangian formalism, where the presence of massive cores can be taken into account as follows where m j represents the point-like mass hosted by the jth vortex core and whereq j := dq j /dt (with q = x, y).
Note that, Lagrangian (3) is formally equivalent to that describing charged particles in a planar domain subject to a transverse magnetic field, where k j s and ρ * play the role of charges and of magnetic field, respectively. As is well known, the dynamics of the vortex cores is generated by the Euler-Lagrange equations which, in the special but interesting case of two equal vortices k 1 = k 2 = k, whose cores host two equal masses m 1 = m 2 = m, takes the form for i, j ∈ {1, 2} and i = j [ u 3 is the unit vector perpendicular to the plane (x, y)].
The resulting system of 4 differential equations admits a notable solution, provided that the two vortices are placed symmetrically with respect to the box-trap center and that their distance d and the angular frequency Ω marking their motion of precession fulfill the following equation: As expected, Eq. (4) shows a pathology when d → 2R, meaning that the vortex pair is approaching the circularbox boundary. Moreover, in the limit of infinite box radius (R → +∞), one can retain only those terms ∝ R 4 and the relation d(Ω) can be expressed in closed form, i.e.: In Sec. IV A, the equilibrium distance d predicted by Eq. (4) and relevant to two equal point-like massive vortices in a circular box will be compared to the one obtained by numerically solving two coupled stationary GPEs. To conclude this Section, we would like to remark that the extension of model (3) to the case of harmonic confinement is far from being trivial, as the unavoidable curvature of the enveloping wavefunction produces non-negligible effective forces acting on the vortices' centers, which distort the usual vortex dynamics already in the case of zero species-b atoms. It is indeed the use of a box-like potential that allows one to bypass the influence of the aforementioned non-negligible effective forces, thus allowing for a cleaner emergence of the genuine phenomenology characterizing vortex/bright-soliton complexes.

III. THE BOSONIC MIXTURE
We consider a bosonic mixture of 23 Na and 39 K [50,51]. Each atomic species is characterized by an order parameter, ϕ a = √ N a ψ a and ϕ b = √ N b ψ b respectively. In a mean-field treatment of the problem, we assume that the system is effectively quasi-2D, as a result of a strong confinement along the z-direction. Because of this, it can be effectively modeled by the following two coupled stationary Gross-Pitaevskii equations where N a (N b ) corresponds to the number of species-a (species-b) atoms, g c = 4π 2 a c /m c , with c = a, b are the intraspecies interaction strengths and g ab = 2π 2 a ab /m ab is the interspecies coupling. Notice that m a (m b ) is the atomic mass of sodium (potassium), while m ab = (m −1 a + m −1 b ) −1 is their reduced mass; similarly a a and a b are the intraspecies scattering lengths, while a ab is the interspecies scattering length. Parameter z is the effective thickness of the disk-like box trap and functions ψ a and ψ b are normalized to 1, since ϕ a and ϕ b are, respectively, normalized to N a and N b .
Vortical solutions of Eqs. (6) are found by moving to a frame rotating with angular velocity Ω (this corresponds to adding the term −ΩL z to the Hamiltonian, whereL z is the operator associated to the third component of the angular momentum) and then employing the imaginary-time method [41,52]. The starting condition for the imaginary-time dynamics is such that species-a hosts a vortex pair while species-b is localized (two narrow Gaussian distributions) at the vortex cores. As the fictitious dynamics advances, the position of the vortex cores, their healing length, together with the spatial distribution of species-b atoms is iteratively self-consistently refined, until convergence is reached. We conclude this Section by noticing that, for our purposes, solutions of Eqs. (6) of the type vortex/bright-soliton pairs (see, e.g. Fig. 1) could be either actual ground states or excited metastable states (i.e. either global or local minima in the energy landscape).

IV. MASSIVE VORTEX PAIRS IN A BINARY MIXTURE OF BECS
Eigensystem (6) was solved sweeping model parameter N b , the number of species-b atoms, which constitute the massive cores of species-a vortices. As explained in Sec. III, a suitable ansatz for the starting condition of the imaginary-time dynamics was chosen. In the whole range of N b that we explored (i.e. N b ∈ [5, 1000]), our numerical simulations [53] converged to a stationary state of the type illustrated in Fig. 1. Basically, condensate a is highly confined by the box-like potential (whose radius is R = 50 µm) and it is marked by the presence of two corotating vortices. The latter are symmetrically positioned with respect to center of the trap, they are such that the density |ψ a | 2 goes to zero in the center of the cores and feature a quantized circulation. On the other hand, component b occupies the vortex cores which, in turn, constitute an effective double-well potential for species-b atoms [43].
One can gain a further insight into the discussed eigensolution by computing the mass current density (with c = a, b) which also corresponds to the momentumper-particle distribution. As illustrated in the second row of Fig. 1, both vortices in species a rotate anticlockwise, thus determining a collective motion of precession which is anticlockwise too. As concerns species-b atoms, they are dragged by condensate a and remain bound within the vortex cores thus featuring their same motion of precession around the center of the trap. With reference to the middle right panel of Fig. 1, one can appreciate that the left (right) peak of |ψ b | 2 is translating downward (upward), i.e. along a direction tangential to the precession orbit.
Eventually, another important information that can be extracted from the eigensolutions ψ a and ψ b of equations (6) concerns the phase fields. The latter, denoted by θ a and θ b , and such that ψ c = |ψ c | 2 e iθc (where c = a, b), have been plotted in the lower panels of Fig. 1. As expected, θ a features singularities in correspondence of the vortices' centers, while, less expectedly, θ b features singularities too. Nevertheless, the latter are found far from the peaks of |ψ b | 2 . According to the basic properties of quantum fluids, the circulation In the first row, yellow (blue) is associated to large (zero) values of the density |ψ| 2 . Third row corresponds to the phasefield associated to ψa and ψ b [blue (yellow) color corresponds to −π (+π)]. Left (right) column corresponds to species a (b). The following parameters have been used: Na = 5 × 10 4 , of the velocity vector field associated to ψ c , that means v c = mc ∇θ c (with c = a, b), is zero if the closed path γ does not encircle singularities of the phase field θ c (with c = a, b). Conversely, due to the Feynman-Onsager quantization rules, it takes values n h/m c , where n ∈ Z, if one or more singularities of the associated field θ c are encircled by γ. This information will be useful to better understand the discussion about the angular momentum of species-b bosons (see Sec. V B).

A. Mass of the cores and equilibrium distance
Increasing the number of species-b atoms (within the investigated range N b ∈ [5, 1000]), the distance d vor between the centers of the vortices increases. Similarly, the distance between the two peaks of |ψ b | 2 , d peak , increases upon increasing N b . Fig. 2 shows how the presence of massive cores deforms and displaces the vortices. Notice that the position of the peaks of |ψ b | 2 does not exactly match that of the minima of |ψ a | 2 due to the centrifugal force on species-b atoms and the finite repulsive coupling (g ab < +∞) between the two fluids which, in turn, allows for a non-zero penetration of fluid b into fluid a. Therefore, observables d vor and d peak , which in the analytical model based on point-like vortices and point-like massive cores (see Sec. II) collapse on the same variable (d), when estimated from the numerical solution of Eqs. (6), despite being closely related, do not necessarily coincide.
The functional dependence of d vor and d peak [extracted from the numerical solutions of system (6)] on N b is illustrated in Fig. 3, together with the relation d(N b ), obtained, in turn, by means of substitutions into Eq. (4). Relations (8) allow one to match the analytical model (4) based on point-like vortices and point-like massive cores with the actual parameters used to model the quantum fluids within the mean-field approach [see system (6)].
The agreement between the analytical prevision (yellow dotted line) and numerical results (blue solid and red dashed lines) is remarkably good, both qualitatively (same quasi-linear functional dependence on N b ) and quantitatively (offset < 2%). Moreover, we would like to mention that numerical results (namely, the slope and the vertical shift of the corresponding lines of Fig. 3) can be shown to further approach the analytical prevision upon increasing N a and/or diminishing g b , two changes that result in narrower cores and, therefore, in a scenario where Eq. (4), based on our point-like approximation, reliably describes the mixture vortex state.

B. Competition between centrifugal force and interspecies repulsion
As already mentioned, d peak , although closely related to d vor , is always slightly bigger than the latter. The motion of precession of the vortices around the center of the trap is responsible, in fact, for a centrifugal force on species-b atoms which are, therefore, pushed outwards. This tendency is only partially opposed by the repulsive interaction between the two quantum fluids and it is the competition between these two forces what determines the exact values of d vor and d peak . Increasing the interspecies repulsion g ab , fluid a gets more impenetrable to Upper panel: the position of markers '+' corresponds to dvor/2, while the distance between markers '×' corresponds to 2ξa. Lower panel: the position of markers '×' corresponds to d peak /2, while the distance between markers ' * ' corresponds to 2σ b . The following parameters have been used: Na = 5 × 10 4 , Ω = 5 rad/s, R = 50 µm, ma = 3.82 × 10 −26 kg, m b = 6.48 × 10 −26 kg, ga = 52 × (4π 2 a0)/ma, g b = 7.6 × (4π 2 a0)/m b , g ab = 24.2 × (2π 2 a0)/m ab , z = 2 µm. species-b atoms, which therefore prove to be more tightly bound within the valleys of |ψ a | 2 . As a result of this increased reaction to the centrifugal force, the difference d peak − d vor is remarkably smaller, as illustrated in Fig.  4 (where g ab has been set 2.5 times bigger than the value used for Fig. 3).

C. Stability beyond the immiscibility condition
We remark that the results illustrated in this article are obtained under the assumption that the two quantum fluids are immiscible, meaning that their intra-and interspecies coupling parameters are such that g ab / √ g a g b > 1.
Actually, we should mention that vortex/bright-soliton complexes prove to be rather robust composite objects, as FIG. 4. Equilibrium distance: comparison between numerical (dvor and d peak ) and analytical (d) results. For these simulations, parameter g ab is 2.5 times bigger than the one used for Fig. 3, all the others being unchanged.
their stability extends also beyond the immiscible regime. We have verified this by numerically solving Eqs. (6), thus extending the results of Ref. [34], derived in the context of dark-bright soliton complexes, to the case of vortex/bright-soliton complexes. The obtained density distributions, |ψ a | 2 and |ψ b | 2 , are qualitatively similar to the ones illustrated in Fig. 1 and obtained in the case of immiscible components. Of course, the difference between the miscible and the immiscible regime is that, in the latter case, bright solitons are more tightly confined within the vortex cores than in the former case. Moreover, due to the outward shift of the bright soliton from the vortex center, ensuing from the centrifugal force, the soliton itself, at a given time, feels an effective elliptical potential which, in turn, deforms its circular shape. Concerning our sample, which is immiscible, the deviation of the shape of bright solitons from a perfectly circular one can be neglected at first approximation, as the associated ellipticity, although being a decreasing function of N b , is always ≈ 0.98 in the whole considered range of N b . When the two quantum fluids are miscible, in fact, bright solitons manage to invade the majority component in a more significant way. As a consequence, they play the role of less rigid (i.e. softer) massive cores. If, on one hand, this circumstance extends the robustness of these composite objects to the case of miscible quantum fluids, on the other hand, in the miscible regime, our point-like model (3) partially looses its validity, as the mass of the bright solitons is not concentrated in the vortices centers any more, but it spreads and occupies more extended spatial regions.

V. ANGULAR MOMENTUM OF VORTICES AND CORES
This section is devoted to the analysis of the angular momentum of each component, an investigation that can offer a deeper insight into the Physics of the system. In particular, we show that the two massive cores (made of species-b atoms) orbit around the center of the trap, being dragged by the motion of precession of the vortices. Nevertheless, they do not rotate, i.e. their orientation remains constant while they revolve.

A. Angular momentum of condensate a
The angular momentum (per particle, in units of ) of condensate a can be computed as This quantity can be evaluated numerically from the solution of Eqs. (6). On the other hand, it can also be estimated by means of a fully-analytical approach. Along the same lines discussed in Ref. [54] (where the authors investigated the case of harmonic confinement), in fact, it is possible to derive the following expression: where r vor := d vor /2 constitutes the orbit radius. As shown in Fig. 5, Eq. (10) well fits the numerical data obtained by means of Eq. (9), the mismatch being < 0.8%. In this regard, it can be shown that the fitting accuracy further increases if one increases N a and/or decreases g b , because, in this case, the point-like approximation of vortices and cores gets increasingly valid. The angular momentum (per particle, in units of ) of component b can be analogously computed as a quantity that can be evaluated numerically on the basis of the solution of Eqs. (6). As already mentioned, the two species-b cores orbit around the center of the trap but they do not rotate around their own centers of mass. To prove this statement, we proceed along three different lines.
a. Mass current density in the rotating frame. In the lab frame, it is possible to compute the mass current density J b associated to ψ b [see Eq. (7)]. The corresponding vector field is illustrated in the middle right panel of Fig. 1. It is clear that the left (right) core is moving downward (upward), dragged by the anticlockwise motion of precession of the vortices. Due to the characteristic magnitude of | J b |, this plot does not allow one to understand whether the cores change their orientation or not along their orbit around the center of the trap. To circumvent this limitation, we have computed the species-b mass current density in a (non-inertial) frame rotating with the same angular velocity Ω distinguishing the motion of precession of the vortices. More specifically, in the rotating frame, J b,rot reads where and V = Ω ∧ r. J b is numerically computed through Eq. (7). Notice also that, at t = 0, |ψ b | 2 ≡ |ψ b,rot | 2 , thus justifying its use in Eq. (12). The result of this procedure is illustrated in the upper panel of Fig. 6 which shows that the two species-b cores, when observed from the rotating frame, rotate around their respective centers of mass, with angular velocity −Ω (the minus sign being due to the clockwise direction). On top of that, we have evidenced how these two cores rotate almost as if they were rigid bodies, meaning that the (absolute value of the) velocity field v b,rot around each center of mass linearly increases with the distance r C from the respective center of mass (this is not in contrast with the irrotational property of quantum gases, as the new reference frame is not inertial). In view of the symmetry of the ground-state configuration depicted in Fig. 1, we refer to the left (right) center of mass when considering the velocity field in the half-plane x < 0 (x > 0). The lower panel of Fig. 6 shows the local angular velocity of species-b cores when observed from the rotating frame. This quantity, defined as [where r C is the distance of point (x, y) from the left (right) center of mass when x < 0 (x > 0)], takes the (almost) constant value ≈ −5 rad/s in the most part of the regions where |ψ b | 2 is non-zero.
As an alternative indicator, in order to investigate the possible rotational properties of bright solitons, one could employ the vorticity distribution in the rotating frame w b,rot = ∇ ∧ v b,rot . Apart from two Dirac-delta-like singularities exactly where phase singularities are (see lower right panel of Fig. 1), one would observe two quasiplateaus at | w b,rot | = 2Ω b,rot ≈ −10 rad/s. This circumstance is in great agreement with the fact that the vorticity distribution associated to a rigid body rotating with angular frequency Ω is uniform and equal to 2Ω.
In conclusion, we have proved that, in the (noninertial) rotating frame, the two species-b cores rotate around their respective centers of mass with (almost uniform) angular velocity −Ω. This allows us to conclude that the they keep their orientation fixed when observed from the lab frame.
b. Analytical estimate of the angular momentum. To corroborate what elucidated in the previous paragraph, we show that the functional dependence of quantity (11) on model parameter N b can be well fitted by the semi-analytical model whereL and where terms are introduced to take into account that, in the lab frame, the two species-b cores revolve but keep their orientation fixed [Heaviside functions Θ(−x) and Θ(x) allow one to select the left and the right core respectively]. The integrals in expressions (15)- (17) represent three moment of inertia corresponding, respectively, to the anticlockwise revolution of the whole system around the center of the trap O, and to the effective clockwise rotation of the left (right) core around its own center of mass C L (C R ) [in this regard, r 2 := x 2 + y 2 , r 2 Cα := (x − x Cα ) 2 + (y − y Cα ) 2 , with α = L, R]. The latter effective motions indeed compensate for the fact that a pure motion of revolution [captured by Eq. (15)] would determine a change in the orientation of the cores along the circular orbit. Figure 7 shows a very good agreement between Eq. (11) and formula (14), the error being always < 4%. c. Observation of the associated phase field. As illustrated in the lower right panel of Fig. 1, the phase field associated to ψ b , θ b , features two singularities far from the centers of bright solitons. If one considers a closed path γ encircling one of the bright solitons, since no singularities of θ b are surrounded by it, the circulation C γ [ v b ] of the velocity field v b along γ must be zero (see Sec. III). This implies that both bright solitons do not rotate around their own axes (in the laboratory inertial frame). On the other hand, these singularities in θ b indeed play an important role because the corresponding velocity field v b = m b ∇θ b is what determines the collective precession motion of species-b bosons.
In summary, we notice that the interspecies repulsive coupling g ab is the interaction underlining the dragging of species-b cores by species-a vortices, which therefore exhibit the same motion of precession. In spite of this precession, species-b cores keep their orientation constant (in the inertial frame of the laboratory). Due to the irrotational properties of quantum fluids, in fact, a solitonlike distribution cannot be put in rotation around its own axis (as if it was a rigid body) without creating a phase singularity at its center. On the other hand, the creation of such a phase singularity would turn the soliton-like original distribution into a vortical object. To conclude this Section, we would like to mention the possible existence of the Andreev-Bashkin effect [55][56][57][58], according to which, the mass current density J i (with i = a, b) in one species will depend, in general, also on the velocity v j (with j = b, a) of the other species. In other words, the condensate density ρ ij is a non-diagonal matrix, a circumstance which implies the relations At the microscopic scale, this drag between mass current densities comes from the formation of quasi-particles with non-zero content of mass for either of the two components [58]. As a consequence, the transport properties of the two quantum fluids turn out to be coupled: the flow of one species influences the mass transport in the other species [55]. This effect is known to be rather elusive [58] and, for our sample condition, we expect it to be negligible, as the off-diagonal matrix elements ρ ab and ρ ba , which depend on the the overlap between ψ a and ψ b , should be small, given that the two discussed fluids are immiscible. In view of its complexity, the possible presence of this effect in the discussed system of orbiting vortex/bright-soliton complexes will be analyzed in a future work.

VI. VORTEX HEALING LENGTHS AND SIZE OF THE MASSIVE CORES
The presence of species-b massive cores within speciesa vortices affects the healing length of the latter. The intraspecies repulsive interaction, in fact, tends to enlarge the cores which, in turn, tend to swell (from the inside) the profile of the vortices because of the interspecies repulsive coupling. Flipping the perspective, the expansion of the cores is dammed by the species-a fluid, which plays the role of an effective confining potential for species-b atoms.
In the attempt to estimate the equilibrium healing length ξ a of species-a vortices and the equilibrium characteristic size of species-b cores, σ b , we present the following heuristic equations 2 2m a 1 ξ 2 a = +g a n a − g ab n a n b π( where n a = N a πr 2 box z , Notice that the the first equation of system (18) reduces, in the case of no interspecies interaction, to ξ a,0 = 2 2m a g a n a , the well-known formula derived in the context of singlespecies vortices [59]. Similarly, the second equation, if g ab = 0, gets structurally similar to formula giving the characteristic size of a soliton in the case of attractive interactions [59] (of course, in this context, n b,0 represents the central density). The extra term g ab n a n b π(σ 2 b − ξ 2 a ) z is introduced to take into account the interspecies repulsion, an interaction that manifests only in those regions where ψ a and ψ b overlap, i.e. only in the two annuli centered in the vortices' centers and whose outer and inner radii are σ b and ξ a respectively.
In order to compare the previsions provided by equations (18) with the values extracted from the numerical solutions of system (6), one has to give the operational definition of "vortex healing length" and "core characteristic size". From the numerical side, with reference to Fig. 2, we agree to measure the half width at half maximum of the valley of |ψ a | 2 , λ a , and the half width at half maximum of the peak of |ψ b | 2 , λ b . From the analytical side, the estimates of quantities λ a and λ b are given by the solutions of system (18), ξ a and σ b , multiplied by two suitable constant conversion factors, 1.30 and 1.15 respectively, which are determined from their numerical counterpart in the case N b = 1, that means in a scenario where species-b cores have a negligible impact on species-a vortices.
As illustrated in Fig. 8, equations (18) well capture the functional dependence of λ a and of λ b on N b .

VII. CONCLUDING REMARKS
In this work, we have investigated a notable class of configurations exhibited by a bosonic immiscible binary mixture loaded in a box-like circular trap, namely, minimum-energy states where species-b atoms are trapped within the vortex cores of species-a fluid. Both within a fully-analytical framework and by means of a systematic analysis of the numerical solutions of the associated two coupled GPEs, we have shown that the presence of massive cores alters the equilibrium distance distinguishing the motion of precession of the vortex pair. Interestingly, for the considered choices of model parameters (repulsive intra-and interspecies interactions such that, in the homogeneous case, the miscibility condition g ab < √ g a g b is not met) the dynamical mean-field picture of the mixture has been shown to reduce to much simpler effective equations exhibiting an evident Lorentzlike magnetic form, where massive vortices play the role of massive charges confined on a plane and subject to a magnetic field. Species-b cores, in turn, are dragged by fluid a and thus follow their same motion of precession around the trap center; nevertheless, while orbiting, they keep their orientation constant, meaning that there is no tangential entrainment between the two fluids. We have also derived, in the context of the Thomas-Fermi approximation, a simple formula to estimate the angular momentum of condensate a and we have shown, by means of a suitable change of reference frame, that species-b cores effectively behave as two rigid bodies. Eventually, we have introduced a system of heuristic but effective equations to estimate the characteristic size of vortices and cores hosted therein.