Beyond superfluidity in driven non-equilibrium Bose-Einstein condensates

The phenomenon of superfluidity in open Bose-Einstein condensates (BEC) is analysed numerically and analytically. It is found that a superfluid phase is feasible even above the speed of sound, when forces due to inhomogeneous non-equilibrium processes oppose the contributions of homogeneous processes. Furthermore a regime of accelerating impurities can be observed for particular pumping/decay strategies. All findings are derived within the complex Gross-Pitaevskii (GP) theory with creation and annihilation terms. Utilising this framework the effective force acting on an impurity as it moves with velocity v through the open condensate can be calculated. The result shows that the force is continuously increasing/decreasing with increasing velocity starting from the state of zero motion at v = 0, a property that can be traced down to the additional homogeneous annihilation/creation term in the extended GP model. Our findings stand in stark contrast to the concept of a topological phase transition to frictionless flow below a critical velocity as observed for equilibrium Bose-Einstein condensates analytically, numerically and for trapped atoms experimentally.


INTRODUCTION
What defines superfluidity of a many-body system? The answer can be given in terms of a statement based on Landau's theory of superfluidity [7][8][9]: Below a certain critical velocity, due to non-existing energetically affordable elementary excitations within the many-body quantum system, an impurity moves dissipationless through the superfluid state of matter. Here the concept of drag force acting on the impurity due to interactions of the impurity with the fluid as it excites the quantum system turns out to be key as measure of dissipation.
Particularly an impurity which is moving with velocity v through a fluid in its quantum mechanical ground state can cause transitions from the fluid's ground state to excited states lying on the line ε = pv in the energymomentum space [9,10]. However if the whole energy spectrum of the fluid is above this line, the motion of the impurity cannot excite the system. This implies the superfluid phase where the impurity moves without resistance through the ensemble of unexcited matter particles. Even when the line ε = pv intersects the energy spectrum of the fluid in its ground state, transition probabilities to these states can be strongly suppressed due to Boson interactions or due to the nature of the external perturbing potential [10]. For all scenarios the drag force experienced by the impurity gives us a quantitative measure of the state of the fluid and below a critical velocity, if the fluid cannot be excited, the impurity experiences no drag [1,2] a phenomenon already envisaged in the classic papers in the beginning of the 20th century [11][12][13][14][15].
Experimentally superfluidity has been unambiguously observed for the condensed state in various weakly interacting and dilute effectively Bose gases of atoms or molecules at ultra-low temperatures in the nano Kelvin range [4,5,[16][17][18] and even in strongly interacting gases of 6 Li fermions [5]. On the other hand, theoretically, the ideal Bose gas, i.e. a quantum gas without interactions at all states of motion obeys dissipation of energy. Only for certain cases of interacting quantum systems such as in the BE condensed phase with its intrinsic nonlinearity due to particle interactions between Bosons, we observe a phase of superfluidity. Here excitations are suppressed when an obstacle moves through the BE condensed phase by the nonlinear self-interactions which results in absence of drag [3,[23][24][25]. More recently it has become clear that a mean-field analysis leaves out subtle quantum fluctuations [19][20][21][22], which, if taken into account e.g. in a linear response framework [10] or by directly considering the quantum correction in a quantum Bogoliubov analysis [22], give rise to non-classical drag forces even below the critical velocity due to mean field theory. However a rigorous analysis shows that the superfluid phase persists even on a quantum level [26,27].
Superfluidity can be tested experimentally and proven by considering the dissipationless flow around impurities [3,[23][24][25]28], even in the presence of quantum fluctuations [19,26,27]. This has been done theoretically in the semiclassical GP framework by the means of a Bogoliubov analysis for point-like [1] and subsequently Gaussian [2] weakly-interacting impurities. This analysis confirmed the existence of a superfluid phase in the leading order contribution -the drag force vanishes below the non-zero critical velocity that is equal the speed of sound for weakly-interacting obstacles, which hold as well quantum mechanically [26], while geometric features significantly alter its magnitude [2]. The analytical insights on nonequilibrium systems presented here build on this type of analysis and will be supported by numerical integration.
Once a superfluid is put in motion other aspects associated with this state of matter emerge. For nonequilibrium quasi-Bose-Einstein condensates such as Polariton condensates in their lowest energy state [29] the low arXiv:1611.03430v4 [cond-mat.quant-gas] 4 Apr 2017 scattering rates from defects moving at velocities below the speed of sound and the generation of Cherenkov-type waves at supersonic velocities have been observed experimentally [30], while reduced drag at subsonic speeds has been noted in [31,32]. This reduced drag has been explained by the finite lifetime of Bogoliubov modes due to drain present in this open system [31]. In addition it has been realised that elementary excitations known from equilibrium condensates exist too -dark solitons are feasible in 1d [34,35] or quantised vortices in 2d nonequilibrium condensates [36,37] above a specific critical velocity or even spontaneously due to purely nonequilibrium dynamics [38]. The mathematical extension of the governing partial differential equation for open BEC in those scenarios implies intrinsic adaptions of the excitations' mathematical form [34,34,39] (see e.g. [40][41][42][43][44][45][46] for rigorous results on equilibrium condensates). However several findings suggest similar response to motion or obstacles in relative motion to the open system [34,36,37,47]. Now in the semiclassical mean-field regime many Bosons in the ground state are described by solutions to the Gross-Pitaevskii equation (GPE) [7,50,51], while non-equilibrium condensates in their simplest form are described by an extended GPE with additional complex terms, which correspond to creation and annihilation operators of these modes [30,[47][48][49].
One goal of this paper is to point out implications of the non-equilibriumness of open systems on the possibility of leading order superfluidity by testing a quantum state in relative motion to an obstacle. Here we use the complex GP framework as model of the coherent manybody system, corresponding to many particles being in the same quantum state, and we study the behaviour of Dirac and finite-sized obstacles. Particles can enter and leave this macroscopically occupied coherent quantum state to which we refer to as the condensate wave function. So the number of particles in this macroscopic state is not preserved and the actual condensate wave function depends on the scattering of particles into this state and the particles' decay. For example if no particles are added to the condensate, the occupation number of particles in this mode eventually goes to zero. On the other hand, due to pumping particles into the condensate patterns might emerge [38] and potentially new properties. While by naively applying the analogy to conserved BE condensed systems a superfluid phase aka absence of drag would be expected for small velocities below the speed of sound of the fluid due to particle interactions, while above a critical velocity a drag force would arise due to the possibility of emission of elementary excitations [1,2,30,31]. The continuous flow of particles into the condensate phase and the corresponding balancing drain may alter the scenario as studied numerically in [31] and thus superfluidity could be suppressed [33]. On the other hand could we scatter particles into the con-densate, such that the condensate wave function implies an effectively vanishing drag force aka superfluidity?

Physical scenario
So far the arrangement of pumping and the presence of decay has been considered in terms of analog behaviour to equilibrium BEC, e.g. the absence or reduction of scattering from an obstacle as it moves relativ to the condensate [52]. Here we particularly will consider an inserted obstacle to be co-moving with a particular pump distribution or equivalently both the pump and the obstacle are stationary and the surrounding Polariton quantumfluid is passing by. We will study the effect of arranging the pumping distribution properly without having additional momentum, when both are moving with the same relative speed to the condensate. This will allow us to observe in what way the obstacle can be affected by the non-equilibrium pumping distribution.
To elucidate those and superfluid aspects of an open interacting and coherent many-body system, let us introduce first the explicit theoretical framework we use as model of the condensate and subsequently we will derive by deduction key statements.

CONDENSATE WAVE EQUATION
The free energy of a BEC, which experiences the external potential V and which has a self-interaction strength g, is given by [7,47,50,51,53,54] Here we use the convention The dispersion of the condensate is ω(|k|) : R d → R and for simplicity can be assumed to be parabolic ω(|k|) ∼ k 2 (with |k| = k) corresponding to GP theory or free particles [58]. Naturally GP theory is considered in 3d, however, reduction to lower dimensions occurs e.g. for strong confinement along a spatial dimension [54]. In 3d g = 4π 2 a/m and max r V (r) = 4π 2 b/m ≡ V 0 are respectively the particle-particle and particle-impurity coupling, where a and b are the corresponding scattering lengths and m is the effective mass [1]. For the 2d case we have g 2d = √ 2π 2 a/a z and here a z = /mω z is the oscillator length with ω z being the trapping strength [7,54]. The condensate wave function ψ(r, t) is the minimiser of (1), i.e. the mode in which all particles of the dilute weakly-interacting closed Bose gas condense [7,54]. Mathematically ψ(r, t) : R d+1 → C and by performing a variation of the energy E[ψ] with respect to ψ * under the norm/particle preserving constraint ψ 2 2 = 1, we get the Euler-Lagrange equation for the minimiser. The so-called Gross-Pitaevskii equation is The chemical potential in Eq. 2, i.e. the energy needed to add another particle, is ∂E[ψ]/∂N = µ [54]. To account for the open non-equilibrium dynamics of the condensate (2) one generally considers an extension of the form which is applicable for atom laser systems [49] as well as for non-equilibrium Polariton condensates [30,47]. In (3) we consider additional physical parameters: Γ d is the homogeneous decay (or pump) rate of condensed particles, while P approximates the inhomogeneous nonequilibrium processes acting on the condensate fraction [47,49], which shall have an integrable Fourier transform. We assume that P (ψ, r, t) P (r, t)ψ(r, t), while we drop the superscript in what follows. Nonlinear density dependent processes can be regarded as adding a constant to the pump P , when considering a first order linear waves analysis. Examples of pump terms of this form are e.g. spatially dependent incoherent scattering of reservoir particles into the condensate phase [47]. These terms have been used to describe experimental results of Polariton condensates in the mean-field regime [55][56][57] and see [38,47] for theoretical works on that matter. Due to the generalisation of the guiding equation to include non-equilibrium processes the norm of the wave function ψ 2 2 = f (t) now varies in time.

IMPURITY WAVES
The impurity moving with velocity v in the stationary fluid frame is modelled by the external potential V (r, t) = V 0 e − 1 2σ 2 (r−vt) 2 [2], which has an amplitude V 0 and a width σ, or simply by a Dirac delta function V Dirac (r, t) = V 0 δ(r − vt) [1]. These potentials model an inserted atom or a laser beam [7,47,59]. We suppose the linear waves generated by these potentials are of the form ψ = φ 0 + δψ [1,2,60], where φ 0 represents the unperturbed part solving the complex GPE (3) without potential and δψ(r, t) denotes a small perturbation due to the presence of an impurity. By inserting this Ansatz in (3) and dropping terms of order δψ 2 we get the Bo-goliubov equation for this perturbation, Here we restrict our consideration to small weaklyinteracting impurities, V δψ. Furthermore we utilise the identity ∂δψ(r − vt)/∂t = −v∇δψ(r − vt) [1,2,60] and switch in the frame moving with the impurity, δψ(r, t) = δψ(r − vt) = δψ(r ), return to the notations without superscripts and in addition consider Eq. 4 in k-space. Thus the wavefunction δψ k = e −ikr δψdr satisfies the Bogoliubov equation in k-space given by We note the phase factor identity, µ = gn with n = |φ 0 | 2 and assume φ 0 to be real valued for the sake of simplicity.
Integrating the potential term e.g.
and so on and thus we get the algebraic equations of the form These are analytically solved, where we make use of the notation S 2d = f 2d (k 2 ) + iF(φ 0 P ). Similar linear waves have been derived in several scenarios and we refer to [1,2,58,60,61] for related results. Here the finite sized impurity is represented by a form factor f 2d (k 2 , σ, V 0 ) and the growth and decay terms alter the linear waves by extending the solution to the imaginary plane. The dimensional form factor is determined by the Fourier transform of the impurity and for general dimensions given by . Moreover for the simple impurity V Delta , the form factor becomes f Dirac = V 0 √ n. The energy spectra for the solutions (7) are in general complex-valued and thus a naiv application of Landau's critical velocity based on real-valued dispersions [62] cannot be applied, but we turn to the drag force as effective measure of friction.

DRAG FORCE
After deriving the form of the perturbed wave we investigate the force the impurity experiences as it moves through the non-equilibrium condensate. The definition of the drag force is [1,19] when assuming the mode of the Bose gas Ψ † is fully described by the complex order parameter equation (3). To calculate the drag force we employ the ansatz ψ = φ 0 +δψ and again neglect the δψ 2 terms. The result is Above calculation includes the special case of the Dirac delta, when V 0 ∼ 1/σ D as σ → 0, (see [1] for similar expressions). The remainder in (9) is calculated as explicitly presented in the appendix. So after some complex algebra we obtain the effective semiclassical force acting upon the impurity in the non-equilibrium condensate, In all cases and independent of dimension, we observe that the denominator of the drag force in Eq. 10 has a pole, if and only if This holds true for real-valued (free particle) dispersion relations ω(k) only when Γ d = 0, a case that corresponds to a BEC without homogeneous leakage, decay or constant gain of particles. The corresponding proofs of leading order existence of the superfluid phase transition for finite impurities and Delta impurities in equilibrium BEC was done in [2] and [1] respectively, where Sokhotsky's formula [63] was utilised to solve the integral with pole, thus implying the topological superfluid phase transition. The absence of a pole gives rise to a continuous (drag) force and correspondingly to the suppression of superfluidity when Γ d = 0 as observed numerically in [33] and analytically for a Dirac impurity in [61]. Physically particles leaving the condensed phase ψ add drag to the superfluid phase.
On the other hand here we point out that for a vanishing numerator in (10) or integrals thereof we observe a superfluid phase and an acceleration regime when the effective sign of the force projected on the direction of motion is positive. In Fig. 1 (a) we numerically show the velocity dependence of the drag force acting on the Gaussian impurity in 2d for various pumping strengths, given a pump of the form P = P 0 √ nδ(r − vt). While for low pumping strengths (triangles) the drag is persistent for all v we note that the absolute drag is decreasing again (when the sign is still negative) for larger velocities as previously observed in equilibrium condensates for finite sized impurities [2,59]. In contrast localised Dirac impurities yield a monotonically increasing amplitude of the drag force [1]. On the other hand for larger pumping strengths (circles and squares) we obtain forces with positive values, hence accelerating the impurity and when F v = 0 at specific v we effectively gain a superfluid phase due to balanced inhomogeneous and homogeneous non-equilibrium processes. Fig. 1 (b) shows the explicit (linear) dependence of the force on P for fixed velocities and thus the critical threshold of superfluidity and the impurity accelerating regime.
Considering a semiclassical picture the local gain of particles into the condensed phase can cause the impurity to be pushed further by the pumping induced density variation, thus opposing/annihilating the drag created by potentially exciting the superfluid. In this sense scattering from and into the condensed phase (of particles even without a momentum as assumed by the mathematical form of the pumping term) act as additional external forces on the impurity due to their spatial distribution close to the impurity. Those add to the forces due to exciting the fluid and due to quantum fluctuations. A higher local density of the condensate at the impurity implies the possibility of larger scattering rates between the impurity and the condensate, which therefore provide the possibility of effectively accelerating the impurity through scattering. By considering the simple mathematical scenario of a 1d Gaussian impurity potential and a step function density |ψ| 2 in (8), one can immediately conclude that adding density behind or in the wake of the impurity as it moves through the condensate will cause gain of momentum due to enhanced scattering probabilities in direction of motion. The scattered particles leave the condensate by gaining the momentum and thus induce via momentum conservation a gain of momentum of the impurity and effectively a force acting upon it.
Furthermore note that even an infinitesimal small addition to the superfluid density implies an infinitesimal acceleration of the impurity. Interestingly we note that for a locally driven condensate without homogeneous leakage or pumping, i.e. Γ d = 0, the impurity can be accelerated by this mechanism within a state of true superfluidity, i.e. below the speed of sound.
We point out that once the obstacle is accelerated by the pump, the pump has to be co-moving with the obstacle to maintain the acceleration of the obstacle. However in each stationary frame the pump does not induce additional momentum on the obstacle but is co-moving with the obstacle and merely modifies the density of the superfluid in the corresponding frame, which in turn accelerates the obstacle. In this sense the observed acceleration is a reversed drag force as discussed in [1]. Now to calculate the continuous force acting upon an impurity analytically, when the various contributions in the numerator (10) do not cancel out for the explanatory case, P = 0, we consider the projection of the drag force F on the velocity of the obstacle v, which we denote F v . By considering a parabolic kinetic dispersion ω and using the symmetry properties of the integrand, complex algebra and by approximating the integral by an expansion of the integrand to the quadratic order in v we obtain a simple expression for the integral. In 2d, which is the natural dimension e.g. for Polariton condensates, the result is a drag force acting in opposite direction of motion on the obstacle, Even for small velocities of the impurity the drag force obeys approximately quadratic behaviour in v, i.e. a nonvanishing drag force, in particular since the integral in (12) is strictly positive. This contrasts the equilibrium BEC results [1,2], where a critical velocity equal the speed of sound has been noted below which there is no drag at all corresponding to the superfluid phase. Logical consistency of the results presented here and those in [1,2] is given when Γ d → 0.
Next we solve the integral in (12) for the special case, Γ d µ, which yields the result where we have used the abbreviation µ ≡ µσ 2 and the notation Chi[z] = γ + log(z) + When neglecting the geometry of the impurity, i.e. V → V Delta , and by considering the quadratic order in v approximation, while again confining our consideration to 2d, we directly obtain the formula for the drag force, given d ≡ Γ 2 d − µ 2 > 0. The drag force is decreasing with increasing d approaching zero when d → ∞, hence suggesting an extremal superfluid regime for this special case, while finite d obey a quadratic in v dependence of its magnitude and when d = 0 we obtain the finite drag force F v = −c Dirac v 2 /(48π). The linear Schrödinger equation in the scenario of a quantum fluid flowing past an impenetrable cylindrical obstacle of radius R, obeys a drag law for high velocity or large object size, i.e. v /mR, that approaches the classical limit as well [3], i.e. F ideal v = −cρ 0 Rv 2 , where ρ 0 is the density of the fluid and c a dimension dependent constant. Although the results (12) and (14) show analog approximatively quadratic behaviour, the drag force increases less than quadratic for v /mR for the ideal Bose gas [3].

INHOMOGENEOUS PUMPING
Finally we turn to estimate the pumping terms in a 2d scenario. We suppose P = P 0 √ nδ(r − vt) and thus when switching in the moving frame we have F(P ) = P 0 √ n ∈ R, i.e. a Dirac function pumping spot moving with the impurity, while we note that the complementary extreme case, i.e. a homogeneous pumping spot would simply adapt Γ d in the previous considerations. For the sake of conciseness we consider a Dirac delta impurity and obtain in the quadratic order of v 2 an inhomogeneous force using d ≡ Γ 2 d − µ 2 = 0 with more details presented in the appendix. We note that the inhomogeneous force can be directed along the direction of motion and thus is opposing the drag created by the homogeneous pumping/decay terms for specific µ and Γ d , hence showing analytically in particular the possibility of a balance between homogeneous drag and local driving and thus a non-equilibrium superfluid phase which in the quadratic order order is independent of velocity. Furthermore increasing P 0 above this critical value implies an acceleration of the impurity that is proportional to P 0 .

DISCUSSION
By utilising Bogoliubov's perturbation theory applied to the scenario of a (finite sized) Gaussian and a Delta function impurity moving through a non-equilibrium Bose-Einstein condensate we have obtained the formulas for the drag forces acting opposite the direction of motion on the impurities and the opposing driving forces. In general movement with velocities larger than the speed of sound leads to a non-zero drag force due to Cherenkov radiation of phonons as previously noted for the Dirac and Gaussian impurities in [1] and [2] respectively. In stark contrast, however, we have observed that the drag force is non-vanishing as soon as the impurity is set in motion due to a constant non-equilibrium term in the governing equation, which confirms the numerics in [33]. Furthermore it has been shown that the force depends on the width and amplitude of the moving obstacle in the stated analytical form. We point out that the presented analysis does not include the effective drag due to nonlinear excitations such as vortices, vortex rings, solitary waves or solitons, which add energy dissipation and thus cause additional drag to the impurity. When nonlinear excitations are absent our mathematical analysis clarifies the linear waves contribution to the (drag) force acting on the weakly-interacting impurity as it moves at any velocity through the open BEC. Now, if a balance between homogeneous and inhomogeneous non-equilibrium contributions is given a superfluid phase, i.e. a vanishing drag force in the stationary reference frame is feasible even above the speed of sound as shown numerically and analytically. The most significant insight however is that the impurity can be accelerated by the inhomogeneous terms indicating a regime of the non-equilibrium condensate beyond that of superfluidity. The physical situation can be understood e.g. by a generalised Lagrangian formalism [39], which treats the complex non-equilibrium terms as external forces acting on the condensate wave function. By direct inspection of (8) one observes that for a Gaussian shaped impurity adding amplitude behind the moving impurity will imply an acceleration due to increased likelihood of scattering between the impurity and the condensate and the possible gain of momentum. The non-equilibrium terms, as presented here, can cause drag as well as an effective acceleration of the impurity when the pump is co-moving with the impurity. Finally using a more subtle approximation of quantum fluids by taking into account quantum fluctuations suggests that there may be additional dissipation [19] even in closed BEC. However, those fluctuations can -on average -be balanced or even revresed by inhomogeneously driven forces as shown in this paper.

Calculating the drag force without constant gain/loss
To calculate the drag force we rely on the following calculation. For the sake of simplicity we just consider the case in 2d here while presenting the formulas for the general dimensions in the main text. Using the definition of the drag force, Eq. 8 in the main text, and by switching into momentum space we get, by neglecting δψ 2 terms. Here we have defined F 0 in the last step. To determine the remainder in (16) for δψ, we use the following calculation: Further we calculate the 2d drag force for P = 0 by considering polar coordinates.
Next we assume a parabolic dispersion and use internal symmetries to simplify the integral, i.e.
Finally we expand the integrand in v to the quadratic order. The result is given by where we use the abbreviation cos θ ≡ x. Integration of x 2 / √ 1 − x 2 gives π/2 and for k we substitute k 2 = ρ and so obtain as presented in the main text.