Superfluid Drain Vortex

Drain vortices are among the most common vortices observed in everyday life, yet their physics is complex due to the competition of vorticity’s transport and diffusion, and the presence of viscous layers and a free surface. Recently, it has become possible to study experimentally drain vortices in liquid helium II, a quantum fluid whose physics is characterised by the absence of viscosity and the quantisation of the circulation in the superfluid component. Using the Gross–Pitaevskii equation, we make a simple model of the problem which captures the essential physics ingredients, showing that the drain vortex of a pure superfluid consists of a bundle of vortex lines which, in the presence of a radial drain, twist, thus strengthening the axial flow into the drain.


I. INTRODUCTION
Quantised vorticity is a distinguishing property of superfluids, indeed hundreds of papers have been written on this subject since Vinen's detection of single quanta of circulation in superfluid helium [1].It is therefore remarkable that, until recently, very little attention has been dedicated in the superfluid context to drain vortices (also called suction vortices or bathtub vortices).Combining azimuthal motion around an axis with radial/axial inflow into a hole, drain vortices are familiar to everybody because they can be easily created in a kitchen or bathroom sink filled with water.Familiarity is not the same as physical understanding, however: boundaries and a free surface create subtle Ekman layers, drainpipe axial flows and upwellings which, even for water drain vortices, have been studied only recently [2,3].
We consider liquid helium ( 4 He) at temperature T < T λ where T λ = 2.17 K is the critical temperature at saturated vapour pressure.In this low temperature regime, liquid helium consists of two inter-penetrating fluid components, an inviscid superfluid and a viscous normal fluid, whose proportions strongly depend on T .The normal fluid fraction tends to one for T → T λ and to zero for T → 0; viceversa, the superfluid fraction tends to zero for T → T λ and to one for T → 0 (in practice, at T = 1 K the superfluid fraction is already more than 99%).
As mentioned, the key property of the superfluid component is that any vorticity is concentrated in thin vortex lines of fixed circulation κ = h/m where h is Planck's constant and m is the mass of one helium atom.Therefore we expect that, in the configuration of a drain vortex, the flow pattern of the normal fluid should be (in the first approximation) similar to that of water, including viscous layers and a continuous vorticity field which fills the system, whereas the vorticity of the superfluid should be confined to a central cluster of vortex lines.However, since the vortex lines scatter the thermal excitations (phonons and rotons) which make up the normal fluid, there should also be a mutual friction force between normal fluid and superfluid components, whose precise effects are difficult to guess (the friction depends on the density of the vortex lines and their velocity difference with respect to the normal fluid).We also expect that the flow pattern should depend on whether the drain vortex is created mechanically (e.g. using a propeller) or thermally (e.g. using a heater), as the induced superfluid and normal fluid velocities will be parallel or antiparallel respectively.Finally, unless the container is very large, details of any fluid reinjection into the system will be important.In summary, the superfluid drain vortex problem contains many physical ingredients which may combine in a nontrivial way, even before considering the presence of a free surface, which may deepen creating a funnelling drainpipe.
Experimentally, the problem has been tackled recently by Yano and collaborators [4][5][6].Using a rotor, they created a drain vortex in helium at T = 1.6 K (corresponding to a superfluid fraction of 83%).By measuring the attenuation of second sound, they showed that the drain vortex consists of a cluster of approximately 10 4 quantised vortex lines which accumulate in a narrow central region near the axis of symmetry above the drain hole.The experiment was followed up by numerical simulations [7] based on the Vortex Filament Model (VFM).These simulations determined the evolution of seeding vortex lines at nonzero temperatures in the presence of a prescribed normal fluid in the shape of a Rankine vortex with a superimposed constant axial flow into the drain hole.
In this work, we consider the superfluid drain vortex problem in its simplest form: a vortex cluster in the presence of a drain flow in a pure superfluid at T = 0, in the absence of a free surface.We model the problem using the Gross-Pitaevskii Equation (GPE) for a weaklyinteracting gas of bosons [9,10].The GPE is idealized for helium (which is a liquid, not a dilute gas), and, since its numerical solution requires the resolution of length scales smaller than the vortex core, our calculation is necessarily limited to small-scale vortex configurations.
Still, our model captures the essential physical ingredients, namely the dynamics of vortex lines in a nontrivial three-dimensional geometry which includes a draining superflow and the presence of boundaries, ingredients which were not accounted in the VFM approach [7].Our work is articulated as follows: Section (II) describes our model and the numerical methods employed and Section (III) presents our results.

