# Stokes drift and impurity transport in a quantum fluid

Umberto Giuriato, Giorgio Krstulovic, Miguel Onorato, Davide Proment

► To cite this version:

Umberto Giuriato, Giorgio Krstulovic, Miguel Onorato, Davide Proment. Stokes drift and impurity transport in a quantum fluid. Physical Review A, 2023, 107 (6), pp.L061303. 10.1103/Phys-RevA.107.L061303. hal-04304916

HAL Id: hal-04304916

https://hal.science/hal-04304916v1

Submitted on 15 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stokes drift and impurity transport in a quantum fluid

Umberto Giuriato and Giorgio Krstulovic

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France

Miguel Onorato

Dipartimento di Fisica, Università degli Studi di Torino and INFN, Sezione di Torino, Via Pietro Giuria 1, 10126 Torino, Italy

Davide Proment

School of Mathematics, University of East Anglia, Norwich Research Park, Norwich NR4 7TJ, United Kingdom; Centre for Photonics and Quantum Science, University of East Anglia, Norwich Research Park, Norwich NR4 7TJ, United Kingdom; and ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum fuer Schwerionenforschung, Planckstrasse 1, 64291 Darmstadt, Germany

<figure><img src="image_2.png" /></figure>

(Received 15 December 2022; revised 19 April 2023; accepted 17 May 2023; published 12 June 2023)

Stokes drift is a classical fluid effect in which traveling waves transfer momentum to tracers of the fluid, resulting in a nonzero 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 finding that, at the leading order, the drift direction and amplitude depend on the initial impurity position with respect to the wave phase, and at the second order, our theoretical model recovers the classical Stokes drift but with a coefficient that depends on the relative particle-fluid density ratio. Our theoretical predictions are obtained for classical impurities using multitime analytical asymptotic expansions. Numerical simulations of a two-dimensional Gross-Pitaevskii equation coupled with a classical impurity corroborate our findings. Our findings are experimentally testable, for instance, using fluids of light obtained in photorefractive crystals.

DOI: 10.1103/PhysRevA.107.L061303

The transport properties caused by waves and, more in general, the interaction of waves with particles has been a long-standing problem in many fields 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 Stokes in 1847 [1], is related to the intrinsic nonlinearity that characterizes 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 than the phase velocity. This drift is important for the mass transfer of any object in a wave field, especially in the area of sediment transport [2–4], and it is responsible for important fluid-mixing processes.

Acoustic waves 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–8]. The first

calculation of this effect was reported in 1934 by 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 modeling the particle dynamics that results from averaging over many wave periods, especially because small particles may be used as tracers to visualize the flow.

The transport of particles in quantum fluids exhibiting both normal (viscous and thermal) and superfluid components 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 originate mainly from the velocity field produced by quantized vortices. Particles become trapped in the core of quantized vortices [12,13] and have been used as tracers to probe, for example, the existence of quantized 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

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Rescaled quantum fluid density $\rho/\rho_0$ at different times in the presence of a density wave of wavelength $\lambda = 128\xi$ moving from left to right and an impurity (green dot) of relative mass $M = 0.1$. Note that the plots show only a portion of the entire (periodic) computational box. (a) Wave amplitude $A_\rho/\rho_0 = 0.01$ and the initial position of the impurity in the trough of the wave. (b) Wave amplitude $A_\rho/\rho_0 = 0.01$ and initial position of the impurity on the crest of the wave.</figcaption></figure>

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–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 analyze how the acceleration caused by a superfluid density wave transports the impurity. The quantum fluid is described by a complex field $\psi(\mathbf{x}, t)$ and the impurity classical degrees of freedom, position, and momentum, are $\mathbf{q}$ and $\mathbf{p} = M_p \dot{\mathbf{q}}$, 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

$$i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + (g|\psi|^2 - \mu)\psi + V_p(|\mathbf{x} - \mathbf{q}|)\psi, \quad (1)$$

$$M_p \ddot{\mathbf{q}} = - \int V_p(|\mathbf{x} - \mathbf{q}|) \nabla |\psi|^2 d\mathbf{x}, \quad (2)$$

