---
abstract: 'We study the hydrodynamic forces acting on a finite-size impurity moving in a two-dimensional Bose-Einstein condensate at non-zero temperature. The condensate is modeled by the damped-Gross Pitaevskii (dGPE) equation and the impurity by a Gaussian repulsive potential giving the coupling to the condensate. The width of the Gaussian potential is equal to the coherence length, thus the impurity can only emit waves. Using linear perturbation analysis, we obtain analytical expressions corresponding to different hydrodynamic regimes which are then compared with direct numerical simulations of the dGPE equation and with the corresponding expressions for classical forces. For a non-steady flow, the impurity experiences a time-dependent force that, for small coupling, is dominated by the inertial effects from the condensate and can be expressed in terms of the local material derivative of the fluid velocity, in direct correspondence with the Maxey and Riley theory for the motion of a solid particle in a classical fluid. In the steady-state regime, the force is dominated by a self-induced drag. Unlike at zero temperature, where the drag force vanishes below a critical velocity, at finite temperatures, the drag force has a net contribution from the energy dissipated in the condensate through the thermal drag at all velocities of the impurity. At low velocities this term is similar to the Stokes’ drag in classical fluids. There is still a critical velocity above which the main drag pertains to energy dissipation by acoustic emissions. Above this speed, the drag behaves non-monotonically with impurity speed, reflecting the reorganization of fringes and wake around the particle.'
author:
- 'Jonas Rønning$^1$, Audun Skaugen$^2$, Emilio Hernández-García$^3$, Cristóbal López$^3$, Luiza Angheluta$^1$'
bibliography:
- 'ref-2.bib'
title: 'Classical analogies for the force acting on an impurity in a Bose-Einstein condensate'
---

Introduction
============

The motion of an impurity suspended in a quantum fluid depends on several key factors such as the superfluid nature and flow regime, as well as the size of the impurity and its interaction with the surrounding fluid [@winiecki2000motion; @wouters2010superfluidity; @astrakharchik2004motion; @shukla2016sticking; @pinsker2017gaussian]. Therefore, it is disputable whether the forces acting on an impurity in a quantum fluid should bear any resemblance to classical hydrodynamic forces. In the case of an impurity immersed in superfluid liquid helium, classical equations of motion and hydrodynamic forces are assumed a priori [@poole2005motion], since impurities are typically much larger than the coherence length and then quantum hydrodynamic effects like the quantum pressure can be neglected. For Bose-Einstein condensates (BEC) in dilute atomic gases, impurities can be neutral atoms [@chikkatur2000suppression], ion impurities [@zipkes2010trapped; @balewski2013coupling] or quasiparticles [@jorgensen2016observation]. The size of an impurity in a BEC is typically of the same order of magnitude or smaller than the coherence length, and quantum hydrodynamic effects cannot be ignored.

There are several theoretical and computational studies of the interaction force between an impurity and a BEC at zero absolute temperature, using different approaches depending on the nature of the particle and its interaction with the condensate. A microscopic approach is used to analyse the interaction of a rigid particle with a BEC by solving the Gross-Pitaevskii equation (GPE) for the condensate macroscopic wavefunction and using boundary conditions such that the condensate density vanishes at the particle boundary [@pham2005boundary]. This methodology allows to study complex phenomena such as vortex nucleation and flow instabilities, but it is more oriented to find the effects of an obstacle on the flow rather than the coupled particle-flow dynamics. In addition, the boundary condition introduces severe nonlinearities which can only be addressed numerically. At a more fundamental level of description, the impurity is treated as a quantum particle with its own wavefunction described by the Schrödinger equation and that is coupled with the GPE for the macroscopic wavefunction of the BEC [@berloff2000capture]. A more versatile model for the interaction of impurities with the BEC has been explored in several papers [@astrakharchik2004motion; @shukla2016sticking; @griffin2017vortex; @shukla2018particles; @pinsker2017gaussian]. Here, an additional repulsive interaction (a Gaussian or delta-function potential) is added to model scattering of the condensate particles with the impurity. The hydrodynamic force on the impurity is determined by this repulsion potential and the superfluid density through the Ehrenfest theorem. The strong-coupling limit of this repulsive potential would be equivalent to the rigid boundary-condition approach. Within this modeling approach, some works have studied the complex motion of particles interacting with vortices in the flow, and the indirect interactions between them arising from the presence of the fluid [@shukla2016sticking; @shukla2018particles]. Another line of research using this type of modeling focused mainly on the superfluidity criterion of a uniform BEC at zero temperature regime. Within the Bogoliubov perturbation analysis for a small impurity and weak interaction, analytical expressions can be derived for the steady-state force exerted by the superfluid as function of the constant velocity of the impurity [@astrakharchik2004motion; @roberts2006force; @sykes2009drag; @pinsker2017gaussian]. At zero temperature, this force vanishes below a critical velocity, the speed of long-wavelength sound waves, at least when we ignore the quantum fluctuations [@roberts2006force], and corresponds to the dissipationless motion. Above this velocity, there is a finite drag force and the motion of the impurity is damped by acoustic excitations. While this is a form of drag, in that the force opposes motion by dissipating energy, it is not the same as the classical Stokes’ drag in viscous fluids. Recent experiments probing superfluidity in a BEC are able to indirectly estimate the drag force by measuring the local heating rate in the vicinity of the moving laser beam and show that there is still a critical velocity even at non-zero temperatures and that the critical velocity is lower for a repulsive potential than for an attractive one [@singh2016probing].

In this paper, we study the forces exerted on an impurity moving in a two-dimensional BEC at finite temperature, using an approach similar to [@astrakharchik2004motion; @shukla2016sticking; @griffin2017vortex; @shukla2018particles; @pinsker2017gaussian], in which a repulsive Gaussian potential is used to describe the interaction of the particle with the BEC, but using a dissipative version of the GPE to model the fluid. Our aim is to bridge this microscopic approach with the phenomenological descriptions [@poole2005motion] that assume that the forces from the superfluid are the same as those from a classical fluid in the inviscid and irrotational case. As in the classical-fluid case, we find that the force is made of two contributions: One of them, dominant for very weak fluid-particle interaction, bears a rather complete analogy with the corresponding force in classical fluids (inertial or pressure-gradient force), which depends on local fluid acceleration and includes the so-called Faxén corrections arising from velocity inhomogeneities close to the particle position [@maxey1983equation]. The difference is that, in a classical fluid, these corrections arise from the finite size of the particle and vanish when the particle size becomes zero. In the BEC, Faxén-type corrections arise both from the particle size (modeled by the range of the particle repulsion potential) and from the BEC coherence length. As fluid-particle interaction becomes more important, a second contribution to the force becomes noticeable, which takes into account the drag on the particle arising from the perturbation of the flow produced by the presence of the particle. Thus it can be called a particle *self-induced* force. We are able to obtain explicit formulae for it in the case of constant-velocity motion of the particle in an otherwise homogeneous and steady BEC. This drag is a dissipative (damping) force due to viscous-like drag of the perturbed BEC with the thermal cloud. It occurs in addition to the drag due to acoustic excitations in the condensate that in the absence of dissipation occurs only above a critical velocity for the particle. Here, it can be compared with the corresponding force in classical fluids, namely the viscous Stokes drag. We find that, as the Stokes force, the self-induced dissipative drag is linear in the particle velocity for small velocities, and we obtain an expression for it also at arbitrary velocities.

The rest of the paper is structured as follows. In Sect. \[sec:model\], we discuss the general modeling setup and in Sect. \[sec:perturbation\] a perturbation analysis is used to derive the linearized equations for the perturbations in the wavefunction related to non-steady condensate flow and the particle repulsive potential. Subsections \[sec:inertial\] and \[sec:drag\] derive analytical expressions within perturbation theory for the two contributions to the force experienced by the particle. In Section \[sec:numerics\], we compare our theoretical predictions with numerical simulations of the dissipative GPE coupled to the impurity, and the final section summarizes our conclusions.

