Spin and mass currents near a moving magnetic obstacle in a two-component Bose-Einstein condensate

We first consider a perturbative regime in which the obstacle moves slowly such that the densities of the spin components are well approximated by the solutions of Eq. (5) for $u=0$, i.e., $n_{i}(r)=n_{0}+s_{i}V(r)/\Delta g$. In a later discussion, it will be clear that the approximation is valid when the kinetic energy of the induced flow is negligible compared to the spin interaction energy, i.e., $mu^{2}\ll mc_{s}^{2}\equiv\Delta gn_{0}$. Here, $c_{s}$ is the speed of spin sound for the unperturbed BEC.

For the density distribution $n_{i}(r)$, the potential function $S_{i}(\bm{r})$ can be directly determined using Eq. (6). Because $n_{i}$ has only an $r$-dependence, we perform separation of variables, i.e., ${S_{i}(r,\phi)}={R_{i}(r)}{\Phi_{i}(\phi)}$, and Eq. (6) is transformed to

with $\delta\tilde{n}_{i}(r)=n_{i}(r)/n_{0}-1$ and $l$ being an integer. The boundary condition of ${S}_{i}\rightarrow-ur\cos{\phi}$ as $r\rightarrow\infty$ imposes ${l}=1$, and without loss of generality, we set $\Phi_{i}(\phi)=-u\cos\phi$. The solution for the radial function, $R_{i}(r)$, can be obtained perturbatively using the small parameter $\alpha=|\delta\tilde{n}_{i}(0)|=V_{0}/(\Delta gn_{0})\ll 1$, which is the maximum magnitude of the relative density variations in each spin component. When the radial function is expanded in a power series with respect to $\alpha$ as $R_{i}(r)=\sum_{k\geq 0}\alpha^{k}R_{i}^{(k)}(r)$, the $k$-th function $R_{i}^{(k)}(r)$ is recursively determined from Eq. (9) as the solution to the following equation:

Up to the second order of $\alpha$, we have

where $\rho=r/\sigma$, thereby yielding the approximate solution of $S_{i}(\bm{r})$ as

which satisfies the boundary condition as $r\rightarrow\infty$.

From Eqs. (7) and (8) with $n_{i}(r)=n_{0}(1+s_{i}\alpha e^{-2\rho^{2}})$ and $S_{i}(r,\phi)$ in Eq. (12), we obtained the analytic expressions of the spin and the mass currents as

respectively. Of note is that $|\bm{J}|\propto\alpha$ and $|\bm{M}|\propto\alpha^{2}$, which indicate that the moving magnetic obstacle dominantly generates a spin current, as expected, as well as a mass current via a nonlinear process. The peak spin and mass currents occur at the obstacle center, and are given by $\bm{J}_{0}=\alpha n_{0}\bm{u}$ and $\bm{M}_{0}=\frac{\ln 2}{2}\alpha^{2}n_{0}\bm{u}$, respectively.

Figures 2(a) and (b) show the spin and the mass current distributions predicted using Eqs. (13) and (14), respectively. We used $g_{\uparrow\downarrow}/g=0.93$, which is the value for a mixture of ²³Na in the $|F=1,m_{F}=\pm 1\rangle$ states R17 ; Knoop11 . As expected, both spin and mass currents appeared locally near the moving obstacle. In the spin current distribution, we observed two low-current holes, one in the front and the other in the back of the obstacle; furthermore, the spin current flowed out from the back hole and into the front one. Meanwhile, in the mass current distribution, we observed two low-current holes on the lateral sides of the obstacle, and the mass current swirled around each hole in opposite directions. The flow patterns of $\bm{J}$ and $\bm{M}$ around the moving obstacle resembled those of an electric field around a charge dipole and a magnetic field around a current loop, respectively. In Figs. 2(c) and 2(d), we present the distributions of $\nabla\cdot\bm{J}$ and $\nabla\times\bm{M}$, respectively, which clearly show the dipole configurations of the source and the sink for the spin current and the vorticity of the mass current, respectively.