II. MODEL A. Governing equations
It is convenient to write the GPE in dimensionless form.We use the healing length ξ = / √ gmn 0 and the speed of sound c = n 0 g/m as units of length and speed respectively, τ = ξ/c as unit of time, the density of a uniform condensate, n 0 , as the unit of density, and express the trapping potential in units of the chemical potential gn 0 , g being the interaction parameter, m the mass of one atom, and = h/(2π) the reduced Planck's constant.The resulting dimensionless GPE is where ψ(r, t), n(r, t) = |ψ(r, t)| 2 , V (r), r and t are the dimensionless wavefunction, density, trapping potential, position and time respectively (hereafter all quantities are meant to be dimensionless).During the evolution, Eq. (1) conserves the total number of atoms, N , in the volume V of the system, and the energy, E, given respectively by We solve the three-dimensional GPE numerically in the computational domain x min ≤ x ≤ x max , y min ≤ y ≤ y max , z min ≤ z ≤ z max using the 4 th order Runge-Kutta method in time and centred differences in space; the spatial discretization consists of N x , N y and N z discretization points in the x, y and z directions respectively.Convergence is checked by monitoring the conservation laws; for typical time evolutions described in the next section, N and E are conserved with relative errors (N − N 0 )/N 0 ≈ 10 −5 % and (E − E 0 )/E 0 ≈ 0.05 % where N 0 and E 0 are the initial values.

B. Geometry of the system
For a system whose geometry is described in Fig. 1, we set V (r) = V trap (r) in Eq. ( 1), where the trapping po-tential V trap (r) is equal to zero inside the following three regions: (i) a cylinder of radius r ad extending from z 1 to z 2 which represents the experimental cell; (ii) a drain hole in the shape of a truncated cone with top radius r ′ sink and bottom radius r sink extending from z min to z 1 ; (iii) an injection annulus of inner radius r ad and outer radius r in extending from z 3 to z 4 .Outside these three regions, we set V trap (r) = 10 (in units of the chemical potential).This value corresponds to a high potential barrier imposing ψ = 0 outside the regions (i) -(iii).In the experiment of Yano and collaborators [4] the injection annulus is at the bottom of the experimental cell, so here we present results for z 3 = z 1 (i.e. the injection annulus is at the bottom of the cylinder).

C. Drain and injection
To model the drain hole, we add a negative imaginary part to the potential V (r) in region (ii) (the truncated cone below the cylinder), setting where V sink (r) is equal to a positive constant V 0 in the drain hole and zero elsewhere.Between t and t + ∆t, this negative imaginary potential removes ∆N (t) atoms from the drain hole, depleting the density in that region.From Eq. (1) it is in fact possible to deduce the following continuity equation via the Madelung transformation [9]: where v = j/n is the velocity of particles and the current j is defined as follows: By integrating Eq. ( 6) over the whole volume, we obtain that the rate of loss of particles into the system due to drain potential is given by − dN dt = 2V 0 N sink (t), where N sink (t) is the number of atoms in the drain hole.The number of atoms ∆N (t) removed from the drain hole in time ∆t is hence, in the first approximation, given by ∆N (t) ≈ 2V 0 N sink (t)∆t.As the quantum fluid described by the GPE is barotropic (the pressure is proportional to the square of the density [9]), the density difference arising from the atoms removal in the drain creates a pressure difference which drives a flow towards the drain (see Section III).In order to conserve the total number of atoms of the system, we add into the injection annulus the same number of atoms ∆N (t) that we have removed from the drain hole.After time-stepping the wavefunction from ψ(r, t) to ψ(r, t+∆t), the naive approach would be to change the density in the injection annulus from n(r, t+∆t) to n ′ (r, t+∆t) = n(r, t+∆t)+∆N/V in , where V in is the volume of the injection annulus.However, we have found that the resulting small radial discontinuity of the density at the edge of the injection annulus (r = r ad ) tends to destabilize the solution.A more stable injection is obtained if the injected density profile is continuous.Therefore we set n ′ (r, t+ ∆t) = n(r, t+ ∆t)+ f (r)∆N (t), where the distribution function f (r) vanishes at r = r ad and r = r in , and its volume integral is normalized to one.We choose Our particle re-injection protocol can physically be interpreted as the inclusion of a source term σ in in the continuity equation Eq. ( 5), i.e.
where the production term, σ in (r, t), is nonzero only inside the injection annulus.Since V sink (r) is equal to a constant, V 0 , inside the drain hole and vanishes outside it, in order that the number of atoms in the system remains the same, σ in (r, t) must satisfy III. RESULTS