where $m$ is the mass of the fundamental boson of the quantum fluid, $\mu$ is its chemical potential, and $g$ is the coupling constant of boson-boson local interaction. The potential $V_p(0) \gg \mu$ is localized around $\mathbf{q}$ and effectively determines the size of the impurity as its presence induces a complete depletion of the quantum fluid about the position $\mathbf{q}$ up to the characteristic distance $a_p$ where $V_p(a_p) = \mu$. The total energy of the system, the quantum fluid mass $M = m \int |\psi|^2 d\mathbf{x}$, and the total momentum $\mathbf{P} = \frac{i\hbar}{2} \int (\psi \nabla \psi^* - \psi^* \nabla \psi) d\mathbf{x} + \mathbf{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 $\psi(\mathbf{x}) = \sqrt{\rho(\mathbf{x})/m} e^{i(m/\hbar)\phi(\mathbf{x})}$ that maps Eq. (1) into the continuity and Bernoulli equations of a fluid of density $\rho$ and velocity $v_s = \nabla\phi$. In absence of the impurity, the GP equation has a simple steady solution corresponding to the uniform state (condensate) $|\psi_0| = \sqrt{\rho_0/m} = \sqrt{\mu/g}$. If Eq. (1) is linearized about $\psi_0$, large wavelength waves propagate with the phonon (sound) velocity $c = \sqrt{g\rho_0/m^2}$ and dispersive effects take place at length scales smaller than the healing length $\xi = \sqrt{\hbar^2/2g\rho_0}$. Using fluid dynamical variables, wave perturbations of the uniform solution along the $x$ direction, for instance, are simply

$$\rho = \rho_0 + A_\rho \cos(kx - \omega t), \quad v_w = \frac{\omega A_\rho}{k\rho_0} \cos(kx - \omega t), \quad (3)$$

where $v_w = \partial\phi/\partial x$, we assume $A_\rho/\rho_0 \ll 1$, and the angular frequency $\omega$ is the celebrated Bogoliubov dispersion relation

$$\omega = c|k|\sqrt{1 + \frac{\xi^2 k^2}{2}}. \quad (4)$$

Please refer to Sec. I of the Supplemental Material (SM) for the detailed derivation [24]. The wave period and wavelength are thus defined as $T = 2\pi/\omega$ and $\lambda = 2\pi/k$, respectively.

We integrate numerically Eqs. (1) and (2) using the standard pseudospectral code FROST [25], 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\xi \times 128\xi$, using $1024 \times 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 - \eta_a^2)/(4\Delta_a^2)]\}$. Note that because of the nonlinearity and dispersion balance of the GP system, the quantum fluid density at the impurity boundary takes a distance $\sim \xi$ to heal to the bulk value. Thus, we can define an effective particle radius $\bar{a}_p > a_p$, estimated by measuring the volume of the displaced fluid $\pi \bar{a}_p^2 = \int (|\psi|^2 - |\psi_p|^2) d\mathbf{x}$, where $\psi_p$ is the steady state with one impurity. We express the nondimensional impurity mass as $M = M_p/M_0$, where $M_0 = \rho_0 \pi \bar{a}_p^2$ is the mass of the displaced fluid. In all the simulations, we fix the impurity potential $V_0 = 20\mu$. We set the
<figure><img src="image_4.png" /><figcaption>FIG. 2. Sketch of the drift. At the leading order, the drift direction and magnitude depends on the initial impurity-wave phase $\varphi$.</figcaption></figure>

hard-core size to $a_p = 1.5\xi$ and the effective size to $\bar{a}_p = 3.1\xi$ by choosing $\eta_a = \xi$ and $\Delta_a = 0.75\xi$. 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 at 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 $\varphi = q_0k$, with the convention that $\varphi = 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 $\varphi$'s. The impurity drift is displayed in Fig. 3(a). In Fig. 3(b), we show the measured drift velocity $v_{\text{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 $\pi$, 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 Sec. II of the SM, if the impurity size is much smaller than the wavelength $\lambda$ and neglecting any active effect of the impurity onto the fluid, the impurity dynamics is driven by the effective equation

$$ \ddot{q} = \epsilon \frac{\omega^2}{k} \sin(kq - \omega t). \qquad (5) $$