To understand the characteristic flow patterns of the spin and the mass currents, we analyzed the general divergence and rotation properties of $\bm{J}$ and $\bm{M}$. From the continuity equation in Eq. (3), we obtained $\nabla\cdot(n_{i}\bm{u}_{i})=\bm{u}\cdot\nabla n_{i}$. Combining the latter with the irrotational property of $\nabla\times\bm{u}_{i}=0$, we obtain the following relations:

The first and the second equations result from spin and mass conservation, respectively, and the third and the fourth ones reveal intriguing nonlinear coupling between the spin and the mass channels in the binary system.

For a weak and slow magnetic obstacle, taking the same level of approximation as in the previous subsection, we have $\nabla n_{t}=0$, $\bm{u}_{\uparrow}+\bm{u}_{\downarrow}=\mathcal{O}(\alpha^{2})$, and $\bm{u}_{\uparrow}-\bm{u}_{\downarrow}=\frac{\bm{J}}{n_{0}}+\mathcal{O}(\alpha^{3})$. Then, up to the first order in $\alpha$, the relations can be expressed as

which immediately explains the observed electric- and magnetic-field-like behaviors of $\bm{J}$ and $\bm{M}$, respectively, near the moving magnetic obstacle. The quantities $Q_{J}$ and $\bm{I}_{M}$ can be regarded as the ‘charge’ and the ‘current’ source densities for generating the spin and the mass currents, respectively. Their expressions are consistent with the previous results of $|\bm{J}|\propto\alpha$ and $|\bm{M}|\propto\alpha^{2}$ for $m_{z}\propto\alpha$. Furthermore, we may consider the ‘electric’ and the ‘magnetic’ dipole moments as

respectively, which allow us to predict the currents in the far distant region of $r\gg\sigma$ to be $\bm{J}\sim\frac{2(\bm{p}_{J}\cdot\hat{r})\hat{r}-\bm{p}_{J}}{2\pi r^{2}}$ and $\bm{M}\sim\frac{2(\bm{\mu}_{M}\cdot\hat{r})\hat{r}-\bm{\mu}_{M}}{2\pi r^{2}}$. We emphasize that the relations in Eq. (16) hold regardless of the potential form of $V(\bm{r})$ once the magnetic obstacle is weak and slow.

The existence of nonzero $\nabla\times\bm{M}$ should be highlighted. The superfluid velocity of the binary superfluid system can be expressed as ${\bm{u}}_{M}=\frac{\bm{M}}{n_{t}}$; therefore, the circulation of ${\bm{u}_{M}}$ can be nonzero for $\nabla\times\bm{M}\neq 0$. This is in stark contrast to the case with a single-component BEC, where the circulation of the superfluid velocity should be quantized with $h/m$ as a topological invariant of the system. Noting that the mass circulation of the binary superfluid system can have a continuous value in conjunction with the spin current is important.

To investigate how the flow patterns evolve with increasing obstacle velocity, we numerically calculated $\{n_{i}(\bm{r}),S_{i}(\bm{r})\}$ from Eqs. (5) and (6) for various $u$. A finite difference method was employed for a $241\times 241$ grid system, and the obstacle width, $\sigma$, was set to 40 grid spacings. The boundary conditions imposed were $n_{i}=n_{0}$ and $\bm{v}_{i}=-u\hat{\bm{x}}$ on the edge of the grid system.

Figures 3(a)–(d) show the numerical results for the spin and the mass currents for $u=0.1c_{s}$ and $0.9c_{s}$ with $V_{0}=0.1mc_{s}^{2}$. The numerical results for low $u$ were confirmed to be in good quantitative agreement with the analytical predictions based on Eqs. (13) and (14). We observe that as $u$ increases, the spatial distributions of $\bm{J}$ and $\bm{M}$ stretch along the lateral and the moving directions, respectively, but the flow patterns maintain their characteristic spatial structures [Figs. 3(c) and (d)]. The peak currents still occur at the center of the moving obstacle. In Figs. 3(e) and (f), we plot $|\bm{J}_{0}|$ and $|\bm{M}_{0}|/|\bm{J}_{0}|$ as functions of $u$, respectively. For low $u$, $|\bm{J}_{0}|$ increases linearly with $u$, as predicted in Eq. (13); however, it begins deviating upwardly as $u$ increases over $\approx 0.6c_{s}$. The ratio $|\bm{M}_{0}|/|\bm{J}_{0}|$ increases quadratically with increasing $u$, departing from the predicted value of $\frac{\ln 2}{2}\alpha$.

