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

TL;DR

This paper analytically and numerically investigates vortices with massive cores in a binary Bose-Einstein condensate mixture, revealing their precession dynamics and core behavior in a rotating trap.

Contribution

It provides a fully analytical framework for understanding vortex-core interactions and compares these predictions with numerical solutions of coupled Gross-Pitaevskii equations.

Findings

Vortices precess with a calculable equilibrium distance.

Massive cores follow the same precession motion as vortices.

No tangential entrainment between the two fluids was observed.

Abstract

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.

Equations

$$H_{\infty }=(z_1,\dots ,z_N)=-\frac{{\rho }_{\ast }}{4\pi }\sum _{i=1}^N\sum _{j\ne i}{k}_i{k}_j\ln \frac{|z_i-z_j|}{\lambda },$$

$$H=\frac{{\rho }_{\ast }}{4\pi }\{k_1{k}_2\log \frac{|R^2-z_1{\bar{z}}_2{|}^2}{|R(z_1-z_2){|}^2}$$

$$+k_1^2\log \left(1-\frac{|z_1{|}^2}{R^2}\right)+k_2^2\log \left(1-\frac{|z_2{|}^2}{R^2}\right)\}.$$

$$\{F,G\}=\frac{1}{{\rho }_{\ast }{k}_j}\sum _{j=1}^2\left[\frac{\partial F}{\partial x_j}\frac{\partial G}{\partial y_j}-\frac{\partial G}{\partial x_j}\frac{\partial F}{\partial y_j}\right]$$

$$L=\sum _{j=1}^2\left[\frac{m_j}{2}({\dot{x}}_j^2+{\dot{y}}_j^2)+\frac{k_j{\rho }_{\ast }}{2}(y_j{\dot{x}}_j-x_j{\dot{y}}_j)\right]-H,$$

$$m{\ddot{\vec{r}}}_j=k{\rho }_{\ast }{\vec{u}}_3\wedge {\dot{\vec{r}}}_j+{\rho }_{\ast }\frac{k^2}{2\pi }[\frac{{\vec{r}}_j-{\vec{r}}_i}{|{\vec{r}}_j-{\vec{r}}_i{|}^2}+\frac{{\vec{r}}_j}{R^2-r_j^2}+$$

$$\frac{R^2{\vec{r}}_i-r_i^2{\vec{r}}_j}{R^4-2R^2{\vec{r}}_i{\vec{r}}_j+r_i^2{r}_j^2}],$$

$$x_1(t)=\frac{d}{2}\cos (\Omega t),y_1(t)=\frac{d}{2}\sin (\Omega t)$$

$$x_2(t)=\frac{d}{2}\cos (\Omega t+\pi ),y_2(t)=\frac{d}{2}\sin (\Omega t+\pi ),$$

$$\frac{\pi d^6\Omega (k{\rho }_{\ast }-m\Omega )+3d^4{k}^2{\rho }_{\ast }}{d-2R}=$$

$$-\frac{16\pi d^2{R}^4\Omega (m\Omega -k{\rho }_{\ast })+16k^2{\rho }_{\ast }{R}^4}{d-2R}.$$

$$d=k\sqrt{\frac{{\rho }_{\ast }}{\pi }}\frac{1}{\sqrt{k{\rho }_{\ast }\Omega -m{\Omega }^2}}.$$

$$-\frac{{\hslash }^2}{2m_a}\left[\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right]{\psi }_a+\frac{g_a{N}_a}{{\ell }_z}|{\psi }_a{|}^2{\psi }_a$$

$$+\frac{g_{ab}{N}_b}{{\ell }_z}|{\psi }_b{|}^2{\psi }_a+V_{ext,a}{\psi }_a={\mu }_a{\psi }_a$$

$$-\frac{{\hslash }^2}{2m_b}\left[\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right]{\psi }_b+\frac{g_b{N}_b}{{\ell }_z}|{\psi }_b{|}^2{\psi }_b$$

$$+\frac{g_{ab}{N}_a}{{\ell }_z}|{\psi }_a{|}^2{\psi }_b+V_{ext,b}{\psi }_b={\mu }_b{\psi }_b$$

$${\vec{J}}_c=-\frac{i\hslash }{2}({\psi }_c^{\ast }\nabla {\psi }_c-{\psi }_c\nabla {\psi }_c^{\ast })$$

$${\mathcal{C}}_{\gamma }[{\vec{v}}_c]={\oint }_{\gamma }{\vec{v}}_c\cdot \mathrm{d}\vec{r},(c=a,b)$$