A. Ground state
To find the ground state of the system, we start from the Thomas-Fermi approximation, imposing that ψ T F (r) = 1 − V trap (r) if V trap (r) < 1 and zero otherwise.Without any drain flow or injection, we integrate Eq. ( 1) in imaginary time, replacing t with −it and we enforce at every time step the condition that the total number of atoms in the system does not change.As healing regions develop near the boundaries, the wavefunction settles down to the desired time-independent state which minimizes the energy.

B. Steady drain flow
We solve the GPE using the ground state as initial condition, imposing drain flow and injection as described in Section II C.After an initial transient, we obtain a steady drain flow into the drain hole which is exactly compensated by the injection of atoms in the annular region, so that the number of atoms in the system remains constantly equal to the initial atom number.
It is interesting to relate V 0 (the amplitude of the imaginary potential in Eq. ( 4) generating the drain flow) to the the quantity which is controlled in the experiment: the flow rate from the injection annulus into the drain hole, Q inj .To do this, we integrate Eq. ( 9) over the volume V s of the drain hole, and define the flow rate out of V s , where Σ s is the surface which encloses V s .Because of the box-trap boundary conditions, the only flow in and out of V s is across the top surface Σ s of the truncated cone, i.e.Q sink = Σs nv • dS.We find In the steady state The steady number of atoms in the drain hole N 0 sink is not constant but depends on the potential V 0 , implying that the relation between Q inj and V 0 , Eq. ( 13), is not linear.The dependence of N 0 sink on V 0 can be, at least qualitatively, understood employing the Bernoulli equation which can be derived [9] from Eq. ( 1) neglecting the quantum pressure effects: where v s , p s and v c , p c are velocities and pressures respectively in the centre of the drain hole just below the cylinder (where the velocity is v = −v s ẑ, ẑ being the unit vector in the z-direction and v s > 0) and in the cylinder far away from the drain hole (where |v| = v c ≈ 0).Employing the Bernoulli equation Eq. ( 14), we obtain where the density in the drain hole is approximately n s = N 0 sink /V s .In first approximation, the flux across the top surface Σ s of the drain hole is Q inj = n s v s A s , where A s is the area of Σ s .Using this last relation, the continuity equation Eq. ( 13) and the Bernoulli relation Eq. ( 15), we obtain the following expression for N 0 sink : leading to This non-linear relation between Q inj and V 0 can be observed in Fig. 2(a) where we plot Eq. ( 17) and its linear approximation at small V 0 .It is clear that our assumptions in deriving Eq. ( 17) are valid when V 0 is not too large compared to the chemical potential.This nonlinear effect is also visible in Fig. 2(b) where we plot the dependence on the potential V 0 of the averaged radial and axial components of the current, j r and j z respectively, flowing into the drain hole, averaged over a horizontal disk of radius r ′ sink at z = (z 3 + z 4 )/2.Fig. 3 shows the steady drain flow pattern in a larger geometry plotted on the xz plane (a) and on the xy plane (b).It is apparent that the radial and axial flows into the drain hole are confined to the bottom of the cylinder, the radial component being stronger than the axial one.This flow pattern is very different from the flow pattern of a viscous drain hole described in Ref. [2], which suffers viscous friction near all boundaries; our drain flow is instead inviscid and irrotational.The plot of the pressure distribution (see Fig. 4) confirms that the flow is driven by pressure gradients, as expected.