The small parameter $\epsilon$ is defined as

$$ \epsilon = \eta \frac{A_{\rho}}{\rho_{0}}, \quad \text{with} \quad \eta = \left( \frac{\gamma_{2}C_{a} + \gamma_{1}}{\gamma_{2}C_{a} + M} \right), \qquad (6) $$

<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Time evolution of the impurity rescaled position (solid lines) for different impurity-wave phases. Dotted lines represent the drift prediction obtained in Eq. (7) at the leading order. (b) Rescaled drift versus the impurity-wave phases for waves of wavelength $\lambda = 128\xi$ (circles) and $32\xi$ (triangles); the dotted line is the prediction (7) at the leading order.</figcaption></figure>

and where we have introduced the added mass coefficient ($C_a = 1$ in two dimensions) and two phenomenological dimensionless parameters $\gamma_1 \simeq 0.69$ and $\gamma_2 \simeq 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 the SM).

We want to establish the behavior of the impurity position $q$ at long times: introducing the slow timescale $\tau = \epsilon t$ and using a standard multiscale expansion (see Sec. III of the SM), we obtain the following expression for the drift velocity, i.e., the impurity velocity averaged over the fast timescale $t$:

$$ v_{\text{drift}} = \langle \dot{q} \rangle_t = -\frac{\omega}{k}\epsilon \cos(\varphi) + \frac{\omega}{k}\epsilon^2 \left( 1 + \frac{1}{4} \cos(2\varphi) \right) + O(\epsilon^3). \qquad (7) $$

This theoretical prediction (dashed lines of Fig. 3) is in very good agreement with the data of GP numerical simulations. Note that the drift velocity averaged over the impurity-wave phases vanishes at order $O(\epsilon)$. However, the next-to-leading order remains finite:

$$ \langle v_{\text{drift}} \rangle_{\varphi} = \frac{\omega}{k} \epsilon^2 + O(\epsilon^3). \qquad (8) $$