$$k=\frac{h}{m_a},{\rho }_{\ast }=\frac{N_a{m}_a}{\pi R^2},m=\frac{N_b{m}_b}{2}$$

$$\frac{⟨L_{z,a}⟩}{N_a\hslash }=-i\int {\psi }_a^{\ast }\left(x\frac{\partial }{\partial y}-y\frac{\partial }{\partial x}\right){\psi }_a\mathrm{d}x\mathrm{d}y$$

$$\frac{{\tilde{L}}_{z,a}}{N_a\hslash }=2{\int }_{r_{vor}}^{R}2\pi r\frac{1}{\pi R^2}\mathrm{d}r=2\left[1-{\left(\frac{r_{vor}}{R}\right)}^2\right],$$

$$\frac{⟨L_{z,b}⟩}{N_b\hslash }=-i\int {\psi }_b^{\ast }\left(x\frac{\partial }{\partial y}-y\frac{\partial }{\partial x}\right){\psi }_b\mathrm{d}x\mathrm{d}y,$$

$${\vec{J}}_{b,rot}=m_b|{\psi }_b{|}^2{\vec{v}}_{b,rot},$$

$${\tilde{\Omega }}_{b,rot}(x,y)=-\frac{|{\vec{v}}_{b,rot}|}{r_C},$$

$$\frac{{\tilde{L}}_{z,b}}{N_b\hslash }\approx \frac{{\tilde{\mathcal{L}}}_{O,b}}{N_b\hslash }+\frac{{\tilde{\mathcal{S}}}_{C_L,b}}{N_b\hslash }+\frac{{\tilde{\mathcal{S}}}_{C_R,b}}{N_b\hslash },$$

$$\frac{{\tilde{\mathcal{L}}}_{O,b}}{N_b\hslash }=\frac{\Omega }{\hslash }\int m_b|{\psi }_b{|}^2{r}^2\mathrm{d}x\mathrm{d}y,$$

$$\frac{{\tilde{\mathcal{S}}}_{C_L,b}}{N_b\hslash }=-\frac{\Omega }{\hslash }\int m_b|{\psi }_b{|}^2\Theta (-x)r_{C_L}^2\mathrm{d}x\mathrm{d}y,$$

$$\frac{{\tilde{\mathcal{S}}}_{C_R,b}}{N_b\hslash }=-\frac{\Omega }{\hslash }\int m_b|{\psi }_b{|}^2\Theta (x)r_{C_R}^2\mathrm{d}x\mathrm{d}y$$

$${\vec{J}}_a={\rho }_{aa}{\vec{v}}_a+{\rho }_{ab}{\vec{v}}_b$$

$${\vec{J}}_b={\rho }_{ba}{\vec{v}}_a+{\rho }_{bb}{\vec{v}}_b$$

Full text

I Introduction

Vortices in quantum fluids are topological excitations characterized by quantized circulation Onsager (1949) which are present in a number of nonlinear field theories and models Pismen and Pismen (1999), ranging from superfluid media Feynman (1955); Donnelly (1991) and quantum optics Brambilla et al. (1991); Carusotto and Ciuti (2013) to superconductivity theories Blatter et al. (1994); Huebener (2013) and Josephson-junction arrays Vortex dynamics and quantum effects in Josephson junction arrays are reviewd in the Proceedings of the ICTP Workshop on Josephson Junction Arrays, edited by H.A. Cardeira and S.R. Shenoy Physica B 222 336 (1996); Fazio and Van Der Zant (2001) and play a key-role in fundamental effects such as superfluid turbulence Barenghi et al. (2001), the Berezinskii-Kosterlitz-Thouless transition Kosterlitz and Thouless (1973), fractional statistics Chiao et al. (1985), and in the development of a fully quantized field theory for topologically complex excitations Fetter (1967); Rasetti and Regge (1975); Penna et al. (1998). Among the plethora of different physical systems where vortices can be experimentally investigated, ultracold quantum gases provide a particularly controllable and versatile platform Fetter (2009); Dum et al. (1998) for the study and the observation of the rich phenomenology associated to their formation Butts and Rokhsar (1999); Dalfovo and Stringari (1996), dynamics Pismen and Rubinstein (1991); Cozzini and Stringari (2003), and interactions Castin and Dum (1999). Vortices in Bose-Einstein condensates (BECs) were first obtained by means of a phase-imprinting method involving two hyperfine spin states of ${}^{\text{87}}$Rb Matthews et al. (1999) but, at present, can be produced also by stirring the BEC above a certain critical velocity Recati et al. (2001); Sinha and Castin (2001); Madison et al. (2001), dragging barriers through the BEC itself Onofrio et al. (2000) or interfering multiple condensate fragments Scherer et al. (2007).

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 Drazin and Johnson (1989). Soon after the achievement of Bose-Einstein condensation, different types of solitons have been described and observed Dum et al. (1998); Burger et al. (2002); Frantzeskakis (2010). To our purposes, of particular importance are those systems where a bosonic binary mixture features dark-bright soliton configurations Becker et al. (2008); Hamner et al. (2011); Middelkamp et al. (2011). These structures, first predicted in Ref. Busch and Anglin (2001), are frequently termed as “symbiotic solitons” Pérez-García and Beitia (2005) 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 vortex-bright soliton complexes. First observed in 2000 by the JILA group Anderson et al. (2000), they represent the topological extension of the dark-bright soliton configuration to the case where a component hosts one or more vortices Law et al. (2010). 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 Mukherjee et al. (2019), the possibility, for the effective potential well corresponding to the vortex core, to support not only bound states Gallemí et al. (2018), but also multi-ring excited radial state complexes Charalampidis et al. (2016), and a rich dynamical scenario for the bright-solitary component Pola et al. (2012).