The nonlinear $u$-dependence of $\bm{J}_{0}$ can be attributed to the additional density variations due to the increased flow velocity for high $u$. When the first-order kinetic energy correction term related to $u^{2}$ in Eq. (5) is considered, the magnetization can be expressed as

Because $\bm{u}_{rel}=\bm{u}_{\uparrow}-\bm{u}_{\downarrow}\approx\frac{\bm{J}}{n_{0}}$, the magnetization is enhanced in the center region where the spin current flows in the direction of the obstacle’s motion. As the first relation of Eq. (16) shows, this enhancement in $m_{z}$ results in an increase in the spin current. This mutual enhancing effect qualitatively explains the observed lateral stretching of the elongated, high-$|\bm{J}|$ region with high $u$.

If the magnetization distribution maintains its Gaussian form for high $u$, i.e., $m_{z}(r)=m_{z,0}e^{-\frac{2r^{2}}{\sigma^{2}}}$, we can infer ${\bm{J}_{0}}=\frac{m_{z,0}}{2}\bm{u}$ from Eq. (13) because of the relation $\alpha=\frac{m_{z,0}}{2n_{0}}$. Substituting $\bm{u}_{rel}\approx\frac{\bm{J}_{0}}{n_{0}}=\frac{m_{z,0}}{2n_{0}}\bm{u}$ into Eq. (17), we obtain the magnetization at the obstacle’s center as $m_{z,0}=\frac{2n_{0}V_{0}}{\Delta gn_{0}-\frac{1}{2}mu^{2}}$. This suggests that the high-$u$ effect in the spin and the mass currents might be captured by replacing $\alpha$ in Eqs. (13) and (14) with its effective value, i.e.,

In fact, we observe that the numerical results for $|\bm{J}_{0}|$ and $|\bm{M}_{0}|/|\bm{J}_{0}|$ can be explained well quantitatively with $\alpha_{\text{eff}}$, i.e., $|\bm{J}_{0}|=\alpha_{\text{eff}}n_{0}u$ and $|\bm{M}_{0}|/|\bm{J}_{0}|=\frac{\ln 2}{2}\alpha_{\text{eff}}$, respectively [Figs. 3(e) and (f)].

When the obstacle’s velocity increases above a certain critical value, energy dissipation will occur in the binary superfluid system. According to the Landau criterion, the critical velocity is expressed as $u_{L}=\min[\varepsilon(\bm{p})/(\bm{p}\cdot\hat{\bm{u}})]$ R1 , where $\varepsilon(\bm{p})$ is the elementary excitation energy of momentum $\bm{p}$, and $\hat{\bm{u}}$ is the unit vector along the direction of the obstacle’s motion. In general, in the long wavelength limit, the superfluid system has a linear dispersion of $\varepsilon(p)=cp$ with $c$ being the speed of sound, and the Landau critical velocity is given as $u_{L}=c$. In this section, we investigate the critical velocity $u_{c}$ of the magnetic obstacle based on the local Landau criterion, i.e., by comparing the obstacle’s velocity to the local speed of sound.

First, we determine the speed of sound for a stationary state, in which the two spin components flow with uniform velocity $\bm{u}_{i}$, having unifrom density $n_{i}$. Linearizing the hydrodynamic equations, Eqs. (1) and (2), with $V_{i}=0$ R28 , we obtain

Furthermore, the coupled wave equations for $\delta n_{\uparrow}$ and $\delta n_{\downarrow}$ can be obtained as follows:

If a traveling wave solution of $\delta n_{i}=A_{i}e^{i(\bm{q}\cdot\bm{r}-\omega t)}$ is to be obtained, the wave velocity $c=\omega/q$ should satisfy

