Stokes drift and impurity transport in a quantum fluid

Stokes drift is a classical fluid effect in which travelling waves transfer momentum to tracers of the fluid, resulting in a non-zero drift velocity in the direction of the incoming wave. This effect is the driving mechanism allowing particles, i.e. impurities, to be transported by the flow; in a classical (viscous) fluid this happens usually due to the presence of viscous drag forces. Because of the eventual absence of viscosity in quantum fluids, impurities are driven by inertial effects and pressure gradients only. We present theoretical predictions of a Stokes drift analogous in quantum fluids for classical impurities obtained using multi-time analytical asymptotic expansions. We find that, at the leading order, the drift direction and amplitude depend on the initial impurity position with respect to the wave phase; at the second order, dominant after averaging over initial conditions, our theoretical model recovers the classical Stokes drift but with a coefficient that depends on the relative particle-fluid density ratio. Numerical simulations of a two-dimensional Gross-Piteaveskii equation coupled with a classical impurity corroborate our findings. Our predictions are experimentally testable, for instance, using fluids of light obtained in photorefractive crystals.

The transport properties caused by waves and, more in general, the interaction of waves with particles has been a long standing problem in many field of physics. A seminal example is given by particles floating under or on the surface of propagating (incompressible) water waves that experience a velocity drift which is known as the Stokes drift. This phenomenon, first described by G. G. Stokes in 1847 [1], is related to the intrinsic nonlinearity that characterises the Lagrangian description of a water particle immersed in a linear (or nonlinear) Eulerian wave field. Water particles move on trajectories that are not closed and, on average, advance in the direction of the propagation of the waves with a velocity that is one order of magnitude smaller that the phase velocity. This drift is important for the mass transfer of any object in a wave field specially in the area of sediment transport [2][3][4] and it is responsible of important fluid-mixing processes.
Acoustic waves do also transport particles: just like in an optical wave field, particles in an acoustic field are affected by the so-called acoustic radiation force that is the result of a transfer of linear momentum from waves to a particle [5][6][7][8]. The first calculation of this effect was reported in 1934 by L. King [9] showing that the radiation pressure is always in the direction of propagation of the wave. Since then, there has been a considerable interest in modelling the particle dynamics that results from averaging over many wave periods, especially because small particles may be used as tracers to visualise the flow.
The transport of particles in quantum fluids exhibiting both normal (viscous and thermal) and superfluid com-ponents has been investigated only in the last couple of decades [10,11]. Particles move due to the presence of pressure gradients and the drag force caused by the viscous (normal ) component. In superfluid liquid helium, density waves have very low amplitude due to the low compressibility of the system and pressure gradients origin mainly from the velocity field produced by quantised vortices. Particles become trapped into the core of quantised vortices [12,13] and have been used as tracers to probe, for example, the existence of quantised vortex filaments [14] and vortex reconnections [15]. Also, note that in the typical liquid helium experiments with particles the normal component cannot be neglected, hence particles also experience a classical Stokes drag.
In the limit of zero temperature, a quantum fluid has no normal component in the flow and particles are driven only by pressure gradients of the superfluid component as viscous dissipation is absent [13]. Differently to liquid helium, weakly interacting quantum fluids are highly compressible and both classical [16] and quantum [17] particles (usually called impurities in this context) can move thanks to the interaction with vortices [18][19][20] but also due to density waves [21]. Examples of such fluids are dilute gaseous Bose-Einstein condensates and quantum fluids of light, both of which can be quantitatively described using the Gross-Pitaevskii (GP) semiclassical model. In this Letter we present the Stokes drift and impurity transport in such setting.
For simplicity in the analytical predictions and numerical calculations, we consider a classical impurity (that is a classical-like particle having a well-defined position and momentum) whose characteristic size is of the order of the healing length of the system, and analyse how the acceleration caused by a superfluid density wave transports the impurity. The quantum fluid is described by a complex field ψ(x, t) and the impurity classical degrees of freedom, position and momentum, are q and p = M pq , respectively, given M p the impurity's mass. The dynamics of the system is governed by the (GP) equation coupled with a classical Newton's equation for the impurity; they read where m is the mass of the fundamental boson of the quantum fluid, µ is its chemical potential and g is the coupling constant of boson-boson local interaction. The potential V p (0) µ is localized around q and effectively determines the size of the impurity as its presence induces a complete depletion of the quantum fluid about the position q up to the characteristic distance a p where V p (a p ) = µ. The total energy of the system, the quantum fluid mass M = m |ψ| 2 dx and the total momentum P = i 2 (ψ∇ψ * − ψ * ∇ψ) dx + p are conserved quantities. This model has been successfully used to describe the interaction between impurities mediated by the superfluid [22] and their interaction with quantum vortices and Kelvin waves [13,20,23].
The system has a hydrodynamical interpretation via the Madelung transformation ψ(x) = ρ(x)/m e i m φ(x) that maps Eq.(1) into the continuity and Bernoulli equations of a fluid of density ρ and velocity v s = ∇φ. In absence of the impurity, the GP equation has a simple steady solution corresponding to the uniform state (condensate) |ψ 0 | = ρ 0 /m = µ/g. If Eq.(1) is linearised about ψ 0 , large wavelength waves propagate with the phonon (sound) velocity c = gρ 0 /m 2 and dispersive effects take place at length scales smaller than the healing length ξ = 2 /2gρ 0 . Using fluid dynamical variables, wave perturbations of the uniform solution along the x-direction, for instance, are simply (3) where we assume A ρ /ρ 0 1 and the angular frequency ω is the celebrated Bogoliubov dispersion relation Please refer to Section I of the Supplemental Material (SM) for detailed derivation. The wave period and wave length are thus defined as T = 2π/ω and λ = 2π/k, respectively.
We integrate numerically Eqs.(1-2) using the standard pseudo-spectral code FROST [24], and setting as initial condition the superposition of a linear wave solution, Eq.(3), with the ground state solution (obtained numerically by imaginary time evolution) of the impurity immersed in the quantum fluid. For simplicity we consider only two spatial dimensions with a double periodic rectangular domain of size 1024ξ × 128ξ, using 1024 × 128 collocation points. The potential used to model the impurity is a smoothed hat-function V p (r) = V 0 /2{1 − tanh (r 2 − η 2 a )/(4∆ 2 a ) }. Note that because of the nonlinearity/dispersion balance of the GP system, the quantum fluid density at the impurity boundary takes a distance ∼ ξ to heal to the bulk value. Thus, we can define an effective particle radiusā p > a p , estimated by measuring the volume of the displaced fluid πā 2 p = (|ψ 0 | 2 − |ψ p | 2 ) dx, where ψ p is the steady state with one impurity. We express the non-dimensional impurity mass as M = M p /M 0 , where M 0 = ρ 0 πā 2 p is the mass of the displaced fluid. In all the simulations, we fix the impurity potential V 0 = 20µ. We set the hard-core size to a p = 1.5ξ and the effective size toā p = 3.1ξ by choosing η a = ξ and ∆ a = 0.75ξ. For brevity, we will indicate the impurity position along the wave direction (x-axis) simply as q ≡ q x .
Initially, we place the impurity with zero velocity at different positions q 0 = q(t = 0) with respect to the phase of the incident wave. If the impurity is placed at the wave trough, we observe a drift in the same direction of propagation of the wave, see Fig.1.a); if it is placed a the crest, the impurity moves in the opposite direction, see Fig.1.b). Note that the figures show only a fraction of the numerical box, zoomed closed to the impurity. Our numerical observations are summarized in the sketch of Fig.2 where we define the initial impurity-wave phase as ϕ = q 0 k, with the convention that ϕ = 0 when q 0 is at the wave crest. In order to quantitatively characterize the drift effect, we monitor the impurity displacement as a function of time, for different initial impurity-wave phase ϕs. The impurity drift is displayed in Fig.3.a). In Fig. 3.b), we show the measured drift velocity v drift , computed averaging the impurity velocity over 13 wave periods, as a function of the impurity-wave phase and for two different wavelengths of the carrier wave. As the phase is changed from 0 to π, we observe a smooth transition from backward to forward drift.
To explain the change of direction of the impurity drift with respect to the impurity-wave phase, we build an effective minimal model to describe the problem. We start by considering that the force acting on the impurity, the right hand side of Eq.(2), is nothing but the fluid density gradient convoluted with the impurity potential. As formally derived in Section II of the SM, if the impurity size is much smaller than the wavelength λ and neglecting any active effect of the impurity onto the fluid, the   impurity dynamics is driven by the effective equation The small parameter is defined as and where we have introduced the added mass coefficient (C a = 1 in 2D) and two phenomenological dimensionless parameters γ 1 0.69 and γ 2 0.25 which account for the presence of a healing layer at the particle boundary; these values were obtained by fitting our theoretical prediction above using a small subset of simulations (see SM).
We want to establish the behaviour of the impurity position q at long times: introducing the slow time scale τ = t and using a standard multi-scale expansion, see Section III of the SM, we obtain the following expression for the drift velocity, i.e. the impurity velocity averaged over the fast timescale t   (7) at the leading order. d) Time evolution of the impurity rescaled position with drift parameter 2 for waves of different wavelengths and same initial impurity-wave phase ϕ = π/2; the second order prediction (9) is displayed in dashed line.
impurity-wave phases vanishes at order O( ). However, the next-to-leading order remains finite: It is interesting to notice that such result is equivalent to the Stokes drift in classical fluids for perfect tracers (see Section IV of the SM) v tracer drift = ω 2k (A ρ /ρ 0 ) 2 , a part from a coefficient that depends in this case on M.
In Fig.4 we further validate the model (5) and predictions (7) and (8) by varying its different parameters. and measuring the impurity displacement, initially set at ϕ = π, i.e, the wave trough, to highlight O( ) effects. First we set ϕ = π, i.e the impurity at the wave trough, to highlight O( ) effects: we consider waves of different wavelengths, Fig.4.a), and amplitudes, Fig.4.b), as well as impurity of different masses, Fig.4.c). In all the cases studied we observe that the motion curves of the particle collapse when the time is normalised by the wave period T and the displacement by λ . Then, in Fig.4.d, we also check the prediction for the drift at the order O( 2 ): since averaging over all the impurity-wave phases is computational demanding, we consider only the initial phase ϕ = π/2, for which the leading order vanishes. Equation which fits well the numerical data. Overall, we conclude that predictions (7) and (8) are robust up to O( 2 ). In summary, we have described and explained how a impurity immersed in a quantum fluid experiences a net transfer of momentum from an incoming density wave due to Stokes drift. At leading-order, the impurity drift depends on the initial impurity position with respect to the phase of the incoming wave: remarkably, it can move in any direction, independently on the direction of the incoming wave. When averaging over the initial impurity position with respect to the wave phase, this first-order effects cancel out and a non-vanishing second-order drift exits along the direction of the wave; this is consistent with the classical Stokes drift, but with a different coefficient.
The Stokes drift and impurity transport theoretical predictions reported in this Letter are derived under the assumption of a classical passive impurities, and were further corroborated by GP numerical results whit active impurities. This constitutes evidence of their quantitative applicability to quantum fluid experiments like superfluids of light made in photorefractive crystals. In a recent experiment [25], the dynamics of a very small defect (∼ 1ξ), obtained by imprinting a depletion in one side of the crystal, was studied to address the breakdown of superfluidity. Although the current propagation distance of the crystal remains relatively short (∼ 10ξ), it should be possible to observe the first-order drift effect. In order to observe the second order correction, a larger system would be needed; such requirement is challenging in current experiments but might be overcome on the future.
Finally, it is worth mentioning that, as manifested in Fig. 1, an interesting effect appears at long times. After about 13 linear periods T of interaction with the impurity (note that our system is periodic), the incoming wave deforms into a localised coherent structure that resembles a grey soliton. We expect the emergence of such nonlinear structures to be enhanced for larger incoming wave amplitudes. In this limit, the theory developed in this Letter might fail, which is consistent with the deviations observed in Fig. 4.b) for large amplitudes. The emergence and interaction of a soliton, or a train of solitons, and an impurity are interesting new lines of investigation that should be addressed in future works.
GK and MO acknowledge the support of the Simons Foundation Collaboration grant Wave Turbulence (Award ID 651471). GK was funded by the Agence Nationale de la Recherche through the project GIANTE ANR-18-CE30-0020-01. DP was supported by the EP-SRC First Grant Number EP/P023770/1. Computations were carried out at the Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur. This research was originally conceived to celebrate the bicentenary of the birth of Sir George Gabriel Stokes occurred the 13th August 1819. DP would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Dispersive hydrodynamics: mathematics, simulation and experiments when the final part of this work was undertaken, supported by EPSRC Grant Number EP/R014604/1.