Within an analytical framework and by means of extensive numerical simulations, our work aims at analyzing the static and the dynamical properties of vortex-bright 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}>\sqrt{g_{a}g_{b}}$ is fulfilled, the dynamical picture of the mixture [in which the order-parameter fields of 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 Saffman (1992).

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 mean-field 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 species-b 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).

II.1 Massless vortices

The Hamiltonian of $N$ point-like massless vortices in an ideal unbounded fluid is given by Saffman (1992)

[TABLE]

where $\rho_{*}$ is the fluid planar density, $z_{j}=x_{j}+iy_{j}\,\in\,\mathbb{C}$ is the position of the $j$th vortex in the ambient plane and $k_{j}=n_{j}h/m_{f}$ is its strength ($n_{j}\,\in\,\mathbb{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. Henderson et al. (2009); Muldoon et al. (2012); Gauthier et al. (2016)), 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 V. Penna, in Quantized Vortex Dynamics and Superfluid Turbulence, edited by C.F. Barenghi, R.J.Donnelly, and W.F. Vinen , pp 445-452() (Springer-Verlag, Berlin 2001)

[TABLE]

[TABLE]

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

[TABLE]

involving, in turn, the canonical brackets $\{x_{i},y_{j}\}=\delta_{i,j}/(\rho_{*}k_{j})$ (see for example Penna (1999)).

II.2 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

[TABLE]

where $m_{j}$ represents the point-like mass hosted by the $j$th vortex core and where $\dot{q}_{j}:=\mathrm{d}q_{j}/\mathrm{d}t$ (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 $\rho_{*}$ 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

[TABLE]

[TABLE]

for $i,\,j\in\{1,2\}$ and $i\neq j$ [$\vec{u}_{3}$ is the unit vector perpendicular to the plane $(x,\,y)$].

The resulting system of 4 differential equations admits a notable solution,

[TABLE]

[TABLE]

provided that the two vortices are placed symmetrically with respect to the box-trap center and that their distance $d$ and the angular frequency $\Omega$ marking their motion of precession fulfill the following equation:

[TABLE]

[TABLE]

As expected, Eq. (4) shows a pathology when $d\to 2R$, meaning that the vortex pair is approaching the circular-box boundary. Moreover, in the limit of infinite box radius ($R\to+\infty$), one can retain only those terms $\propto R^{4}$ and the relation $d(\Omega)$ can be expressed in closed form, i.e.:

[TABLE]

In Sec. IV.1, 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}\mathrm{Na}$ and ${}^{39}\mathrm{K}$ Schulze et al. (2018); Richaud et al. (2019). Each atomic species is characterized by an order parameter, $\varphi_{a}=\sqrt{N_{a}}\psi_{a}$ and $\varphi_{b}=\sqrt{N_{b}}\psi_{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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $N_{a}$ ($N_{b}$) corresponds to the number of species-a (species-b) atoms, $g_{c}=4\pi\hbar^{2}a_{c}/m_{c}$, with $c=a,\,b$ are the intraspecies interaction strengths and $g_{ab}=2\pi\hbar^{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_{a}^{-1}+m_{b}^{-1})^{-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 $\ell_{z}$ is the effective thickness of the disk-like box trap and functions $\psi_{a}$ and $\psi_{b}$ are normalized to $1$, since $\varphi_{a}$ and $\varphi_{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 $\Omega$ (this corresponds to adding the term $-\Omega\hat{L}_{z}$ to the Hamiltonian, where $\hat{L}_{z}$ is the operator associated to the third component of the angular momentum) and then employing the imaginary-time method Coppens et al. (2017); Gallemí et al. (2018). 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}\,\in\,[5,\,1000]$), our numerical simulations Computational resources provided by HPC@POLITO () (http://www.hpc.polito.it) 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\,\mu\text{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 $|\psi_{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 Pola et al. (2012).

One can gain a further insight into the discussed eigensolution by computing the mass current density

[TABLE]

(with $c=a,\,b$) which also corresponds to the momentum-per-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 $|\psi_{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 $\psi_{a}$ and $\psi_{b}$ of equations (6) concerns the phase fields. The latter, denoted by $\theta_{a}$ and $\theta_{b}$, and such that $\psi_{c}=\sqrt{|\psi_{c}|^{2}}e^{i\theta_{c}}$ (where $c=a,\,b$), have been plotted in the lower panels of Fig. 1. As expected, $\theta_{a}$ features singularities in correspondence of the vortices’ centers, while, less expectedly, $\theta_{b}$ features singularities too. Nevertheless, the latter are found far from the peaks of $|\psi_{b}|^{2}$. According to the basic properties of quantum fluids, the circulation

[TABLE]

of the velocity vector field associated to $\psi_{c}$, that means $\vec{v}_{c}=\frac{\hbar}{m_{c}}\nabla\theta_{c}$ (with $c=a,\,b$), is zero if the closed path $\gamma$ does not encircle singularities of the phase field $\theta_{c}$ (with $c=a,\,b$). Conversely, due to the Feynman-Onsager quantization rules, it takes values $n\,h/m_{c}$, where $n\in\mathbb{Z}$, if one or more singularities of the associated field $\theta_{c}$ are encircled by $\gamma$. This information will be useful to better understand the discussion about the angular momentum of species-$b$ bosons (see Sec. V.2).

IV.1 Mass of the cores and equilibrium distance

Increasing the number of species-b atoms (within the investigated range $N_{b}\,\in\,[5,\,1000]$), the distance $d_{vor}$ between the centers of the vortices increases. Similarly, the distance between the two peaks of $|\psi_{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 $|\psi_{b}|^{2}$ does not exactly match that of the minima of $|\psi_{a}|^{2}$ due to the centrifugal force on species-b atoms and the finite repulsive coupling ($g_{ab}<+\infty$) 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

[TABLE]

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.

IV.2 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 species-b atoms, which therefore prove to be more tightly bound within the valleys of $|\psi_{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).

IV.3 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}/\sqrt{g_{a}g_{b}}>1$. Actually, we should mention that vortex/bright-soliton complexes prove to be rather robust composite objects, as their stability extends also beyond the immiscible regime. We have verified this by numerically solving Eqs. (6), thus extending the results of Ref. Hamner et al. (2011), derived in the context of dark-bright soliton complexes, to the case of vortex/bright-soliton complexes. The obtained density distributions, $|\psi_{a}|^{2}$ and $|\psi_{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 $\approx 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.

V.1 Angular momentum of condensate a

The angular momentum (per particle, in units of $\hbar$) of condensate a can be computed as

[TABLE]

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. Guilleumas and Graham (2001) (where the authors investigated the case of harmonic confinement), in fact, it is possible to derive the following expression:

[TABLE]

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.

V.2 Angular momentum of condensate b

The angular momentum (per particle, in units of $\hbar$) of component b can be analogously computed as

[TABLE]

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.

Mass current density in the rotating frame.

In the lab frame, it is possible to compute the mass current density $\vec{J}_{b}$ associated to $\psi_{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 $|\vec{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 $\Omega$ distinguishing the motion of precession of the vortices. More specifically, in the rotating frame, $\vec{J}_{b,rot}$ reads

[TABLE]

where $\vec{v}_{b,rot}=\vec{v}_{b}-\vec{V}$ with $\vec{v}_{b}=\frac{\vec{J}_{b}}{m_{b}|\psi_{b}|^{2}}$ and $\vec{V}=\vec{\Omega}\wedge\vec{r}$. $\vec{J}_{b}$ is numerically computed through Eq. (7). Notice also that, at $t=0$, $|\psi_{b}|^{2}\equiv|\psi_{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 $-\Omega$ (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 $\vec{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

[TABLE]

[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 $\approx-5$ rad/s in the most part of the regions where $|\psi_{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 $\vec{w}_{b,rot}=\vec{\nabla}\wedge\vec{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 quasi-plateaus at $|\vec{w}_{b,rot}|=2\tilde{\Omega}_{b,rot}\approx-10\,\mathrm{rad/s}$. This circumstance is in great agreement with the fact that the vorticity distribution associated to a rigid body rotating with angular frequency $\Omega$ is uniform and equal to $2\Omega$.

In conclusion, we have proved that, in the (non-inertial) rotating frame, the two species-b cores rotate around their respective centers of mass with (almost uniform) angular velocity $-\Omega$. This allows us to conclude that the they keep their orientation fixed when observed from the lab frame.

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

[TABLE]

where

[TABLE]

and where terms

[TABLE]

[TABLE]

are introduced to take into account that, in the lab frame, the two species-b cores revolve but keep their orientation fixed [Heaviside functions $\Theta(-x)$ and $\Theta(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_{C_{\alpha}}^{2}:=(x-x_{C_{\alpha}})^{2}+(y-y_{C_{\alpha}})^{2}$, with $\alpha=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\%$.

Observation of the associated phase field.

As illustrated in the lower right panel of Fig. 1, the phase field associated to $\psi_{b}$, $\theta_{b}$, features two singularities far from the centers of bright solitons. If one considers a closed path $\gamma$ encircling one of the bright solitons, since no singularities of $\theta_{b}$ are surrounded by it, the circulation $\mathcal{C}_{\gamma}[\vec{v}_{b}]$ of the velocity field $\vec{v}_{b}$ along $\gamma$ 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 $\theta_{b}$ indeed play an important role because the corresponding velocity field $\vec{v}_{b}=\frac{\hbar}{m_{b}}\nabla\theta_{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 soliton-like 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 Andreev and Bashkin (1976); Babaev (2004); Shevchenko and Fil (2007); Nespolo et al. (2017), according to which, the mass current density $\vec{J}_{i}$ (with $i=a,\,b$) in one species will depend, in general, also on the velocity $\vec{v}_{j}$ (with $j=b,\,a$) of the other species. In other words, the condensate density $\rho_{ij}$ is a non-diagonal matrix, a circumstance which implies the relations

[TABLE]

[TABLE]

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 Nespolo et al. (2017). 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 Andreev and Bashkin (1976). This effect is known to be rather elusive Nespolo et al. (2017) and, for our sample condition, we expect it to be negligible, as the off-diagonal matrix elements $\rho_{ab}$ and $\rho_{ba}$, which depend on the the overlap between $\psi_{a}$ and $\psi_{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 species-a 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 $\xi_{a}$ of species-a vortices and the equilibrium characteristic size of species-b cores, $\sigma_{b}$, we present the following heuristic equations

[TABLE]

[TABLE]

where

[TABLE]

Notice that the the first equation of system (18) reduces, in the case of no interspecies interaction, to

[TABLE]

the well-known formula derived in the context of single-species vortices Pitaevskii and Stringari (2016). Similarly, the second equation, if $g_{ab}=0$, gets structurally similar to formula

[TABLE]

giving the characteristic size of a soliton in the case of attractive interactions Pitaevskii and Stringari (2016) (of course, in this context, $n_{b,0}$ represents the central density). The extra term $g_{ab}n_{a}n_{b}\pi(\sigma_{b}^{2}-\xi_{a}^{2})\ell_{z}$ is introduced to take into account the interspecies repulsion, an interaction that manifests only in those regions where $\psi_{a}$ and $\psi_{b}$ overlap, i.e. only in the two annuli centered in the vortices’ centers and whose outer and inner radii are $\sigma_{b}$ and $\xi_{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 $|\psi_{a}|^{2}$, $\lambda_{a}$, and the half width at half maximum of the peak of $|\psi_{b}|^{2}$, $\lambda_{b}$. From the analytical side, the estimates of quantities $\lambda_{a}$ and $\lambda_{b}$ are given by the solutions of system (18), $\xi_{a}$ and $\sigma_{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 $\lambda_{a}$ and of $\lambda_{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}<\sqrt{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 Lorentz-like 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.

Acknowledgements

The authors wish to thank M. Pi, A. Gallemí and Y. Bidasyuk for their valuable help in setting up the numerical simulations. V. P. is indebted with F. Minardi for useful discussions about trapping techniques. M. G. and R. M. acknowledge financial support from the Spanish MINECO and Fondo Europeo de Desarrollo Regional (FEDER, EU) under Grant No. FIS2017-87801-P, and from Generalitat de Catalunya Grant No. 2017SGR533.