C. Drain vortex flow
In the next numerical experiment, we compute the time evolution of a lattice of vortex lines in the presence of a drain flow.We add N v = 10 parallel vortex lines aligned along the z-direction which form a small lattice (bundle) of initial radius r b 25.This vortex imprint- ing is done in a standard way by multiplying the ground state wavefunction with the approximate wavefunction of a vortex in a homogeneous condensate and let the system then relax in imaginary time.By integrating in time the GPE, we follow the evolution of the vortex bundle, which is shown in Fig. 5.We find that, since vortex lines are advected by the superflow (Helmholtz's theorem) and since near the bottom of the cylinder the inward radial flow is much stronger than near the top, the lower part of the bundle is sucked radially inwards, towards the center of the drain hole, while the top part is basically unaffected (this effect is also apparent in Fig. 6  radius at two different heights).The quantisation of the superfluid's angular momentum [9,11] implies that the angular velocity Ω(z, t) of the bundle at height z and time t is given by the relation Ω(z, t) = N v κ/[2πr 2 b (z, t)].This implies that the bottom parts of the vortex lines which are near the drain hole rotate at larger angular velocity Ω (as r b is smaller), moving hence "ahead" with respect to the top parts of the vortex lines close to the top of the cylinder: the vortex bundle becomes more and more twisted, as shown in Fig. 5 (b).As this twist develops, the projection of the total vortex length on a plane orthogonal to the cylinder axis, Λ ⊥ , increases, leading to an increase of the total length, Λ, as shown in Fig. 6 (b).It takes a finite time for this initial twist to propagate upwards (see Fig. (5)) carried by Kelvin waves on vortices.Once the Kelvin waves reach the top of the cylinder, the bundle's radius, r b , starts decreasing near the top as well, decreasing Λ ⊥ .However, given the enduring asymmetry of the radial flow (stronger close to the drain hole, negligible near the top of the cylinder), Kelvin waves persist on the vortices, leading Λ ⊥ to settle to a finite, non-zero value.The dynamics of this rotating, twisted bundle is characterised by vortex reconnections which occasionally scramble the vortex lines, leading to a moderately disordered vortex configuration which is still strongly polarised in the z direction.
An important consequence of the twisted nature of vortex lines is that it induces a downwards axial superflow into the drain hole, stronger than the axial superflow caused by the drain hole without vortex lines.Fig. 7  temporal evolution of the total length Λ of the vortex bundle (blue) and its orthogonal projection Λ ⊥ on a plane perpendicular to the cylinder's axis (yellow).The parameters are as in Fig. 3.
axial flow in the upper part of the bundle.Probably this is the mechanism which is responsible for the formation of a central drainpipe funnel when the helium has a top free surface, as observed by Yano and collaborators in their experiments [4].Further numerical simulations with N v = 25 vortices confirm the scenario which we have described.

IV. DISCUSSION AND CONCLUSION
We have used the Gross-Pitaevskii equation to model the superfluid drain vortex in its simplest form: at T = 0 (i.e.without any normal fluid and associated viscous effects) and without complications arising from a free surface.We have found that the superfluid drain vortex consists of a moderately disordered twisted bundle of quantised vortex lines.The twist is generated by the radial inflow into the drain, which is stronger near the drain which brings the vortex lines closer to each other.We also found that the twist of the vortex lines induces a strong central axial flow into the drain; in the presence of a free surface at the top, this axial flow is probably the origin of the funnel or drainpipe which has been observed [4].
The drain flow twisted bundle is similar to other twisted vortex states which have been observed in superconductors [12], and, more relevantly, in rotating 3 He [13] and large-scale 4 He vortex rings [14][15][16].Finally, by combining order (the strong axial polarization) with a controlled amount of disorder (induced by the radial drain flow), our results show that this bathtub vortex flow is a promising configuration to develop theoretical tools based on the HVBK equations and Vinen equation [17,18] to model problems of rotating quantum turbulence [19], ranging from liquid helium counterflow [20,21] to neutron star glitches [22,23] to atomic gases [24].
We are grateful to Sam Patrick and Silke Weinfurtner for discussions and acknowledge the financial support of UKRI grant Quantum simulators for fundamental physics (ST/T006900/1).LG acknowledges the support of Istituto Nazionale di Alta Matematica (INdAM).

FIG. 1 .
FIG. 1.The geometry of our numerical experiments.Top: schematic cross section on the y = 0 plane.The yellow and pink regions are the drain hole and the injection annulus, respectively.Bottom: schematic three-dimensional view for z3 = z1.

1 FIG. 4 .
FIG. 4. Steady drain flow pattern plotted on the xy plane at z = (z3 + z4)/2.The arrows indicate the direction of the current and the colours correspond to pressure values.Parameters as in Fig. 3.

6 FIG. 7 .
FIG. 7. Top (a): Steady drain flow pattern in the presence of vortex lines plotted on the xz plane for y = 0.The arrows indicate the direction of the current and the colours represent the current's magnitude.For clarity, only the flow in the cylinder is plotted, ignoring the drain hole and the injection annulus.Bottom (b): Steady drain flow pattern plotted on the xy plane at z = (z3 + z4)/2.The parameters are as in Fig. 3.