I. QUANTUM FLUID DENSITY WAVES
We consider the Hamiltonian of a quantum fluid with an active classical impurity of mass M p immersed in it [1,2]: where q and p are the impurity position and momentum, while V p is the repulsive potential localized around q which models the impurity shape. In absence of impurities, the Hamiltonian models the Gross-Pitaevskii equation that reads: Here the uniform ground state is fixed by the chemical potential ψ 0 = µ/g. Given the Madelung transformation (x,t) , the corresponding hydrodynamic equations are where v = ∇φ. The pressure term is composed by a classical and a quantum term p = p cl + p q : We consider a slightly perturbed ground state density ψ w ( We consider a 1D perturbation along the x-axis. Linearizing around the ground state ρ 0 and separating real and imaginary parts, we get the linearized hydrodynamic equations: where v w = ∂φ/∂x. The system (6-7) admits a traveling density wave solution with density amplitude A ρ , provided that A ρ ρ 0 : with the dispersion relation where c = gρ 0 /m 2 is the speed of sound and ξ = / √ 2gρ 0 is the healing length.

II. EFFECTIVE MODEL FOR IMPURITY-WAVE INTERACTION
When the impurity is present, the ground state ψ p corresponds to a uniform condensate far from the impurity and a density depletion at places where |x − q| < a p , with a p the distance at which V p (a p ) = µ. At the boundary of the impurity the density goes from zero to ρ 0 in a distance of the order of the healing length ξ. The mass of the displaced quantum fluid is defined as from which we can define the effective size of the impurity asā p = M 0 /πρ 0 (in 2D, assuming polar symmetry). Figure 1 displays the quantum fluid density depletion generated by an impurity, the corresponding hard-core size a p and effective sizeā p . A good approximation when a p ξ is given by the Thomas-Fermi ground state that is obtained neglecting the kinetic term in the GP model. It reads with θ H the Heaviside function. If one neglects the healing layer, the impurity can be considered as a homogeneous disc, whose main effect is to exclude a ball B(q,ā p ) centered at the impurity position with radiusā p from the quantum fluid domain. In first approximation, given a density wave perturbation ρ w the presence of the impurity modifies the density as ρ = ρ p (ρ 0 + ρ w )/ρ 0 . By variating the Hamiltonian (1) with respect to q we obtain the equation of motion for the impurity which can be rearranged as where we used that 1 − Vp µ ρp ρ0 ∼ ρp ρ0 . Last expression can be further split into where the first integral vanishes assuming periodic boundary conditions. Finally we recognize the gradient of the pressure: where we assumed p ∼ p cl . Now we can apply the Bernoulli equation −p/ρ = ∂φ ∂t + 1 2 |∇φ| 2 , where the potential flow φ = φ w + φ p + φ BC takes into account the flow generated by the wave φ w and the flow due to the presence of the impurity φ p : The potential φ BC is in principle determined by the condition at the impurity boundary ∇φ BC ·n = [v w (q +ā p n) − v w (q)]· n. It can be neglected for small impurity (ka p 1). The equation of motion for a small impurity can be derived, using the standard procedure of impurity dynamics in classical fluids [3]. We treat separately the two contributions of the flow to the force F = F w + F p . The force generated by the wave flow is: where we introduced the phenomenological parameter γ 1 to account for the compressible effects due to the healing layer. The force generated by the flow φ p is computed using the Stokes theorem, assuming that the effective impurity boundary is placed at a distanceā p from the center. In 2D: where C a is the added mass coefficient (in 2D C a = 1) and γ 2 is another phenomenological adjustable parameter, due to the non-sharpness of the impurity boundary. The effective equation of motion for the impurity is then: Note that the contribution of the quantum pressure in (15)  k and for long waves ξk 1 it can be safely neglected. The values γ 1 0.69 and γ 2 0.25 used in the Letter were obtained by fitting subsets of the simulation data. Specifically, we first set γ 1 = 1 as for an impurity with sharp boundaries. Then, we obtained γ 2 = 0.25 by fitting the leading of our prediction using the measured drift velocity for the lightest impurity (M = 0.1), smallest wavelength (32ξ) and smallest amplitude (A = 0.01). Finally, once γ 2 was fixed, we made γ 1 free and determined its value by minimizing the root-mean-square error between the (normalized) displacement curves for the two extreme masses (M = 0.1 and M = 2 in Fig.4c of the main text).

III. MULTI-SCALE EXPANSION FOR THE EVALUATION OF THE DRIFT
By rescaling the variablesq = kq andt = ωt, equation (20) reads In the following we drop the tilde for simplicity in the notation. We use a multi-scale expansion to establish the behaviour of q(t) at long times rescaled by −1 . We look for a solution q(t, ) = Q(t, t, ), where Q(t, τ, ) is a function of the fast time variable t and the slow time variable τ that gives q(t, ) when τ = t. Equation (21) becomes: Plugging the expansion Q(t, τ, ) = q 0 (t, τ ) + q 1 (t, τ ) + 2 q 2 (t, τ ) + O( 3 ) into (22) we obtain the hierarchy: The order O( 0 ) implies q 0 (t, τ ) = C 00 (τ )t+C 01 (τ ), with C 00 = 0 in order to cancel the secular term in the fast variable t. The order O( 1 ) gives q 1 (t, τ ) = − sin (q 0 (τ ) − t) + C 10 (τ )t + C 11 (τ ), and also C 10 (τ ) must vanish. Substituting these results in the equation (25), we obtain: The solution q 2 (t, τ ) will not contain secular terms in t only if d 2 q 0 /dτ 2 = 0, which implies where B 1 and B 2 are constant at all time scales. Therefore, the solution at order O( 2 ) is Substituting this in the equation (26), we obtain: The solution q 3 (t, τ ) will not contain secular terms in t only if it is a periodic function of t. This implies that the average over t of each derivative ∂ n q3 ∂t n must vanish. Therefore, averaging Eq. (30) over t and imposing ∂ 2 q3 we have the condition d 2 C11 dτ 2 = B 1 , which implies where D 1 and D 2 are constant at all time scales. Therefore, the solution at order O( 1 ) is The velocity at order O( 2 ) is: The constants can be fixed from the initial condition. From Eq. (32) we have: From Eq. (33) and settingq(0) = 0 for simplicity we have: The drift velocity is finally obtained averaging (33) over the fast timescale v drift = q t = − cos(q(0)) + 2 1 + 1 4 cos (2q(0)) + O( 3 ).
Now we can come back to the original variables and define ϕ = kq(0) the the initial impurity-wave phase, with the convention that ϕ = 0 when q(0) is at the wave crest. Therefore, the drift becomes Notice that, even if the initial velocity of the impurity is zero, at the leading order there is still a drift proportional to the parameter that depends on the initial position of the impurity with respect to the wave. In particular, if the impurity is initially placed in a trough of the density wave, it will have a drift in the same verse of the wave, while if it is placed in the peak it will move backwards. If we average over random initial conditions, we see that there is still a net drift at the second order, that is independent on the initial position of the impurity: Note that the same effect is recovered placing the impurity at an initial phase ϕ = (2n + 1) π/2. The only difference is that the initial condition dependence of the second order does not vanish, resulting in a different prefactor: It is interesting that the phase averaged drift (40) resembles, up to a constant, to the drift v tracer drift obtained for a perfect tracer (see next section).

IV. STOKES DRIFT FOR A PERFECT TRACER IN A CLASSICAL VISCOUS FLUID
In a classical fluid, perfect (small) impurity tracers follow perfectly the flow under the assumption of infinite Stokes drag. The dynamic of such impurity tracer would then bė q = v w (q(t), t) = ε ω k cos (kq − ωt), ε = A ρ ρ 0 For completness, we derive here the Stokes drift in the case of a tracer in a classical fluid. Using the rescaled variablesq = kq andt = ωt, equation (42) reads: This can be treated again using a multi-scale approach in the small parameter ε. Dropping the tilde for simplicity, we look for a solution q(t, ε) = Q(t, εt, ε), where Q(t, τ, ε) is a function of the fast time variable t and the slow time variable τ that gives q(t, ε) when τ = εt. Equation (43) becomes: Plugging the expansionq(t, τ, ε) = q 0 (t, τ ) + εq 1 (t.τ ) + ε 2 q 2 (t, τ ) + O(ε 3 ) inside (22) we obtain the hierarchy: The order O( 0 ) implies q 0 (t, τ ) = q 0 (τ ). The order O( 1 ) gives q 1 (t, τ ) = − sin (q 0 (τ ) − t) + dq0 dτ t + C 1 (τ ), therefore dq0 dτ = 0 in order to cancel the secular term in the fast variable. We have q 0 constant at all timescales, and the order O( 2 ) becomes: The function q 2 (t, τ ) is periodic in t, and therefore without secular terms in t, only if ∂q2 ∂t t = 0. Averaging Eq. (48) over the fast variable we get: Now we can write the solution for the velocity up to the order O(ε 2 ) and obtain the Lagrangian drift velocity for the tracer that in the dimensional variables is