It is interesting to notice that such result is equivalent to the Stokes drift in classical fluids for perfect tracers (see Sec. IV of the SM) $v_{\text{drift}}^{\text{tracer}} = \frac{\omega}{2k}(A_{\rho}/\rho_{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 $\varphi = \pi$, i.e., the wave trough, to highlight $O(\epsilon)$ effects. First we set $\varphi = \pi$, i.e., the impurity at the wave trough, to highlight $O(\epsilon)$ effects: we consider waves of different wavelengths [Fig. 4(a)] and amplitudes [Fig. 4(b)], as well as the 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 normalized by the wave period $T$ and the displacement by $\lambda\epsilon$. Then, in Fig. 4(d), we also check the prediction for the drift at the order $O(\epsilon^2)$: since averaging over all the impurity-wave phases is computationally demanding, we consider only the initial phase $\varphi = \pi/2$, for which the leading order vanishes. Equation (7) results in

$$v_{\text{drift}}(\varphi = \pi/2) = \frac{3\omega}{4k}\epsilon^2 + O(\epsilon^3), \quad (9)$$

which fits well with the numerical data. Overall, we conclude that predictions (7) and (8) are robust up to $O(\epsilon^2)$.

In summary, we have described and explained how an 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 of the direction of the incoming wave. When averaging over the initial impurity position with respect to the wave phase, this first-order effect cancels out and a nonvanishing 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 impurity, and were further corroborated by GP numerical results with 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 [26], 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 in 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 localized coherent structure that resembles a gray soliton. We expect the emergence of such nonlinear structures to be

<figure><img src="image_6.png" /><figcaption>FIG. 4. Time evolution of the impurity rescaled position with drift parameter $\epsilon$ for (a) waves of different wavelength, (b) waves of different amplitude, and (c) impurities of different mass. Dotted lines represent the drift prediction (7) at the leading order. (d) Time evolution of the impurity rescaled position with drift parameter $\epsilon^2$ for waves of different wavelengths and same initial impurity-wave phase $\varphi = \pi/2$; the second-order prediction (9) is displayed in a dashed line.</figcaption></figure>

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.

G.K. and M.O. acknowledge the support of the Simons Foundation Collaboration grant Wave Turbulence (Award ID 651471). G.K. was funded by the Agence Nationale de la Recherche through the project GIANTE ANR-18-CE30-0020-01. D.P. was supported by EPSRC First Grant No. EP/P023770/1. Computations were carried out at the MÉsocentre SIGAMM hosted at the Observatoire de la Côte d'Azur. D.P. acknowledges the support of the Université Côte d'Azur for funding his visit to the Laboratoire Lagrange via the Campagne Professeurs Invités (IFA) 2021-2022. D.P. would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program Dispersive hydrodynamics: Mathematics, simulation and experiments when the final part of this work was undertaken, supported by EPSRC Grant No. EP/R014604/1.

[1] G. G. Stokes, Trans. Cambridge Philos. Soc. 8, 441 (1847); Mathematics and Physics Papers (Cambridge University Press, Cambridge, 1880), Vol. 1, p. 197.

[2] M. S. Longuet-Higgins, Philos. Trans. R. Soc. London, Ser. A 245, 535 (1953).

[3] G. Besio, P. Blondeaux, M. Brocchini, and G. Vittori, J. Geophys. Res.: Oceans 109, C04018 (2004).

[4] F. Santamaria, G. Boffetta, M. M. Afonso, A. Mazzino, M. Onorato, and D. Pugliese, Europhys. Lett. 102, 14003 (2013).

[5] H. Bruus, Lab Chip 12, 1014 (2012).
[6] I. D. Toftul, K. Y. Bliokh, M. I. Petrov, and F. Nori, Phys. Rev. Lett. **123**, 183901 (2019).

[7] J.-H. Xie and J. Vanneste, Phys. Fluids **26**, 102001 (2014).

[8] J. Cleckler, S. Elghobashi, and F. Liu, Phys. Fluids **24**, 033301 (2012).

[9] L. V. King, Proc. R. Soc. London, Ser. A **147**, 212 (1934).

[10] D. R. Poole, C. F. Barenghi, Y. A. Sergeev, and W. F. Vinen, Phys. Rev. B **71**, 064514 (2005).

[11] J. I. Polanco and G. Krstulovic, Phys. Rev. Fluids **5**, 032601(R) (2020).

[12] Y. A. Sergeev and C. F. Barenghi, J. Low Temp. Phys. **157**, 429 (2009).

[13] U. Giuriato and G. Krstulovic, Sci. Rep. **9**, 4839 (2019).

[14] G. A. Williams and R. E. Packard, Phys. Rev. Lett. **33**, 280 (1974).

[15] G. P. Bewley, M. S. Paoletti, K. R. Sreenivasan, and D. P. Lathrop, Proc. Natl. Acad. Sci. USA **105**, 13707 (2008).

[16] T. Winiecki and C. S. Adams, Europhys. Lett. **52**, 257 (2000).

[17] R. C. Clark, Phys. Lett. **16**, 42 (1965).

[18] N. G. Berloff and P. H. Roberts, Phys. Rev. B **63**, 024510 (2000).

[19] A. Villois and H. Salman, Phys. Rev. B **97**, 094507 (2018).

[20] U. Giuriato, G. Krstulovic, and S. Nazarenko, Phys. Rev. Res. **2**, 023149 (2020).

[21] U. Giuriato and G. Krstulovic, Phys. Rev. B **103**, 024509 (2021).

[22] V. Shukla, M. Brachet, and R. Pandit, Phys. Rev. A **94**, 041602(R) (2016).

[23] U. Giuriato and G. Krstulovic, Phys. Rev. B **102**, 094508 (2020).

[24] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevA.107.L061303 for a detailed derivation of the quantum fluid density waves, the effective model for impurity-wave interaction, the multiscale expansion for the evaluation of the drift, and a comparison with the classical Stokes drift.

[25] G. Krstulovic, A theoretical description of vortex dynamics in superfluids. Kelvin waves, reconnections and particle-vortex interaction, Habilitation à diriger de recherches (HDR) thesis, Université Côte d'Azur, France, 2020.

[26] C. Michel, O. Boughdad, M. Albert, P.-E. Larre, and M. Bellec, Nat. Commun. **9**, 2108 (2018).