with $\hat{\bm{q}}=\frac{\bm{q}}{|\bm{q}|}$. In general, four solutions for $c$ are provided, but because $c(-\hat{\bm{q}})=-c(\hat{\bm{q}})$, we consider only the two positive solutions for the propagation direction of $\hat{\bm{q}}$ and denote them by $c^{+}$ and $c^{-}$ with $c^{+}\geq c^{-}$. For $n_{i}=n_{0}$ and $\bm{u}_{i}=0$, the fast (slow) sound speed is given by $c^{\pm}=c_{n(s)}=\sqrt{\frac{(g\pm g_{\uparrow\downarrow})n_{0}}{m}}$, and the sound propagates with $A_{\uparrow}=A_{\downarrow}$ ($A_{\uparrow}=-A_{\downarrow}$), corresponding to phonon (magnon) excitations in a symmetric binary superfluid system.

Figure 4 shows the sound speeds $c^{\pm}(\hat{\bm{q}})$ for various flow conditions of $\bm{u}_{\uparrow}=\frac{u_{\text{rel}}}{2}\hat{\bm{x}}$, $\bm{u}_{\downarrow}=-\frac{u_{\text{rel}}}{2}\hat{\bm{x}}$, $n_{\uparrow}=n_{0}+\frac{m_{z}}{2}$, and $n_{\downarrow}=n_{0}-\frac{m_{z}}{2}$. We observe that $c^{+}$ is not significantly affected by changes in $u_{\text{rel}}$ and $m_{z}$, implying the strong phonon characteristics of the fast sound. Meanwhile, $c^{-}$ is sensitive to them. This decreases with increasing $u_{\text{rel}}$ and $m_{z}$, and the reduction rate is the fastest along the spin current direction. In our moving-obstacle situation, the relative velocity and the density imbalance between the two spin components are maximum at the obstacle’s center. Additionally, the obstacle’s direction of motion is the same as the direction of the spin current. Therefore, as the obstacle’s velocity increases, the stability of the induced superflow will break first in the region of the obstacle’s center according to the local Landau criterion.

Figure 5(a) shows the speed of sound, $c^{-}(\hat{\bm{x}})$, at the obstacle’s center as a function of $u$ for various $\alpha$ from $0.1$ to $0.9$. The local flow condition of $\{n_{i},\bm{u}_{i}\}$ in the center region was numerically obtained for each set of $\{u,\alpha\}$, and $c^{-}(\hat{\bm{x}})$ was determined from Eq. (22). As $u$ increases, $c^{-}$ decreases and eventually becomes equivalent to $u$, which defines the obstacle’s cirtical velocity, $u_{c}$. Note that the countersuperflow instability corresponds to an imaginary solution of $c$ in Eq. (22) and is irrelevant to our current study.

The critical velocity $u_{c}$ decreases with increasing potential magnitude, $V_{0}$ [Fig. 5(b)], which is attributable to the reduction in $c^{-}$ due to the increased $m_{z,0}$. In the limit of $V_{0}\rightarrow 0$, $u_{c}$ approaches $c_{s}$ as $m_{z,0}\rightarrow 0$ and $\bm{u}_{i}\rightarrow 0$. When $V_{0}$ reaches $mc_{s}^{2}$, i.e., $\alpha=1$, $u_{c}$ vanishes because the system is fully polarized at the obstacle’s center and $c^{-}=0$. Interestingly, our numerical results indicate that $u_{c}$ decreases almost linearly with increasing $V_{0}$, suggesting an empirical critical line of $u_{c}(V_{0})=c_{s}(1-\frac{V_{0}}{mc_{s}^{2}})$. We also scrutinized how $u_{c}$ was affected by the intercomponent interaction strength and observed that when $\Delta{g}/g$ increased from the ²³Na value of 0.07, $u_{c}$ decreased for the same obstacle condition [Fig. 5(c)]. Because $\Delta{g}/g=1$ corresponds to a non-interacting two-component case, we may conclude that the observed linear dependence of $u_{c}$ on $V_{0}$ is driven by the interactions between the two spin components.