Modeling approach {#sec:model}
=================

We model the interaction between the impurity and a two-dimensional BEC through a Gaussian repulsive potential which can be reduced to a delta-function limit similar to previous studies [@astrakharchik2004motion; @pinsker2017gaussian]. The BEC itself, which is at a finite temperature, is described by a macroscopic wavefunction $\psi({\boldsymbol{r}},t)$ that evolves according to the damped Gross Pitaevskii equation (dGPE) [@Bradley_2012; @Reeves_2013; @skaugen2016vortex]: $$\begin{aligned}
\label{eq:GPe}
&& i\hbar\partial_t\psi = \nonumber \\
&& (1-i\gamma)\left(-\frac{\hbar^2}{2m}\nabla^2+g|\psi|^2-\mu+ V_{ext}+g_p{\cal U}_p\right)\psi,\end{aligned}$$ where $g$ is an effective scattering parameter between condensate atoms. $V_{ext}$ is any external potential used to confine the condensate atoms or to stir them. The damping coefficient $\gamma>0$ is related to finite-temperature effects due to thermal drag between the condensate and the stationary thermal reservoir of excited atoms at fixed chemical potential $\mu$. This damping $\gamma$ is very small at low temperatures and can be expressed as function of temperature $T$, chemical potential $\mu$ and the energy of the thermal cloud [@bradley2008bose]. The dGPE can be derived from the stochastic projected Gross-Pitaevskii equation in the low-temperature regime [@blakie2008dynamics] and has been used to study different quantum turbulence regimes and vortex dynamics [@Reeves_2014; @billam2015spectral; @Reeves_2013; @Bradley_2012] that are also observed in recent experiments [@neely2013characteristics]. An hydrodynamic description in terms of density and velocity of the BEC can be developed using the Madelung transformation of the wavefunction: $\psi=|\psi|e^{i\phi}$. The macroscopic number density is $\rho({\boldsymbol{r}},t)=|\psi({\boldsymbol{r}},t)|^2$ and the condensate velocity is ${\boldsymbol{v}}({\boldsymbol{r}},t)=(\hbar/m) \nabla
\phi({\boldsymbol{r}},t)$. This velocity can also be obtained from the superfluid current ${\boldsymbol{J}}({\boldsymbol{r}},t)$ as $${\boldsymbol{J}} = \frac{\hbar}{2 m i} \left(\psi^*\nabla\psi - \psi\nabla\psi^* \right) =\rho {\boldsymbol{v}}  \ .
\label{eq:current}$$ $\psi^*$ denotes the complex conjugate of $\psi$. In addition to damping the BEC velocity, the presence of $\gamma\neq 0$ in the dGPE also singles out the value $\rho_h=|\psi|^2=g/\mu$ as the steady homogeneous density value when the phase is constant and $V_{ext}=0$.

The interaction potential ${\cal U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p)$ between the condensate and the impurity is modeled by a Gaussian potential ${\cal U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p)=
\mu/(2\pi \sigma^2) e^{-({\boldsymbol{r}}-{\boldsymbol{r}}_p)^2/(2\sigma^2)}$. The parameter $g_p>0$ is the weak coupling constant for repulsive impurity-condensate interaction, ${\boldsymbol{r}}_p={\boldsymbol{r}}_p(t)$ denotes the center-of-mass position of the impurity, and $\sigma$ its effective size. Here we consider an impurity of size $\sigma$ of the order the coherence length $\xi=\hbar/\sqrt{m\mu}$ of the condensate. The impurity is too small to nucleate vortices in its wake [@reeves2015identifying]. Instead, the impurity will create acoustic excitations which at supersonic speeds correspond to shock waves. Similar acoustic fringes in the condensate density have been reported numerically in [@wouters2010superfluidity] for a different realization of non-equilibrium Bose-Einstein condensates. In the limit of a point-like impurity, the Gaussian interaction potential converges to a two-body scattering potential ${\cal U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p,t)= \mu \delta({\boldsymbol{r}}-{\boldsymbol{r}}_p(t))$ that has been used in previous analytical studies [@astrakharchik2004motion; @shukla2016sticking; @griffin2017vortex; @shukla2018particles]. Note that we are modeling only the interaction of the particle with the BEC, so that the viscous-like drag we will obtain arises from the indirect coupling to the thermal bath via the BEC. Any direct interaction of the particle with the thermal cloud of normal atoms will lead to additional forces that we do not consider here.

In order to gain insight into the forces and their relationship with the classical case, we keep the set-up as simple as possible. We consider a strictly two-dimensional condensate, which can be obtained by plane optical traps. We assume that the size of the condensate is large enough so that we can neglect inhomogeneities in the confining part of $V_{ext}$ in the region of interest. Also, we consider a neutrally buoyant impurity so that effects of gravity can be neglected. This would imply $V_{ext}=0$ except if an external forcing is introduced to stir the system, in which case we assume the support of this external force is sufficiently far from the impurity.

The impurity and the condensate will exert an interaction force on each other that is determined by the Ehrenfest theorem for the evolution of the center-of-mass momentum of the particle. The potential force $-g_p\nabla {\cal U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p)$ is the force exerted by an impurity on a condensate particle at position ${\boldsymbol{r}}$. By space averaging over condensate density, we then determine the force exerted by the impurity on the condensate as $-g_p\int d{\boldsymbol{r}} |\psi({\boldsymbol{r}},t)|^2 \nabla {\cal
U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p)$  [@shukla2018particles]. Hence, the force acting on the impurity has the opposite sign and is equal to $$\begin{aligned}
\label{eq:fp0}
{\boldsymbol{F}}_p(t) &=& + g_p\int d^2{\boldsymbol{r}} |\psi({\boldsymbol{r}},t)|^2 \nabla {\cal U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p)\end{aligned}$$ which, through an integration by parts, is equivalent to $$\begin{aligned}
\label{eq:fp1}
{\boldsymbol{F}}_p(t) &=& -g_p\int d^2{\boldsymbol{r}}  ~{\cal U}_p({\boldsymbol{r}}-{\boldsymbol{r}}_p ,t)  \nabla |\psi({\boldsymbol{r}},t)|^2 .
\label{eq:forcegradrho}\end{aligned}$$ Note that this last expression can also be used, reversing the sign, to give the force exserted on the BEC by a laser of beam profile given by ${\cal U}_p$.

At zero temperature, i.e. $\gamma=0$, and neglecting the effect of quantum fluctuations [@roberts2006force; @sykes2009drag], the impurity moves without any drag through a uniform condensate below a critical velocity, which is the low-wavelength speed of sound $c=\sqrt{\mu/m}$, as determined by the condensate linear excitation spectrum, in agreement with Landau’s criterion of superfluidity [@astrakharchik2004motion]. Above the critical speed, the impurity will create excitations, and depending on the size of the impurity these excitations range from acoustic waves (Bogoliubov excitation spectrum) to vortex dipoles and to von-Karman street of vortex pairs [@reeves2015identifying]. Previous studies focused on the theoretical investigations of the self-induced drag force and energy dissipation rate in the presence of Bogoliubov excitations emitted by a pointwise  [@astrakharchik2004motion; @roberts2006force; @sykes2009drag] or finite-size [@pinsker2017gaussian] particle, or numerical investigations of the drag force due to vortex emissions [@winiecki2000motion; @griffin2017vortex; @shukla2018particles]. The energy dissipation rate depends on whether the impurity is heavier, neutral or lighter with respect to the mass of the condensate particles [@shukla2018particles]. The dependence on the velocity of the self-induced drag force above the critical velocity changes with the spatial dimensions [@astrakharchik2004motion]. This means that the energy dissipation rate is also dependent on the spatial dimensions. If instead of a single impurity one considers many of them there will be, besides direct inter-particle interactions, additional forces between the impurities mediated by the flow, leading to a much more complex many-body dynamics even in an otherwise uniform condensate, as discussed in [@shukla2016sticking]. Here we neglect all these effects and consider a single impurity in a two-dimensional BEC.

We rewrite the dGPE in dimensionless units by using the characteristic units of space and time in terms of the long-wavelength speed of sound $c=\sqrt{\mu/m}$ in the homogeneous condensate and the coherence length $\xi=\hbar/(m
c)=\hbar/\sqrt{m\mu}$. Space is rescaled as ${\boldsymbol{r}}\rightarrow
\tilde{{\boldsymbol{r}}} \xi$ and time as $t\rightarrow \tilde t \xi/c$. In addition, the wavefunction is also rescaled $\psi
\rightarrow \tilde \psi\sqrt{\mu/g}$, where $g/\mu$ is the equilibrium particle-number density corresponding to the solution with constant phase if $V_{ext},\mathcal U_p=0$. The external potential, $V_{ext}=\mu \tilde V_{ext}$, and the interaction potential, $g_p \mathcal U_p = \mu \tilde g_p
\tilde{\mathcal U_p}$, are measured in units of the chemical potential $\mu$ with $\tilde{\mathcal U_p} = 1/(2\pi a^2)
e^{-({{\boldsymbol{\tilde}} r}-{{\boldsymbol{\tilde}} r}_p)^2/(2a^2)}$, and $a =
\sigma/\xi$, $\tilde{g}_p = g_p/(\xi^2\mu)$. Henceforth, the dimensionless form of the dGPE reads as $$\begin{aligned}
\label{eq:GPe_dimless}
\tilde \partial_t \tilde\psi =
(i+\gamma)\left(\frac{1}{2}\tilde\nabla^2+1-\tilde V_{ext}-\tilde g_p\tilde{\mathcal U_p}-|\tilde\psi|^2\right)\tilde\psi  .\end{aligned}$$ We use these dimensionless units and express the force (\[eq:forcegradrho\]) exerted on an impurity as ${\boldsymbol{F}}_p=
(\mu^2\xi/g) \tilde{{\boldsymbol{F}}}_p$, where $$\begin{aligned}
\tilde{{\boldsymbol{F}}}_p(t) =
-\tilde g_p\int d^2\tilde{{\boldsymbol{r}}} \tilde{\mathcal U_p}({\boldsymbol{r}}-{\boldsymbol{r}}_p)\tilde\nabla |\tilde\psi(\tilde{{\boldsymbol{r}}},t)|^2. \label{eq:fp2}\end{aligned}$$ For the rest of the paper, we will now omit the tildes over the dimensionless quantities.

In the limit of a point-like particle, $\mathcal U_p=
\delta({\boldsymbol{r}}-{\boldsymbol{r}}_p)$, the force from Eq. (\[eq:fp2\]) becomes $$\begin{aligned}
{\boldsymbol{F}}_p(t) = -g_p\nabla|\psi({\boldsymbol{r}},t)|^2 |_{{\boldsymbol{r}}={\boldsymbol{r}}_p(t)}.
\label{eq:fp3}\end{aligned}$$

Perturbation analysis {#sec:perturbation}
=====================

For a weakly-interacting impurity, the condensate wavefunction $\psi$ can be decomposed into an unperturbed wavefunction $\psi_0({\boldsymbol{r}})$ describing the motion and density of the fluid in the absence of the particle and the perturbation $\delta\psi_1({\boldsymbol{r}})$ due to the impurity’s repulsive interaction with the condensate, hence $\psi = \psi_0+g_p
\delta\psi_1$. The unperturbed wavefunction $\psi_0({\boldsymbol{r}},t)$ can be spatially-dependent, if it is initialized in a nonequilibrium configuration, or if external forces characterized by $V_{ext}$ are at play. Here, we consider deviations with respect to the steady and uniform equilibrium state (which in our dimensionless units is $\psi_h=1$). As stated before, we do not consider large extended inhomogeneities produced by a trapping potential, and assume that any stirring force acting on the BEC is far from the particle. Thus, we treat inhomogeneities close to the particle as small perturbations to the uniform state $\psi_h=1$: $\psi_0({\boldsymbol{r}},t)=1+\delta\psi_0({\boldsymbol{r}},t)$. Combining the two types of perturbations, and using the relationships of the wavefunction to the density, velocity and current (Eq. (\[eq:current\]), which in dimensionless units reads $\rho{\boldsymbol{v}}=(\psi^*\nabla\psi-\psi\nabla\psi^*)/(2i)$) we find $$\begin{aligned}
\label{eq:density_perturb}
\psi &=& 1 +\delta \psi_0 +g_p \delta\psi_1   \\
\rho &=& 1 +\delta \rho_0+g_p \delta\rho_1,   \\
{\boldsymbol{v}} &=& \delta{\boldsymbol{v}}^{(0)} + g_p \delta {\boldsymbol{v}}^{(1)},\end{aligned}$$ where $$\begin{aligned}
&&\delta \rho_0 =
\delta\psi_0+\delta\psi_0^*, \quad  \delta\rho_1 =
\delta\psi_1+\delta\psi_1^* ,\label{eq:density_linear} \\
&&\delta {\boldsymbol{v}}^{(0)} = \frac{1}{2 i}\nabla \left(\delta \psi_0 - \delta \psi_0^*\right), \quad
\delta{\boldsymbol{v}}^{(1)} =\frac{1}{2 i}\nabla (\delta\psi_1 - \delta\psi_1^*).\nonumber\\
\label{eq:velocity_linear}\end{aligned}$$

Combining Eq. (\[eq:fp2\]) with the expressions for the density perturbations, we have that the total force can be split into the contribution from the density variations in the BEC by causes external to the particle (initial preparation, stirring forces in $V_{ext}$, ...), and the density perturbations due to the presence of the particle ${\boldsymbol{F}}_p=
{\boldsymbol{F}}^{(0)}+{\boldsymbol{F}}^{(1)}$: $$\begin{aligned}
{\boldsymbol{F}}^{(0)}(t) = - \frac{g_p}{2\pi a^2} \int d^2{\boldsymbol{r}} e^{- \frac{({\boldsymbol{r}}-{\boldsymbol{r}}_p(t))^2}{2a^2}}
\nabla \delta\rho_0({\boldsymbol{r}},t),
\label{eq:fp_unperturb}
\\
{\boldsymbol{F}}^{(1)}(t) = - \frac{g_p^2}{2\pi a^2} \int d^2{\boldsymbol{r}} e^{- \frac{({\boldsymbol{r}}-{\boldsymbol{r}}_p(t))^2}{2a^2}}
\nabla \delta\rho_1({\boldsymbol{r}},t).
\label{eq:fp_perturb}\end{aligned}$$ The perturbative splitting of the force in these two contributions is completely analogous to the corresponding classical-fluid case in the incompressible [@maxey1983equation] and in the compressible [@parmar2012equation] situations. The ${\boldsymbol{F}}^{(0)}$ contribution is the equivalent to the classical inertial or pressure-gradient force on a test particle, which does not disturb the fluid, in a inhomogeneous and unsteady flow. In the following it will be called the *inertial* force. The ${\boldsymbol{F}}^{(1)}$ contribution takes into account perturbatively the modifications on the flow induced by the presence of the particle, and it will be called the *self-induced drag* on the particle. To complete the comparison with the classical expressions [@maxey1983equation; @parmar2012equation], we need to express Eqs. (\[eq:fp\_unperturb\]) and (\[eq:fp\_perturb\]) in terms of the unperturbed velocity field ${\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)=\delta{\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)$ and of the particle speed ${\boldsymbol{V}}_p(t)={\boldsymbol{\dot r}}_p(t)$. We are able to do so in a general situation for the inertial force ${\boldsymbol{F}}^{(0)}$. For ${\boldsymbol{F}}^{(1)}$, we obtain analytical expressions in the simple case where the impurity is moving with a constant velocity in an otherwise uniform BEC.

The desired relationships between $\nabla \delta\rho_0$ and $\nabla \delta\rho_1$ in Eqs. (\[eq:fp\_unperturb\])-(\[eq:fp\_perturb\]), and $\delta{\boldsymbol{v}}^{(0)}$ and ${\boldsymbol{V}}_p$ will be obtained from the linearization of the dGPE Eq. (\[eq:GPe\_dimless\]) around the uniform steady state $\psi_h=1$: $$\begin{aligned}
\partial_t \delta\psi_0  &=& (i+\gamma)\left(\frac{1}{2}\nabla^2 -1\right)\delta\psi_0 \nonumber \\
&-& (i+\gamma)\delta\psi^*_0,
\label{eq:psi0lin} \\
\partial_t \delta\psi_1  &=&
(i+\gamma)\left(\frac{1}{2}\nabla^2-1\right)\delta\psi_1 \nonumber \\
&-&(i+\gamma)\delta\psi_1^* - (i+\gamma) \mathcal U_p({\boldsymbol{r}}-{\boldsymbol{r}}_p) \ .
\label{eq:psi1lin}\end{aligned}$$ Terms containing $V_{ext}$ are not included in Eq. (\[eq:psi0lin\]) because of our assumption of sufficient distance between possible stirring sources and the neighborhood of the particle position, the only region that–as we will see– will enter into the calculation of the forces. In the next sections we solve these linearized equations to relate density perturbations to undisturbed velocity field and particle velocity.

Inertial force {#sec:inertial}
--------------

To convert Eq. (\[eq:fp\_unperturb\]) for the inertial force into an expression suitable for comparison for the corresponding term in classical fluids, we need to express $\nabla\delta\rho_0$ in terms of the undisturbed velocity field ${\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)=\delta{\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)$. To this end, we substract Eq. (\[eq:psi0lin\]) from its complex conjugate, obtaining: $$\begin{aligned}
\left(\nabla^2-4\right)\nabla\delta\rho_0=
4 \left(\partial_t-\frac{\gamma}{2}\nabla^2\right)\delta{\boldsymbol{v}}^{(0)},
\label{eq:drhodv0}\end{aligned}$$ where we have used Eqs. (\[eq:density\_linear\]) and (\[eq:velocity\_linear\]). Since the force formulae require to obtain the condensate density in a neighborhood of the particle position, it is convenient to move to a coordinate frame with center always at the (possibly moving) particle location ${\boldsymbol{r}}={\boldsymbol{r}}_p(t)$. Thus we change variables from $({\boldsymbol{r}}, t)$ to $({\boldsymbol{z}}, t)$, with ${\boldsymbol{z}}={\boldsymbol{r}} -{\boldsymbol{r}}_p(t)$, and the velocity field will be now referred to the particle velocity ${\boldsymbol{V}}_p(t)={\boldsymbol{\dot r}}_p(t)$: $\delta {\boldsymbol{w}}^{(0)}({\boldsymbol{z}},t)=\delta{\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)- {\boldsymbol{V}}_p(t)$. Equation (\[eq:drhodv0\]) becomes: $$\begin{aligned}
\left(\nabla_z^2 - 4\right)&&\nabla_z\delta\rho_0 = \nonumber \\
&4& \left(\partial_t-{\boldsymbol{V}}_p\cdot\nabla_z-\frac{\gamma}{2}\nabla_z^2\right)\delta{\boldsymbol{w}}^{(0)} + {\boldsymbol{\dot V}}_p(t),
\label{eq:drhodv}\end{aligned}$$ which has the corresponding equation for its Green’s function given by $$\left(\nabla_z^2 - 4\right)G({\boldsymbol{z}}) = \delta({\boldsymbol{z}})$$ with the boundary condition $G(|{\boldsymbol{z}}|\rightarrow\infty)\rightarrow 0$ (corresponding to vanishing $\nabla_z \delta\rho_0({\boldsymbol{r}})$ at $|{\boldsymbol{r}}|=\infty$). The solution is given by the zeroth order modified Bessel function $G({\boldsymbol{z}})=-K_0(2|{\boldsymbol{z}}|)/(2\pi)$. Hence, the gradient of the density perturbation can be written as the convolution with the Green’s function:

$$\begin{aligned}
\nabla_z\delta\rho_0({\boldsymbol{z}},t)=
-\frac{2}{\pi}\int d{\boldsymbol{z}}' K_0(2|{\boldsymbol{z}}-{\boldsymbol{z}}'|)
\left[\left(\partial_t-{\boldsymbol{V}}_p\cdot\nabla_{{\boldsymbol{z}}'}-\frac{\gamma}{2}\nabla^2_{{\boldsymbol{z}}'}\right)\delta{\boldsymbol{w}}^{(0)}({\boldsymbol{z}}',t) + {\boldsymbol{\dot V}}_p(t) \right] ,
\label{eq:nablarho}\end{aligned}$$

and the expression for the force (\[eq:fp\_unperturb\]), using the comoving variables $({\boldsymbol{z}},t)$, becomes: $$\begin{aligned}
{\boldsymbol{F}}^{(0)}(t)=
-\frac{g_p}{\pi^2a^2}\int d{\boldsymbol{z}} e^{-\frac{z^2}{2 a^2}} \int d{\boldsymbol{z}}' K_0(2|{\boldsymbol{z}}-{\boldsymbol{z}}'|)
\left[\left(\partial_t-{\boldsymbol{V}}_p\cdot\nabla_{{\boldsymbol{z}}'}-\frac{\gamma}{2}\nabla^2_{{\boldsymbol{z}}'}\right)\delta{\boldsymbol{w}}^{(0)}({\boldsymbol{z}}',t) + {\boldsymbol{\dot V}}_p(t) \right] .
\label{eq:F0full}\end{aligned}$$

The above expression is a weighted average of contributions from properties of the fluid velocity in a neighborhood of the impurity center-of-mass position (${\boldsymbol{z}} = 0$ in the comoving frame). The size of this neighborhood is given by the combination of the range of the Bessel function kernel, which in dimensional units would be the correlation length $\xi$, and the range of the Gaussian potential, $a$, giving an effective particle size. In classical fluids, the analogous force on a spherical particle involves the average of properties of the undisturbed velocity field within the sphere size [@parmar2012equation], and there is no equivalent to the role of $\xi$.

As in the classical case [@maxey1983equation; @parmar2012equation], if fluid velocity variations are weak at scales below $a$ and $\xi$, we can approximate the condensate velocity by a Taylor expansion near the impurity, i.e.: $$\begin{aligned}
\delta w_i^{(0)}(\mathbf z',t)&\approx& \delta
w_i^{(0)}(t)+\sum_j
e_{ij}(t) z'_j \nonumber\\
&+&\frac{1}{2}\sum_{jk}e_{ijk}(t)  z'_j z'_k + \ldots,\end{aligned}$$ where the indices $i,j,k=x,y$ denote the coordinate components. $e_{ij}(t)=\partial_j \delta w_i^{(0)}(\mathbf z,t)|_{{\boldsymbol{z}}=0}$ and $e_{ijk}(t)=\partial_j
\partial_k \delta w_i^{(0)}(\mathbf z,t)|_{{\boldsymbol{z}}=0}$ are gradients of the unperturbed condensate relative velocity. Inserting this expansion into Eq. (\[eq:F0full\]), and performing the integrals of the Gaussian and of the Bessel function (using for example $\int K_0(2|{\boldsymbol{z}}|)d{\boldsymbol{z}}=\pi/2$ and $\int z_i z_j K_0(2|{\boldsymbol{z}}|)d{\boldsymbol{z}}
=(\delta_{ij}/2)\int_0^\infty 2\pi z^3
K_0(2z)dz=\delta_{ij}\pi/4$), we obtain: $$\begin{aligned}
{\boldsymbol{F}}^{(0)}(t) &\approx& g_p {\boldsymbol{\dot V}}_p(t)
+ g_p\left[\partial_t
    - {\boldsymbol{V}}_p(t)\cdot\nabla_z +\frac{a^2}{2}\partial_t\nabla^2_z    \right. \nonumber \\
&-& \left.   \frac{\gamma}{2}\nabla_z^2 +
    \frac{1}{4}\partial_t\nabla^2_z  \right] \delta{\boldsymbol{w}}^{(0)}({\boldsymbol{z}},t)|_{{\boldsymbol{z}}={\boldsymbol{0}}} \ .
    \label{eq:Fi_th}\end{aligned}$$ The terms containing Laplacians are analogous to the Faxén corrections in classical fluids [@maxey1983equation] which arise for particles with finite size. Here, they arise from a combination of the finite effective size of the particle, $a$, and of the quantum coherence length, $\xi=1$. This last effect remains even in the limit of vanishing particle size $a\rightarrow 0$. Interestingly, one of the two terms in these quantum corrections depend on $\gamma$ hence indirectly on the presence of the thermal cloud.

As in the classical case, if flow inhomogeneities are unimportant below the scales $a$ and $\xi$, we can neglect the Laplacian terms in Eq. (\[eq:Fi\_th\]). Returning to the variables $({\boldsymbol{r}},t)$ in the lab frame of reference, the terms containing ${\boldsymbol{V}}_p$ cancel out, showing that the inertial force is mainly given by the local fluid acceleration: $$\begin{aligned}
\label{eq:pd1}
{\boldsymbol{F}}^{(0)}(t) = g_p \partial_t \delta{\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t) \big|_{{\boldsymbol{r}}={\boldsymbol{r}}_p(t)}.
 \end{aligned}$$ We have assumed a small non-uniform unperturbed velocity field ${\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)= \delta{\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)$. To leading order in this small velocity, the partial derivative $\partial_t\delta{\boldsymbol{v}}^{(0)}$ and the material derivative $D\delta{\boldsymbol{v}}^{(0)}/Dt=\partial_t\delta{\boldsymbol{v}}^{(0)}+\delta{\boldsymbol{v}}^{(0)}\cdot \nabla \delta{\boldsymbol{v}}^{(0)}$ are identical. In classical fluids the same ambiguity occurs and it has been established, on physical grounds and by going beyond linearization, that using the material derivative is more correct [@maxey1983equation]. After all, using this material derivative in the equation of motion simply means that, under the above approximations and in places where stirring and other external forces are absent, the local acceleration on the impurity arises from the corresponding acceleration of the condensate. Since for $a\rightarrow 0$ the condensate-impurity interaction has a similar scattering potential (delta function) as that for the interaction between condensate particles, similar accelerations would be experienced by a condensate particle and by the impurity, just modulated by a different coupling constant. Thus, replacing $\partial_t$ by $D/Dt$ in (\[eq:pd1\]) the approximate inertial force becomes: $$\begin{aligned}
\label{eq:pd}
{\boldsymbol{F}}^{(0)}(t)  &=&
g_p \left.\frac{D {\boldsymbol{v}}^{(0)}}{Dt}\right|_{{\boldsymbol{r}}={\boldsymbol{r}}_p(t)}\ ,
\label{eq:InertialSimpleDimless}
 \end{aligned}$$ or, if we return back to dimensional variables: $$\begin{aligned}
\label{eq:pd_dim}
{\boldsymbol{F}}^{(0)}(t)  &=&
\frac{g_p}{g}  m \left.\frac{D {\boldsymbol{v}}^{(0)}}{Dt}\right|_{{\boldsymbol{r}}={\boldsymbol{r}}_p(t)}\ .
\label{eq:InertialSimpleDim}
 \end{aligned}$$ This is equivalent to the equation for the inertial force in classical fluids [@maxey1983equation] except that the coefficient of the material derivative in the classical case is the mass of the fluid fitting in the size of the impurity. In the comoving frame, replacement of the partial by the material derivative amounts to replace $(\partial_t-{\boldsymbol{V}}_p\cdot
\nabla_{{\boldsymbol{z}}'})\delta{\boldsymbol{\omega}}^{(0)}$ in Eq. (\[eq:F0full\]) by $D\delta{\boldsymbol{\omega}}^{(0)}/Dt$. Eq. (\[eq:pd\]) is expected to be valid for small values of $g_p$ and in regions where fluid velocity and density inhomogeneities are both small and weakly varying. At this level of approximation neither compressibility nor dissipation effects appear explicitly in the inertial force, in analogy with classical compressible fluids [@parmar2012equation]. But these effects are indirectly present by determining the structure of the field ${\boldsymbol{v}}^{(0)}({\boldsymbol{r}},t)$.

Self-induced drag force {#sec:drag}
-----------------------

In classical fluids, consideration of the self-induced force on a particle moving at arbitrary time-dependent velocity, coming from the perturbation in produces in the flow, leads to different terms, namely [@maxey1983equation; @parmar2012equation] the viscous (Stokes) drag, the unsteady-inviscid term that in the incompressible case becomes the added-mass force, and the unsteady-viscous term that in the incompressible case becomes the Basset history force. They are expressed in terms of the undisturbed velocity flow ${\boldsymbol{v}}^{(0)}$ and the particle velocity ${\boldsymbol{V}}_p(t)$. Here, for the BEC case, we are able to obtain the self-induced force only for a particle moving at constant speed on the condensate. For the classical fluid case, in this situation the only non-vanishing force is the Stokes drag, so that this is the force we have to compare our quantum result with. We note that the condensate itself in the absence of the particle perturbation can be in any state of (weak) motion since in our perturbative approach summarized in Eqs (\[eq:psi0lin\])-(\[eq:psi1lin\]), the inhomogeneity $\delta\psi_0$ and the $g_p$-perturbation $\delta\psi_1$ are uncoupled.

It is convenient to transform the problem to the frame of reference moving with the particle $({\boldsymbol{r}},t)\rightarrow ({\boldsymbol{z}},t)$ with ${\boldsymbol{z}}={\boldsymbol{r}} - {\boldsymbol{r}}_p(t)$, so that Eq. (\[eq:psi1lin\]) becomes $$\begin{aligned}
\partial_t \delta\psi_1 - {\boldsymbol{V}}_p \cdot \nabla \delta\psi_1 &=&
(i+\gamma)\left(\frac{1}{2}\nabla^2-1\right)\delta\psi_1 \nonumber \\
&-&(i+\gamma)\delta\psi_1^* - (i+\gamma) \mathcal U({\boldsymbol{r}}) \ .
\label{eq:psi1linComoving}\end{aligned}$$ Note that such Galilean transformations of the GPE using a constant ${\boldsymbol{V}}_p$ are often accompanied by a multiplication of the transformed wavefunction by a phase factor $\exp(i{\boldsymbol{V}}_p \cdot {\boldsymbol{z}} + \frac i 2 V_p^2 t)$, in order to transform the condensate velocity (see below) to the new frame of reference, and account for the shift in kinetic energy. Indeed, such a combined transformation leaves the GPE unchanged at $\gamma = 0$ [@pismen1999] (but not for $\gamma>0$). The density perturbation $\delta\rho_1$ is already given correctly by $\delta\psi_1+\delta\psi_1^*$, where $\delta\psi_1({\boldsymbol{z}},t)$ is the solution of (\[eq:psi1linComoving\]), without the need of any additional phase factor. The velocity in the comoving frame would need to be corrected as $\delta\omega^{(1)}({\boldsymbol{z}},t)=\delta {\boldsymbol{v}}^{(1)}-{\boldsymbol{V}}_p$, with $\delta v^{(1)}$ given by expression (\[eq:velocity\_linear\]) in terms of he solution of (\[eq:psi1linComoving\]).

Eq. (\[eq:psi1linComoving\]) in the steady-state can be solved by using the Fourier transform $\delta\psi_1({\boldsymbol{z}})
=1/(2\pi)^2 \int d^2{\boldsymbol{k}} e^{i{\boldsymbol{k}}\cdot {\boldsymbol{z}}}
\delta\hat\psi_1({\boldsymbol{k}})$. It follows that the linear system of equations for $\delta\hat\psi_1({\boldsymbol{k}})$ and $\delta\hat\psi_1^*(-{\boldsymbol{k}})$ is given by $$\begin{aligned}
\left[-2i {\boldsymbol{k}}\cdot{\boldsymbol{V}}_p +(i+\gamma) (k^2+2)\right]&\delta\hat\psi_1& + 2(i+\gamma)\delta\hat\psi_1^* = \nonumber \\
&-& 2(i+\gamma) e^{-\frac{a^2k^2}{2}} ,\nonumber\\
\left[-2i {\boldsymbol{k}}\cdot{\boldsymbol{V}}_p +(-i+\gamma) (k^2+2)\right]&\delta\hat\psi_1^*& + 2(-i+\gamma)\delta\hat\psi_1 =  \nonumber \\
&-& 2(-i+\gamma) e^{-\frac{a^2k^2}{2}} .\nonumber\\\end{aligned}$$ By solving these equations, we find $\delta\hat\psi_1({\boldsymbol{k}})$ and $\delta\hat\psi_1^*(-{\boldsymbol{k}})$, and the Fourier transform of the density perturbation $\delta\rho_1 =
\delta\psi_1^*+\delta\psi_1$ then follows as $$\delta\hat\rho_1 = \frac{e^{-\frac{k^2a^2}{2}} (4 k^2(1+\gamma^2)-8i\gamma{\boldsymbol{k}}\cdot {\boldsymbol{V}}_p)}{4{\boldsymbol{k}}\cdot{\boldsymbol{V}}_p({\boldsymbol{V}}_p \cdot {\boldsymbol{k}}+i\gamma k^2 +2i\gamma)- k^2(4+k^2)(1+\gamma^2)}.
    \label{Rho1Delta}$$ Using the convolution theorem, we can express the self-induced force (\[eq:fp\_perturb\]) (in the co-moving frame, i.e. with $\mathbf r_p=0$) in terms of $\delta\hat\rho_1$ as $${\boldsymbol{F}}^{(1)}= - \frac{g_p^2}{(2\pi)^2}\int d^2{\boldsymbol{k}} i{\boldsymbol{k}}\delta\hat{\rho}_1({\boldsymbol{k}}) e^{-\frac{k^2a^2}{2}} .$$ This force can be decomposed into the normal and tangential components relative to the particle velocity ${\boldsymbol{V}}_p$: $ {\boldsymbol{F}}^{(1)} = F_{\|} {\boldsymbol{e}}_{\|}+F_{\perp} {\boldsymbol{e}}_{\perp}$. Due to symmetry, the normal component vanishes upon polar integration, and we are left with the tangential, or drag, force

$$\begin{aligned}
  F_{\|}&=& -\frac{g_p^2}{(2\pi)^2} \int_0^\infty dk\int_0^{2\pi} d\theta e^{-k^2a^2}\frac{i k^2\cos(\theta) \left[4 k^2(1+\gamma^2) - 8i\gamma kV_p\cos(\theta)\right]}{4kV_p\cos(\theta)(kV_p\cos(\theta)+ i\gamma k^2 + 2i\gamma) - k^2(4 + k^2)(1+\gamma^2)}.
\label{eq:force_int}\end{aligned}$$

$V_p$ is the modulus of ${\boldsymbol{V}}_p$. At zero temperature, i.e. when $\gamma=0$, the drag force reduces to the one that has also been calculated for a point particle in Refs. [@astrakharchik2004motion] and in [@pinsker2017gaussian] for a finite-$a$ particle: $$F_\parallel = - \frac{g_p^2}{\pi^2}\int_0^\infty dk\int_0^{2\pi} d\theta\frac{i k^2 \cos{\theta}e^{-k^2a^2}}{4V_p^2 \cos^2{\theta}- (4+k^2)} \ ,$$ which is zero for particle speed smaller than the critical value given by the long-wavelength sound speed, $V_p < c=1$. Above the critical speed, the integral has poles and acquires a non-zero value given by $$F_\parallel = -\frac{g_p^2 k_{max}^2}{4V_p}  e^{-\frac{a^2 k_{max}^2}{2}} \left[I_0\left(\frac{a^2 k_{max}^2}{2}\right) - I_1\left(\frac{a^2 k_{max}^2}{2}\right) \right]
    \label{eq:forcegamma0}$$ in terms of the modified Bessel functions of the first kind $I_n(x)$ and where $k_{max} =2\sqrt{V_p^2 -1}$. For vanishing $a$ the dominant term is proportional to $(V_p^2-1)/V_p$  [@astrakharchik2004motion]. This drag is pertaining to energy dissipation by radiating sound waves in the condensate away from the impurity. We stress again that we assume $a$ small enough such that emission of other excitations, such as vortex pairs, does not occur. It is important to note [@astrakharchik2004motion; @pinsker2017gaussian] that in order to obtain a real value for the force in Eq. (\[eq:forcegamma0\]) one has to consider that it has been obtained from the limit $\gamma\rightarrow 0^+$ in (\[eq:force\_int\]), which implies that an infinitesimal positive imaginary part needs to be considered in the denominator to properly deal with the poles in the integral.

![image](Figure_1.pdf){width="\textwidth"}

In general, for a non-zero $\gamma$, Eq. (\[eq:force\_int\]) simplifies upon an expansion in powers of $V_p$ to the leading order. For the linear term in $V_p$, we can perform the polar integration and arrive at $$F_\parallel = -\frac{2}{\pi}V_p \frac{\gamma}{1+\gamma^2}g^2_p \int \frac{k^3 e^{-a^2k^2}}{(4 +k^2)^2}dk \ .$$ Substituting $u = a^2(k^2 + 4)$, we find $$\begin{aligned}
  F_{\|}  =-V_p \frac{\gamma}{1+\gamma^2}g_p^2 \frac{1}{\pi} \left[ e^{4a^2 }E_1(4a^2)(1+4a^2) - 1\right], \nonumber\\
\label{eq:drag_coeff_full}\end{aligned}$$ where $E_1(x)$ denotes the positive exponential integral. When $a \to 0$, the expression inside the bracket diverges as $-\gamma_E - 1 - \ln(4a^2)$ with $\gamma_E$ begin the Euler-Mascheroni constant. It is therefore necessary to keep a finite size $a$.

This drag force is analogous to the viscous Stokes drag force in classical fluids since it is due to an effective interaction of the impurity with the normal fluid through the thermal drag on the BEC. The effective drag coefficient depends on the thermal drag such that it vanishes at zero temperature. But it also depends non-trivially on the size of the impurity and it diverges in the limit of point-like particle. Faxén corrections involving derivatives of the unperturbed flow are not present here because of the decoupling between $\delta\psi_0$ and $\delta\psi_1$ arising in the perturbative approach leading to (\[eq:psi0lin\])-(\[eq:psi1lin\]).

Numerical results {#sec:numerics}
=================

![image](Figure_2.pdf){width="\textwidth"}

To test the analytical predictions of the inertial force and the self-induced drag deduced above from the total force expression Eq. (\[eq:fp2\]), we performed numerical simulations of the dGPE. Actually, our simulations are done in the co-moving frame of the impurity moving at constant velocity ${\boldsymbol{V}}_p$, so that the equation we solve is (see numerical details in the Appendix): $$\begin{aligned}
&\partial_t \psi - {\boldsymbol{V}}_p \cdot \nabla \psi  =(i+\gamma)\left[\frac{1}{2}\nabla^2 \psi +\left(1-g_p\mathcal U_p - |\psi|^2\right)\psi\right], \nonumber \\
     \label{eq:ComovingdGPE}\end{aligned}$$ where the impurity is described by the Gaussian potential of intensity $g_p=0.01$ and effective size $a=\xi=1$, and is situated in the middle of the domain with the coordinates $x/\xi= 128$ and $y/\xi= 64$. As an initial condition, we start with the condensate being at rest and in equilibrium with the impurity. This is done by imaginary time integration of Eq. (\[eq:ComovingdGPE\]) for $V_p=0$ and $\gamma=0$. Then, at $t=0$, we solve the full Eq. (\[eq:ComovingdGPE\]), and as a consequence, sound waves are emitted from the neighborhood of the impurity (the size of the impurity is below the critical size for vortex nucleation). Their speed is determined by the dispersion relation $\omega({\boldsymbol{k}})$ giving the frequency as a function of the wavenumber and can be obtained by looking for plane-wave solutions to Eq. (\[eq:psi0lin\]). If $\gamma=0$, $\omega({\boldsymbol{k}})$ is given by the Bogoliubov dispersion relation [@BEC-review] $\omega({\boldsymbol{k}}) =
k\sqrt{1+k^2/4}$ (with $c=\xi=1$). Note that the smallest velocity, $c=1$, is that of long waves, and that waves of smaller wavelength travel faster. For $\gamma>0$, the planar waves are dampened out and the dispersion relation becomes $\omega({\boldsymbol{k}})=-i\gamma (k^2/2+1)+\sqrt{k^2+k^4/4-\gamma^2}$. The damping rate is determined by $\gamma$ and increases quadratically with the wavenumber. Also, in this case all waves have a group velocity faster than a minimum one that for small $\gamma$ is close to $c=1$.

When $V_p<1$ all the waves escape the neighborhood of the impurity (see an example in Fig. \[fig:density\_gamma\](a)) and are dissipated in a boundary buffer region that has large $\gamma$ (see numerical details in the Appendix and Supplemental Material [@SupplMat]). After a transient the condensate achieves a steady state when seen in the frame comoving with the impurity. Fig. \[fig:density\_gamma\](b) shows a steady spatial configuration for $\gamma=0$ and $V_p=0.9$. Figs. \[fig:density\_gamma\](d-e) show different profiles of the condensate density along the $x$ direction across the impurity position for $V_p$. The condensate density is depleted near the impurity due to the repulsive interaction, and its general shape depends on the speed $V_p$ and thermal drag $\gamma$. If $\gamma=0$ and $V_p\leq 1$ the density of this steady state has a rear-front symmetry with respect to the particle position (see specially Fig. \[fig:density\_gamma\](d)), so that under integration in Eq. (\[eq:fp2\]) the net force is zero. The presence of dissipation ($\gamma>0$) breaks this symmetry even if $V_p<1$ so that a net drag will appear in agreement with the calculation of Sect. \[sec:drag\]. When $V_p>1$ there are waves that can not escape from the neighborhood of the impurity, forming fringes in front of it and a wake behind it. (see Fig. \[fig:density\_gamma\](c) and \[fig:density\_gamma\](f) and Supplemental material [@SupplMat]). Similar modulations in the condensate density around an obstacle in supersonic flows has also been observed experimentally [@carusotto2006bogoliubov]. As a consequence of this rear-front asymmetry, drag is nonvanishing even when $\gamma=0$.

Movies showing the transient and long-time density behavior for several values of $V_p$ and at $\gamma=0$ are included as Supplemental Material [@SupplMat]. They are presented in the comoving impurity frame, so that the impurity appears static. The fluid suddenly stars to move towards the negative $x$ direction, and density approaches a steady state after the transient. Note that during all the dynamics, the density deviation with respect to the equilibrium value $\rho=1$ is very small, justifying the perturbative approach of Sect. \[sec:perturbation\] for this situation. The time evolution for $\gamma>0$ is qualitatively similar to the $\gamma=0$ one shown in the movies, except that the waves become damped and that there is a front-rear asymmetry in the steady state.

Our numerical setup is well suited to measure the force produced by the perturbation of the impurity on the fluid, i.e. the self-induced drag. Nevertheless, in the absence of the impurity the unperturbed state is the trivial $\psi=1$, so that $\delta\psi_0=0$ and the inertial force is identically zero. In order to test the accuracy of our expressions for the inertial force without the need of additional simulations under a different set-up, we still use the computed condensate density and velocity dynamics, produced by the impurity introduced in the system at $t=0$, but we evaluate the inertial force exerted by this flow on another test particle located at a different position. In fact, there is no need to think on the flow as being produced by an impurity: it can be produced by a moving laser beam that can modeled by an external potential $V_{ext}$ and the only impurity present in the system is the test particle on which the force is evaluated. In the following we evaluate the inertial and the self-induced drag forces on the different particles from the general expressions Eqs. (\[eq:fp\_unperturb\])-(\[eq:fp\_perturb\]) and from the approximate expressions of Sects. \[sec:inertial\] and \[sec:drag\].

Numerical evaluation of the inertial force {#sec:numerical-inertial}
------------------------------------------

We consider a test particle traveling at the same speed $V_p$ as the impurity or laser beam producing the flow, but located at a distance of 10 coherence lengths in front of it, and 20 coherence lengths in the $y$ direction apart from it. This distance is sufficient to avoid inclusion of $\mathcal U_p$ or $V_{ext}$ in Eq. (\[eq:psi0lin\]) for the neighborhood of the test particle. Condensate and test particle interact via a coupling constant $g_p'$ sufficiently small so that the full force on the later, Eq. (\[eq:fp2\]), is well approximated by the inertial part Eq. (\[eq:fp\_unperturb\]), being the perturbation the particle induces on the flow, and thus the force (\[eq:fp\_perturb\]) completely negligible.

Figure (\[fig:inertial\_force\]) shows, for different values of $V_p = 0.1, 0.8, 1$ at $\gamma=0$, the $x$ component of the time-dependent force produced by the transient flow inhomogeneities hitting the test particle in the form of sound waves. The size of the test particle, taking several values, is called $a'$ to distinguish it from the size $a$ of the particle producing the flow perturbation. Blue lines are computed from the exact Eq. (\[eq:fp2\]) or equivalently from Eq. (\[eq:fp\_unperturb\]) to which it reduces for sufficiently small $g_p'$. Because of the rather explicit appearance of the interaction potential in this formula, we label the blue lines in Fig. (\[fig:inertial\_force\]) as ‘potential force’. High frequency waves arrive before low-frequency ones, because its larger sound speed. We also see how the frequencies become Doppler-shifted for increasing $V_p$. We have derived in Sect. \[sec:inertial\] several approximate expressions for the inertial force. First, Eq. (\[eq:F0full\]) is obtained with the sole assumption (besides $g_p$ sufficiently small) of smallness of the unsteady and/or inhomogeneous part $\delta\psi_0$ of the wavefunction, which allows linearization. Eq. (\[eq:Fi\_th\]) assumes in addition weak inhomogeneities below scales $a$ and $\xi$, and finally Eqs. (\[eq:InertialSimpleDimless\]) and (\[eq:InertialSimpleDim\]) (equivalent under the previous linearization approximation) completely neglects such inhomogeneities (or equivalently, they correspond to $a,\xi\rightarrow 0$). We show as black lines in Fig. (\[fig:inertial\_force\]) the prediction of this last approximation, similar to the most standard classical expressions. Since we have computed the wavefunction $\psi=1+\delta\psi_0$ in the comoving frame from Eq. (\[eq:ComovingdGPE\]), we actually use expression (\[eq:Fi\_th\]) without the Faxén Laplacian terms, with $\delta{\boldsymbol{\omega}}^{(0)}=\nabla(\delta \psi_0-\delta
\psi_0^*)/(2i) - {\boldsymbol{V}}_p$, and ${\boldsymbol{\dot V_p}}=0$. Fig. (\[eq:ComovingdGPE\]) shows that the full force computed from Eq. (\[eq:fp2\]) is well-captured by the approximate expression of the inertial force for small test-particle size $a'$. Accuracy progressively deteriorates for increasing $a'$, and also for increasing $V_p$, but this classical expression remains a reasonable approximation until $a'\approx 1$. The accuracy can be improved by considering higher-order Faxén corrections, Eq. (\[eq:Fi\_th\]), or even better, by considering the integral form in Eq. (\[eq:F0full\]). We have explicitly checked that keeping the full Gaussian integration in Eq. (\[eq:F0full\]) but approximating the integrand in the Bessel integral by its value at the particle position gives a very good approximation to the exact force even for $a'=1$.

Numerical evaluation of the drag force {#sec:numerical-drag}
--------------------------------------

![Plot of the drag force in the steady-state regime as function of the constant speed $V_p$. Dashed lines are the analytical predictions based on Eq. (\[eq:forcegamma0\]) (for $\gamma=0$) and Eq. (\[eq:force\_int\]) (for $\gamma>0$). The symbols correspond to the numerically computed force from Eq. (\[eq:fp2\]) based on direct simulations of the dGPE Eq. (\[eq:ComovingdGPE\]). The inset figure shows the small $V_p$ behavior, with solid straight lines giving the linear dependence of the drag force on the speed for $\gamma>0$, in the small-$V_p$ approximation given by Eq. (\[eq:drag\_coeff\_full\]). We use $a=1$ and $g_p=0.01$.[]{data-label="fig:force_drag"}](Figure_3.pdf){width="50.00000%"}

We now return to the situation in which there is a single impurity in the system, with size $a=\xi=1$ and $g_p=0.01$. It moves in the positive $x$ direction with speed $V_p$ producing a perturbation on the uniform and steady condensate state $\psi=1$. We compute it in the comoving frame, in which the particle is at rest and fluid moves with speed $-V_p$, by using Eq. (\[eq:ComovingdGPE\]). Since in the absence of the impurity there is no inhomogeneity nor time dependence, $\delta\psi_0=0$ and the exact force on the impurity, Eq. (\[eq:fp2\]), is also given by the self-induced drag expression given by Eq. (\[eq:density\_perturb\]). After a transient, that in analogy with the results for compressible classical fluids [@Longhorn1952; @parmar2012equation] we expect to be of the order of the time needed by the sound waves to cross a region of size $a$ or $\xi$, the condensate density near the particle achieves a steady state in the comoving frame, and we then measure the steady drag on the particle. Figure (\[fig:force\_drag\]) shows this force, for several values of $V_p$ and $\gamma$, as dots. The approximate value of the drag force that is obtained under the assumption of small perturbation (small $g_p$) that allows linearization is shown as dashed lines. It is computed from Eq. (\[eq:forcegamma0\]) for $\gamma=0$ and Eq. (\[eq:force\_int\]) for $\gamma>0$. the agreement is excellent. As shown in the inset figure, in the regime of small velocities, the self-induced drag is indeed linearly dependent on the speed with an effective drag coefficient that is well captured by Eq. (\[eq:drag\_coeff\_full\]). This Stokes-like drag at small speeds is due to the thermal drag of the condensate by the normal fluid, quantified by $\gamma$. We notice that the dependence of the drag force on $V_p$ is consistent with having a critical velocity for superfluidity even at $\gamma>0$, in the sense that there is still a relatively abrupt change in the force (sharper for smaller $\gamma$) around a particular impurity speed. The superfluidity of BECs at finite temperature is still an open question. Recent experiments [@sing2016probingsuperfluidity; @Weimer2015_criticavelocity] report superfluid below a critical velocity which is related to the onset of fringes [@wouters2010]. In the dGPE, the steady state drag is always nonzero. Nonetheless, there is a critical velocity above which acoustic emission significatively increases that can be associated to the breakdown of superfluidity. This is the regime where the drag force is dominated by the interaction of the impurity with the supersonic shock waves that are reminiscent of the Kelvin wake in classical fluids. This is seen in Fig. \[fig:density\_gamma\](c) and observed experimentally [@carusotto2006bogoliubov]. The maximum drag force occurs near the velocity for which the cusp lines forming the wake still retain an angle close to $\pi$. With increasing speed, this angle becomes more acute (Fig. \[fig:density\_gamma\](f)), and this lowers the density gradient around the impurity.

Conclusions
===========

We have studied, from analytic and numerical analysis of the dGPE, the hydrodynamic forces acting on a small moving impurity suspended in a 2D BEC at finite temperature. In the regime of small coupling constant $g_p$ and thermal drag $\gamma$, the force arising from the gradient of the condensate density can be decomposed onto the inertial force that is produced by the inhomogeneities and time-dependence of the condensate in the absence of the particle, and the self-induced force which is determined by the perturbation produced by the impurity on the condensate. When the unperturbed flow can be considered homogeneous on scales below the particle size and the condensate coherence length, the classical Maxey and Riley expression [@maxey1983equation], giving the inertial force in terms of the local or material fluid acceleration, is a good description of the force. When inhomogeneities become relevant below these scales, Faxén-type corrections arise, similar to the classical ones in the presence of a finite-size particle, but here the coherence length plays a role similar to the particle size. In addition, the condensate thermal drag enters into these expressions, at difference with the classical viscous case. We also determined the self-induced force in the steady-state regime and shown that it is non-zero at any velocity $V_p$ of the moving impurity if $\gamma>0$. For small $V_p$, this force is given as a Stokes drag which is linearly proportional to $V_p$ with a drag coefficient dependent on the thermal drag $\gamma$. With increasing velocities, there are corrections to the linear drag and above a critical speed of the order of $V_p=c=1$, the self-induced drag is dominated by the interactions of the impurity with the emitted shock waves. Following the recent experimental progress on testing the superfluidity in BEC at finite temperature [@singh2016probing], it would be interesting to test experimentally our prediction of the linear drag on the impurity due to the condensate thermal drag at small velocities by using measurements of the local heating rate.

We have checked our analytical expressions with numerical simulations in the situation in which the impurity moves at constant velocity, possibly driven by external forces different from the hydrodynamic ones analyzed here. When the coupling constant $g_p$ is sufficiently small so that only the inertial force is relevant, the equation of motion of the impurity under the sole action of the inertial force would be $m_p d{\boldsymbol{V}}_p(t)/dt={\boldsymbol{F}}^{(0)}(t)$, with $m_p$ the mass of the particle and ${\boldsymbol{F}}^{(0)}(t)$ one of the suitable approximations to the inertial force given in Sect. \[sec:inertial\]. For larger $g_p$, when the condensate becomes distorted by the impurity, we have computed the self-induced drag only in the steady case, so that can not write a general equation of motion for the impurity in interaction with a time-dependent perturbed flow. In analogy with classical compressible flows [@Longhorn1952; @parmar2012equation], we expect history-dependent forces in this unsteady situation. The dependence on the thermal drag, however, would be quite different from that of viscous classical fluids, because of the lack of viscous boundary layers in the BEC case.

In this study, we have focused on a small impurity that can only shed acoustic waves. Another interesting extension of this would be to further investigate the drag and inertial forces for larger impurity sizes, which can emit vortices, and study the effect of vortex-impurity interactions on the hydrodynamics forces.

Acknowledgements {#acknowledgements .unnumbered}
================

We are thankful to Vidar Skogvoll, Kristian Olsen, Zakarias Laberg Hejlesen and Per Arne Rikvold for stimulating discussions. This work was partly supported by the Research Council of Norway through its centers of Excellence funding scheme, Project No. 262644, and by Spanish MINECO/AEI/FEDER through the María de Maeztu Program for Units of Excellence in R&D (MDM- 2017-0711).

Appendix: Numerical integration of dGPE {#appendix-numerical-integration-of-dgpe .unnumbered}
=======================================

![Simulation domain showing the buffer region, outside the main simulation region, in which thermal drag is greatly enhanced to eliminate the emitted waves sufficiently far from the moving particle (which is at $x/\xi=128$, $y/\xi=64$). The density shown is the steady state (in the comoving frame, hence the direction of the arrows indicating the flow velocity in this frame) for $V_p=1.6$ and $\gamma=0$. []{data-label="fig:BufferRegion"}](Figure_4.pdf){width="45.00000%"}

Numerical simulations of dGPE Eq. (\[eq:ComovingdGPE\]) are run for a system size of $128\times 256$ (in units of $\xi$) corresponding to the grid size $dx=0.25\xi$, and $dt=0.01\xi/c$. To simulate an infinite domain where the density variations emitted by the impurity do not recirculate under periodic boundary conditions, we use the fringe method from [@reeves2015identifying]. This means that we define buffer (fringe) regions around the outer rim of the computational domain (see Fig. \[fig:BufferRegion\]) where the thermal drag $\gamma$ is much larger than its value inside the domain, such that any density perturbation far from the impurity is quickly damped out and a steady inflow is maintained. The thermal drag becomes thus spatially-dependent and given by $\gamma({\boldsymbol{r}}) = \max[\gamma(x),\gamma(y)]$, where $$\begin{aligned}
&\gamma(x)= \frac{1}{2}\big(2 + \tanh{[(x-x_p-w_x)/d]}\nonumber\\
&-\tanh{[(x-x_p+w_x)/d]}\big) + \gamma_0,\end{aligned}$$ and similarly for $\gamma(y)$. Here ${\boldsymbol{r}}_p = (x_p,y_p)=
(128\xi,64\xi)$ is the position of the impurity and $\gamma_0$ is the constant thermal drag inside the buffer regions (bulk region). We set the fringe domain as $w_x=100\xi$, $w_y=50\xi$ and $d=7\xi$ as illustrated in Figure \[fig:BufferRegion\].

By separating the linear and non-linear terms in Eq. (\[eq:ComovingdGPE\]), we can write the dGPE formally as [@audunsthesis] $$\partial_t \psi = \hat\omega(-i\nabla)\psi + N({\boldsymbol{r}},t),$$ where $\hat\omega(-i\nabla) = i[\frac 1 2 \nabla^2+1]+{\boldsymbol{V}}_p\cdot \nabla$ is the linear differential operator and $N({\boldsymbol{r}},t)=-(i+\gamma)(\mathcal U_p+|\psi|^2)\psi + \gamma
\psi +\frac{1}{2}\gamma\nabla^2\psi$ is the nonlinear function including the spatially-dependent $\gamma$ and $\mathcal U_p$. Taking the Fourier transform, we obtain ordinary differential equations for Fourier coefficients $\psi({\boldsymbol{k}},t)$ as $$\partial_t\hat\psi({\boldsymbol{k}},t) = \hat\omega({\boldsymbol{k}})\hat\psi({\boldsymbol{k}}, t) + \hat N({\boldsymbol{k}} ,t),
    \label{eq:DecompositionGPE}$$ which can be solved by an operator-splitting and exponential-time differentiating method [@cox2002exponential]. It means that we exploit the fact that the linear part of Eq. (\[eq:DecompositionGPE\]) can be solved exactly by multiplying with the integrating factor $e^{-\hat\omega({\boldsymbol{k}})t}$. This leads to $$\partial_t \left(\hat\psi({\boldsymbol{k}},t) e^{-\hat\omega({\boldsymbol{k}})t}\right) = e^{-\hat\omega({\boldsymbol{k}})t}\hat N({\boldsymbol{k}}, t).
    \label{eq:ETD}$$ The nonlinear term $\hat N({\boldsymbol{k}},t)$ is linearly approximated in time for a small time-interval $(t,t+\Delta t)$, i.e $$\hat N({\boldsymbol{k}}, t+\tau) = N_0 + \frac{N_1}{\Delta t}\tau$$ where $N_0 = \hat N(t)$ and $N_1 = \hat N(t+\Delta t) -N_0$. Inserting this into Eq. (\[eq:ETD\]) and integrating from $t$ to $t+\Delta t$ we get $$\begin{aligned}
&\hat{\psi}({\boldsymbol{k}},t+\Delta t) = \hat{\psi}({\boldsymbol{k}},t) e^{\hat{\omega}({\boldsymbol{k}})\Delta t} + \frac{N_0}{\hat\omega({\boldsymbol{k}})}\left(e^{\hat \omega({\boldsymbol{k}}) \Delta t} -1\right) \nonumber\\
&+ \frac{N_1}{\hat\omega({\boldsymbol{k}})} \left[\frac{1}{\hat\omega({\boldsymbol{k}}) \Delta t}(e^{\hat\omega({\boldsymbol{k}}) \Delta t} -1) -1 \right].\end{aligned}$$ Since computing the value of $N_1$ requires knowledge of the state at $t+\Delta t$ before we have computed it, we start by setting it to zero and find a value for the state at $t +\Delta
t$ given that $\hat N(t)$ is constant in the interval. We then use this state to calculate $N_1$, and add corrections to the value we got when assuming $N_1=0$.
