A new self-consistent approach of quantum turbulence in superfluid helium

We present the Fully cOUpled loCAl model of sUperfLuid Turbulence (FOUCAULT) that describes the dynamics of finite temperature superfluids. The superfluid component is described by the vortex filament method while the normal fluid is governed by a modified Navier–Stokes equation. The superfluid vortex lines and normal fluid components are fully coupled in a self-consistent manner by the friction force, which induces local disturbances in the normal fluid in the vicinity of vortex lines. The main focus of this work is the numerical scheme for distributing the friction force to the mesh points where the normal fluid is defined (stemming from recent advances in the study of the interaction between a classical viscous fluid and small active particles) and for evaluating the velocity of the normal fluid on the Lagrangian discretisation points along the vortex lines. In particular, we show that if this numerical scheme is not careful enough, spurious results may occur. The new scheme which we propose to overcome these difficulties is based on physical principles. Finally, we apply the new method to the problem of the motion of a superfluid vortex ring in a stationary normal fluid and in a turbulent normal fluid.


Introduction
Quantum turbulence [1][2][3] is the disordered motion of quantum fluids-fluids which are governed by the laws of quantum mechanics rather than classical physics. Examples of quantum fluids are atomic Bose-Einstein condensates [4,5], the low temperature liquid phases of helium isotopes 3 He and 4 He, polariton condensates, and the interior of neutron stars. The underlying physics is the condensation of atoms which obey Bose-Einstein statistics. In this work, we focus on the superfluid phase of 4 He, frequently referred to as helium-II or He-II for short. Helium II can be described as an intimate mixture of two fluid components [6,7]: the viscous normal fluid, which is similar to ordinary viscous fluids obeying the classical Navier-Stokes equation, and the inviscid superfluid, whose vorticity is confined to vortex lines of atomic thickness and quantised circulation (also referred to as quantised vortices or superfluid vortices). At nonzero temperatures, quantum turbulence thus manifests itself as a a e-mail: luca.galantucci@newcastle.ac.uk (corresponding author) disordered tangle of vortex lines interacting with either a laminar or a turbulent normal fluid depending on the circumstances.
Despite the two-fluid nature and the quantisation of superfluid vorticity, several surprising analogies between classical turbulence and quantum turbulence have been established in helium II in the last years. One example is the same temporal decay of vorticity [8] in unforced turbulence. Another example is the energy spectrum of forced homogeneous isotropic turbulence, which, according to experiments [9][10][11], numerical simulations [12][13][14][15] and theory [16,17], displays the same Kolmogorov energy spectrum (the distribution of kinetic energy over the length scales) which is observed in classical turbulence. However, at sufficiently small length scales (smaller than the average inter-vortex distance), or when superfluid and normal fluid are driven thermally in opposite directions [18,19], fundamental differences emerge between classical and quantum turbulence [20,21].
The comparison between quantum and classical turbulence is therefore a current topic of lively discussions in low-temperature physics community. Essential progress in the understanding of the underlying physics is provided by the recent development of new experimental techniques for the visualisation of superfluid helium flows, including micron-sized tracers, polymer particles [22], solid hydrogen/deuterium flakes [23,24] and, more recently, smaller (thus less intrusive) metastable helium molecules [25] which fluoresce when excited by a laser.
Alongside experiments, numerical simulations have played an important role in understanding the physics of quantum turbulence, from interpreting experimental data to proposing new experiments. However, most numerical simulations have determined the motion of vortex lines as a function of prescribed normal fluid profiles without taking into account the backreaction of the vortex lines [26] onto the normal fluid. Only a small number of investigations have addressed this back-reaction, but most efforts have suffered from some shortcomings: They were restricted to two-dimensional channels [27,28], had too limited numerical resolution to address fully developed turbulence [29], considered only simple vortex configurations [30,31], were limited to decaying turbulence [32]), or did not reach the statistically steady state which is necessary to make direct comparison to classical turbulence.
This work presents a fully coupled, three-dimensional numerical model which tackles all these shortcomings and contains innovative features concerning the numerical architecture and the physical modelling of the interaction between normal fluid and vortex lines. We shall refer to our model as FOUCAULT, Fully-cOUpled loCAl model of sUperfLuid Turbulence. Compared to past studies, our model consists of an efficiently parallelised pseudo-spectral code capable of distributing the calculation of the normal fluid velocity field and the temporal evolution of the superfluid vortex tangle amongst distinct computational cluster nodes. This feature allows the resolution of a wider range of flow length scales compared to the previous literature, spanning large quasi-classical length scales as well as smaller quantum length scales. This wide range is of fundamental importance for the numerical simulation of turbulent flows of helium II with realistic experimental parameters. Clearly, the importance of the dynamically self-consistent approach varies depending on the application. Our motivation is that we need a computational tool like FOUCAULT to find the answer. In certain applications, we expect our approach to be relevant, for example the recent experiments with micron-size, hydrogen tracer particles. The particles are affected by the normal fluid (mainly via the Stokes drag) and by the superfluid flow which is generated by the vortex lines (via inertial effects). If the vortex lines perturb the normal fluid locally, the motion of the tracers will change (as preliminary addressed very recently [33]), changing, for instance, the probability that they are trapped or not in the vortex cores. The correct interpretation of experimental data depends on having an appropriate model of this physics.
From the physical point of view, FOUCAULT includes a new approach to determine the friction force between superfluid vortices and the normal fluid, developed from progress in the study of classical creeping flows [34], but modified in order to avoid the unphysical dependence of the friction on the numerical discretisation on the vortex lines. In addition, the force between the vortices and the normal fluid (ideally a Dirac delta function centred on the vortex lines) is regularised exactly employing a method recently developed in classical turbulence [35,36] for the consistent modelling of the two-way coupling between a viscous fluid and small active particles. This method stems from the small-scale viscous diffusion of the normal fluid disturbances generated by the vortex motion, thus regularising the fluid response to vortex forcing in a physically consistent manner. Our fully coupled local model is therefore a significant improvement with respect to previous algorithms which employed arbitrary numerical procedures for the distribution of the vortex forcing on the Eulerian computational grid of the normal fluid [37].
The organisation of the paper is the following. In Sect. 2, we describe the equations of motion of the superfluid, the normal fluid and the vortex lines with emphasis on the different existing models for the calculation of the friction force. In Sect. 3, we outline our numerical method, including details on the Vortex Filament Method, the Navier-Stokes solver, interpolation schemes and the procedure employed for regularising the friction force consistently. Next, in Sect. 4, we show two applications of our model. Finally, in Sect. 5 we summarise the main points and outline future work.

Helium II
If the temperature of liquid helium is lowered below T λ = 2.17 K at saturated vapour pressure, a second order phase transition takes place to another liquid phase, known as helium II, which exhibits quantum-mechanical effects (unlike the high-temperature phase above T λ , called helium I, which can be adequately described as an ordinary Newtonian fluid). In the absence of vortex lines, the motion of helium II is well described by the two-fluid model of Landau and Tisza. This model describes helium II as the intimate mixture of two copenetrating fluid components [6,7,[38][39][40], which can be accelerated by temperature and/or pressure gradients: the viscous normal fluid and the inviscid superfluid. Each component is characterised by its own kinematic and thermodynamic fluid variables. In this description, the total density of helium II, ρ, is the sum of the partial densities ρ n and ρ s of the normal fluid and the superfluid: ρ = ρ n + ρ s (hereafter we use the subscripts 'n' and 's' to refer to normal fluid and superfluid components, respectively). While ρ is approximately independent of temperature for T < T λ , the superfluid and normal fluid densities are strongly temperature dependent: in the limit T → T λ , helium II becomes entirely normal (ρ n /ρ → 1), whereas in the limit T → 0 it becomes entirely superfluid (ρ s /ρ → 1). In practice, as the normal fluid density decreases quite rapidly with decreasing temperature, a helium sample at T ≤ 1 K is effectively a pure superfluid (ρ s /ρ ≥ 0.99).
Loosely speaking, the superfluid component corresponds to the quantum ground state of the system governed by a macroscopic complex wavefunction , and the normal fluid corresponds to thermal excitations (strictly speaking, this identification of the superfluid component with the condensate is not correct: helium II is a liquid of interacting bosons, not a weakly interacting gas). At temperatures of the most experimental interest, the mean free path of the thermal excitations is short enough that the normal fluid behaves like an ordinary (classical) viscous fluid with non-vanishing dynamic viscosity and entropy per unit mass, η and s respectively. In contrast, the superfluid component is inviscid and incapable of carrying entropy (hence heat). Physically, the existence of two distinct velocity fields v n and v s signifies that, locally, two simultaneous distinct movements are possible, even though individual helium atoms cannot be separated into two components.
However, if helium II rotates, or if the relative speed v ns = |v n − v s | between normal fluid and superfluid exceeds a critical value, superfluid vortex lines nucleate, coupling the two fluids [39,40] in a way that invalidates the model of Landau and Tisza. It is instructive to write the superfluid wavefunction in terms of its amplitude and phase, = | |e iφ , and use the quantum mechanical prescription relating the velocity to the gradient of the phase: v s = (h/m)∇φ. One finds [5] that a vortex line is a topological defect: a one-dimensional region in three-dimensional space where the phase φ is undefined, and the magnitude | | vanishes exactly. Vortex lines are thus nodal lines of the wavefunction. The single valuedness of the wavefunction implies that the circulation integral Γ following a closed path C is either zero (if C does not encircle the vortex line), or takes the fixed value if C encircles the vortex line, where κ = h/m is the quantum of circulation. A multiplycharged vortex line carrying more than one quantum of circulation is usually unstable, breaking into many singly charged vortex lines. Equation (1) also implies that the superfluid velocity around a straight vortex line has the form v s = κ/(2πr ), where r is the distance to the vortex axis. Since the superfluid component has zero viscosity, this azimuthal velocity field persists forever. Notice that the corresponding superfluid vorticity is a Dirac delta function defined on the vortex axis. Physically, we can think of a vortex line as a thin tubular hole (centred on the vortex axis) around which the circulation has the fixed value κ. The radius of the hole is also fixed by quantum mechanics and has the value a 0 ≈ 10 −10 m, which is the characteristic distance over which the amplitude of drops from its bulk value (at infinity, away from the vortex) to zero (on the vortex axis).
Vortex lines act as scattering centres for the thermal excitations (phonons and rotons) constituting the normal fluid [41,42]. This interaction produces an exchange of momentum, hence a mutual friction force F ns , between the superfluid component and the normal fluid component. The quantisation of the superfluid circulation and the friction force between superfluid and normal fluid make helium II a more complex, richer system than the original two-fluid scenario of Landau and Tisza.
This aim of this work is to develop a suitable numerical framework to address simultaneously the coupled temporal evolution of (i) the normal fluid velocity field and (ii) the superfluid vortex lines, taking the friction into account. In the following Sects. 2.2-2.5, we outline the equations of motions of the fluid components, the superfluid vortex lines and the expression of the friction force employed in our numerical algorithm.

The normal fluid
The normal fluid is a gas of elementary excitations called phonons and rotons. In the range of temperatures which we are interested in, namely 1.5 K < T < T λ , the normal fluid density is mostly due to rotons [42]. The roton mean free path is [43] is the roton group velocity, μ = 0.16m 4 is the roton effective mass [44] and m 4 = 6.65 × 10 −27 kg is the helium mass; therefore, λ m f p varies from 13 × 10 −10 m at T = 1.5K to 2.2 × 10 −10 m at T = 2.15 K. These values are much smaller than the typical intervortex distance, ≈ 10 −6 m to 10 −4 m, at the flow's small scales, and the current experimental facilities at the flow's large scales, which typically range from 10 −3 m to 10 0 m [33,45]. The normal component can hence be effectively described as a fluid with its own velocity v n , density ρ n , entropy per unit mass s and dynamic viscosity η. Given the small temperature gradients observed in experiments [46], ρ n , s and η can be treated as uniform and constant properties of the fluid. Furthermore, the normal fluid can be adequately treated as a Newtonian fluid, i.e. the viscous stress tensor is linear in velocity gradients. Stemming from these physical characteristics of the normal component, a long series of studies have focused on the derivation of the equations of motion of helium II in the presence of vortex lines, at length scales Δ much larger than the average inter-vortex spacing [42,47,48]: these equations of motion are nowadays referred to as the Hall-Vinen-Bekarevich-Khalatnikov (HVBK) equations. In the cited circumstance where the normal fluid undergoes incompressible isoentropic motion with constant and uniform dynamic viscosity, the HVBK equations of motion for the normal component in presence of vortex lines coincide with the Navier-Stokes equations for a classical incompressible viscous fluid with the addition of an extra term F ns representing the friction force: where p is pressure. Considering hereafter only isothermal helium II and introducing characteristic units of length and time, respectively λ and τ , Eqs. (2) and (3) can be made nondimensional as follows: where · indicates non-dimensional quantities and ν n = η/ρ n is the kinematic viscosity of the normal fluid. The numerical schemes for the integration of Eqs. (4) and (5) and the numerical handling of the friction force F ns are discussed in Sects. 3.2 and 3.4, respectively.

The superfluid
The complete set of HVBK equations consists of Eqs. (2) and (3) supplemented with the equations of motion of the superfluid component which, in the incompressible approximation and neglecting the tension force, are as follows: ∇ · v s = 0.
As already said, the superfluid vorticity ω s = ∇ × v s is confined to the vortex lines which can effectively be described as one-dimensional objects because the vortex core size a 0 ≈ 10 −10 m is several orders of magnitude smaller than the flow length scales probed in the present work (≈ 10 −5 m to 10 −6 m). The vortex lines can hence be treated as parametrised space curves s(ξ, t) in a three-dimensional domain, where ξ is arc-length and t is time. Within this mathematical description of vortex lines, the superfluid vorticity ω s can be expressed in terms of δ-distributions, as follows where is the unit tangent vector to the curve s(ξ, t) and L indicates the whole vortex configuration. The superfluid velocity may hence be expressed in terms of the spatial distribution of vortex lines s(ξ, t) via the Biot-Savart integral [49] where ∇φ(x, t) is the potential flow arising from the macroscopic boundary conditions. The corresponding non-dimensional equation is straightforwardly obtained and reads as follows: The numerical computation of Eq. (11) is addressed in Sect. 3.1.

The friction force
Historically, three distinct theoretical frameworks have been employed in the literature to model the interaction between vortex lines and normal fluid: (i) the coarse-grained philosophy, (ii) the local approach, and (iii) the fully-coupled local approach. The first approach probes the flow at length scales Δ much larger than the average inter-vortex spacing ; it was originally derived in the pioneering studies of Hall and Vinen [41,42] and successively led to the HVBK theoretical formulation [47,48]. The advantages of such approach are that it is dynamically self-consistent, and it makes possible to cover a wide range of length scales as for classical DNS of turbulence; the drawbacks are that it does not give information at length scales smaller than the inter-vortex distance, and, by relating the vortex line density to the magnitude of the superfluid vorticity, assumes perfect polarisation of the vortex lines, ignoring unpolarised vortex lines which contribute to the dissipation but not to the energy. The second approach addresses helium II dynamics at the scale of individual vortex line elements of length δ < ; it was first employed by Schwarz [49][50][51] and later reformulated by Kivotides, Barenghi and Samuels [30]. The benefits of this approach are that it models both polarised and unpolarised vortex lines equally well, but the results strongly depend on the imposed profile of the normal fluid, which is poorly determined even in the experiments. The third approach is a modification of the second that includes the effects of the local backreaction of the vortex lines on the normal fluid. The main advantage of the third approach is that it is dynamically self consistent (it does not require a priori information about the normal fluid) but has the same computational disadvantage of the second model: the evaluation of Biot-Savart integrals is computationally expensive. In addition, this third model also requires the numerical integration of the Navier Stokes equations for the normal fluid.
(i) Coarse-grained friction At length scales Δ , the discrete, singular nature of the superfluid vorticity field is lost. As a result, the averaged superfluid velocity and vorticity fields, indicated hereafter by v s and ω s respectively, are continuous fields, smoothly varying on the macro-scale Δ. At these scales, Hall and Vinen deduced the following expression for the friction force F ns acting per unit volume on the normal fluid in a bucket of rotating helium II: where v n is the coarse-grained normal fluid velocity and B and B are friction coefficients determined experimentally at scales Δ via second sound attenuation measurements in a rotating cryostat [41,42]. Equation (12) was generalised for non-straight vortex configurations by Bekarevich and Khalatnikov [47]. In this form (which contains a vortex tension term), the HVBK equations were applied to Couette flow [52][53][54], predicting with success the transition from Couette flow to Taylor vortex flow and the weakly nonlinear regime at higher velocities beyond the transition. Being laminar, these flows satisfy the assumption behind the derivation of the HVBK equations that the vortex lines are locally polarised. The application of the HVBK equations to turbulent flows is less straightforward, as the vortex lines are polarised only partially (most vortex lines are at random directions with respect to each other); hence, the relation | ω s | = κ L between the vortex line density L and the coarse-grained superfluid vorticity is only approximate. Nevertheless, the HVBK equations have been used to model the turbulent cascade in helium II [55].

(ii) Local friction
Developing the vortex filament method, Schwarz derived an expression for the force per unit length f sn exerted by the normal fluid on a single vortex line element of length δ < , position s(ξ, t) and unit tangent vector s [49]: . The field v n is thus determined by the macroscopic boundary conditions of the flow which is investigated and is entirely decoupled from the evolution of the superfluid vortex tangle. As a consequence, in the framework elaborated by Schwarz, the normal fluid velocity field may be prescribed a priori, depending on the fluid dynamic characteristics of the system studied: the explicit presence of the prescribed flow V n in Eq. (13) underlies this aspect. The second assumption, strongly linked to the previous, is the use of macroscopic coefficients α and α for the calculation of the friction acting on a single vortex element: α and α are in fact simply redefinitions of the friction coefficients B and B determined at large scales Δ in a very particular vortex line configuration (a lattice of straight vortex lines in a rotating bucket) [41,42]. Eq. (13), on the contrary, is intended to describe the force experienced by a vortex line element in a turbulent tangle. The measured values of α and α are displayed in Fig. 1. This decoupling of the normal fluid flow from the vortex lines motion (and, hence, the possibility of imposing arbitrarily a priori the normal fluid velocity field V n (x, t) felt by the vortices) confers to the theoretical framework pioneered by Schwarz a kinematic character: the evolution of the vortex tangle is determined for a given imposed normal flow. This local kinematic approach has been extensively employed in past studies to shed light on fundamental aspects of superfluid turbulence. In particular, various models of imposed normal flow V n have been studied: uniform [51,[56][57][58], parabolic [59][60][61][62][63], Hagen-Poiseuille and tail-flattened flows [64], vortex tubes [65], ABC flows [66], frozen normal fluid vortex tangles [67], random waves [58], time-frozen snapshots of turbulent solutions of Navier-Stokes equations [58,60,63] and time-dependent homogeneous and isotropic turbulent solutions of linearly forced Navier-Stokes equations [68]. (iii) Self-consistent local friction The self-consistent local approach was formulated in order to take into account the back reaction of the vortex lines onto the normal fluid [30,31,69], by modelling the dragging of the roton gas of excitations constituting the normal fluid by the vortex lines moving relative to it. Here the velocity field varies at length scales smaller than the average inter-vortex spacing and depends on the evolution of the vortex tangle: it cannot be prescribed, but has to be determined via the integration of the incompressible Navier-Stokes equations (2) and (3). This mismatch between the local normal fluid velocity v n perturbed by superfluid vortices and the macroscopic flow V n implies the necessity of re-determining the friction coefficients at these scales, now related to the local cross sections between the roton gas and the vortex lines [69]. This self-consistent local approach was first employed to investigate the normal fluid velocity field induced by simple vortex configurations, namely vortex rings [30] and vortex lines [31]; later, it was also used to study both the forcing of the normal fluid flow by decaying superfluid vortex tangles [70,71] and, vice versa, the stretching of an initially small superfluid vortex tangles by either a decaying turbulent normal fluid [32,37,72] or decaying normal fluid structures such as rings [73] and Hopf links [34].
The present work employs the most recent formulation of the third approach [34], slightly revisited. We model each vortex line element as a cylinder of radius a 0 and length δ (coinciding with the spatial discretisation of vortex lines, see Sect. 3.1 for details), where δ a 0 . The very small vortex core size a 0 implies that when a superfluid vortex is in relative motion with respect to the normal fluid, it generates a low Reynolds number normal flow around itself. Typical values of the Reynolds number associated to this vortex-generated normal flow in helium experiments are 10 −5 to 10 −4 . From the theory of low Reynolds number flows [74][75][76][77], the force per unit length f D which the normal fluid exerts on the vortex line is where the coefficient D (not to be confused with the friction coefficient discussed in Ref. [43]) is and γ = 0.5772 is the Euler-Mascheroni constant. Furthermore,ṡ = ∂s(ξ, t)/∂t is the velocity of the vortex line (see next paragraph 2.5), v n is evaluated on the vortex line element, that is to say v n = v n (s(ξ, t), t) using the interpolation schemes described and tested in Sect. 3.3, and the quantity v n ⊥ indicates the component of the normal fluid velocity lying on a plane orthogonal to s .
The expression (14) for the viscous force f D acting on the vortex line differs from the most recent approach [34] as it is independent of the discretisation δ on the lines. Including further the Iordanskii force f I = − ρ n κ s × (v n −ṡ) [78,79], the total force per unit length f sn acting on the superfluid vortices stemming from the interaction with the normal fluid is as follows In conclusion of this section, for the sake of completeness, it is important to mention that also Eq. (13) has been employed to take into account the back reaction of the superfluid vortices on the normal fluid, by averaging f sn on a normal fluid grid cell containing many vortex line elements [27][28][29].

The motion of vortex lines
The derivation of the equations of motion of the vortex lines is straightforward once Eq. (16) is taken into account. Since the vortex core is much smaller then any other scales of the flow, the vortex inertia can be neglected. As a consequence, the sum of all the forces acting on a vortex vanishes. Since the vortex line is effectively like a small cylinder in an inviscid fluid (the superfluid) surrounded by a circulation and in relative motion with respect to a background flow, it suffers a Magnus force per unit length which is where Setting the sum of all forces acting on the vortex line equal to zero, f sn + f M = 0, and employing Eqs. (16) and (17), we have Assuming that each vortex line element moves orthogonally to its unit tangent vector, i.e. s · s = 0, Eq. (18) leads [34,69] to the following equation of motion forṡ(ξ, t) where From the physical point of view, the motion of the vortices is hence governed only by temperature, which determines ρ n /ρ s and ν n /κ.

The self-consistent model
Our model of fully self-consistent motion of helium II at finite temperature comes from gathering the dimensionless equations (4), (5), (11), (16) and (19). The dimensional scaling The dimensionless viscosity ν 0 n is set to properly resolve the small scales of the normal fluid. Note that Γ is dimensionless and depends on temperature. When comparing the simulations to experiments, physical time and length scales are recovered from Eq. (24) by choosing the unit of length of the system, λ. With this notation, our self-consistent model, written in dimensionless form (after dropping tildes), reads: with In Eq. (25), we also added the external force F ext in order to sustain turbulence in the normal fluid. We recall that the physical parameters β, β and Γ depend only on temperature. Their behaviour is displayed on Fig. 1 using the tabulated values of ρ n , ρ s and ν n from reference [80]. For the sake of completeness, in Fig. 1 we also report the temperature dependence of the mutual friction coefficients α and α introduced by Schwarz [49] while modelling locally the mutual friction force in presence of a prescribed normal fluid flow (cfr. Eq. (13)). Finally note that for a given temperature, there are two more dimensionless parameters where v rms n is the root-mean-square or the characteristic normal fluid velocity fluctuation and L is its integral scale. The Reynolds number Re tells how turbulent is the normal fluid and the turbulent intensity I turb measures the strength of thermal counterflows with respect to normal fluid turbulent fluctuations. Note that, alternatively, we can also use a Reynolds number Re λ = v rms n λ T ν n based on the Taylor microscale of the flow λ T , which provides a precise definition in terms of velocity gradients and fluctuations [81] (the relation between the Taylor micro scale λ T and the turbulent kinetic energy dissipation in homogeneous and isotropic turbulence is as follows = 15ν n (v rms n /λ T ) 2 leading to Re λ ∼ √ 15 Re 1/2 ).

Numerical method
The numerical integration of the self-consistent model Eqs. (25)-(29) demands special care.
In particular, the coupling between the vortex filament method and the Navier-Stokes equation requires interpolations to determine the values of the normal fluid at the filament positions and a redistribution of the force f ns onto the mesh where the normal fluid is defined. As we shall see at the end of this section, a careless treatment of the coupling may lead to spurious numerical artefacts.
In this section, we describe the numerical method used to solve Eqs. (25)- (29) and provide physical and numerical justifications for the choice of the numerical scheme used to couple superfluid and normal fluid. We start by briefly describing the standard numerical schemes employed for the Vortex Filament Method and for the Navier-Stokes equations and then proceed to illustrate the novel friction coupling based on the method recently developed for the modelling of the two-way coupling between a classical viscous fluid and small active particles.

Vortex filament method
Here we briefly describe the Vortex Filament Method (VFM) to determine the time evolution of vortex lines. For a more in-depth overview of the VFM, we point the reader to a recent review article [26] and references within. The underlying assumption of the VFM is that vortex lines in the superfluid component can be considered one-dimensional space curves around which the circulation is one quantum of circulation κ. This is a reasonable assumption in helium II as the vortex core size a 0 is much smaller than any other characteristic length scale of the flow. However, during the time evolution of a turbulent vortex tangle, there are important events such as vortex reconnections, during which this assumption breaks down. We shall discuss how these events can be accounted for in the VFM at the end of the section.

(i) Lagrangian discretisation
At each time, we discretise the vortex tangle L in N p Lagrangian vortex points ..,N p the distance between neighbouring points is kept in the range [δ/2, δ] by removing or inserting additional points [82]. As the vortex configuration evolves, so does this Lagrangian discretisation -N p depends on the total length of the tangle. The vortex points evolve in time according to Eq. (19) written in dimensionless form, namelẏ where spatial derivatives along the vortex lines are performed employing fourth-order finite difference schemes which account for the varying mesh size along the vortex lines [83], and time integration is performed using the third-order Adams-Bashforth method.
For a detailed description of the standard numerical schemes employed and the illustration of the corresponding discretised equations, we refer to previous work in literature [82]. The interpolation of the normal fluid velocity v n on each vortex point s i presents some numerical issues; the different schemes which we have tested are outlined and discussed in Sect. 3.3. The evaluation of the superfluid velocity v s on the vortex points s i via Eq. (11) must be dealt with caution as the Biot-Savart integral diverges when x → s i in Eq. (11). This singularity is removed in a standard fashion by splitting the Biot-Savart integral into local and nonlocal contributions [50], that is (omitting time dependency to ease notation) by writing . In order to speed-up the calculation of Biot Savart integrals for dense tangles of vortices, Eq. (33) is approximated using a tree algorithm [84] which scales as N p log N p rather than N 2 p . (ii) Vortex reconnections In the limit of zero temperature, Eq. (10) expresses the superfluid's underlying incompressible, inviscid Euler dynamics in integral form [85]. Hence, on the basis of Helmholtz theorem, the topology of the superflow ought to be frozen, i.e. reconnections of vortex lines are not envisaged. We know however from experiments [86,87] and from more microscopic models [88][89][90][91][92][93] that when vortex lines come sufficiently close to each other, they reconnect, exchanging strands, as envisaged by Feynman [94]. In order to model vortex reconnections within the VFM, we supplement Eq. (32) with an ad-hoc algorithmical reconnection procedure. This strategy was originally proposed by Schwarz [49], and since then a number of alternative algorithms have been proposed. Whilst this procedure is essentially arbitrary, a number of different implementations have been extensively tested and compared [95], finding no significant physical difference between these implementations.

Navier Stokes solver
We solve Eq. (25) using a standard pseudo-spectral code in a three dimensional periodic domain of size L x × L y × L z with n x , n y and n z collocation points in each direction. Derivatives of the fields are directly computed in spectral space, whereas nonlinear terms are evaluated in physical space. The code is de-aliased using the standard 2/3-rule [96], and therefore, the maximum wavenumber is k max = 2π min [n x /L x , n y /L y , n z /L z ]/3. The main advantage of pseudo-spectral codes is that nonlinear partial differential equations are solved with spectral accuracy; this means that spatial approximation errors decrease exponentially with the number of collocation points. The drawbacks of pseudo-spectral codes are that the computational box must be periodic and Fourier transforms are intensively used.
The external forcing in Eq. (25) consists of a superposition of random Fourier modes: where F k is a Gaussian vector of zero mean and unit variance satisfying F k · k = 0 (incompressibility condition), and F * k = F −k . N f is a normalisation constant to impose that |F ext (x)| 2 = f 2 0 . A second choice of forcing is the so-called frozen forcing. It is simply obtained (after each time step) by setting the velocity field equal to a prescribed field for wave-vectors in a predefined range. This second forcing mimics a physical forcing working at constant velocity.
The pressure in the Navier-Stokes equation ensures the incompressible condition. The condition is easily satisfied by inverting the Poisson equation Let us denote by P the projector into the subspace of divergence-free functions, and define The equation of motion for the normal fluid simplifies to ∂v n ∂t In practice, since the projector P takes a trivial form in Fourier space, we directly solve Eq. (36). In addition, we employ a second-order Runge-Kutta time-advancement scheme. Note that the force acting on the normal fluid due to the interaction with the vortex lines also needs to be projected. Finally, we remark that the mean normal fluid velocity evolves in time due to the friction force as a consequence of momentum transfer between components: We refer to established literature for further details concerning the standard pseudo-spectral algorithm employed for the numerical integration of the Navier Stokes equations [96].

Interpolation
The equation of motion of the vortex lines, Eq. (28), needs the values of the normal fluid velocity at the Lagrangian discretisation points along the vortex lines, where v n is not known. In principle, as the normal fluid is periodic, any inter-mesh value can be computed exactly (within spectral accuracy) by using a Fourier transform. We call this interpolation scheme Fourier interpolation and use it as benchmark for comparisons. The Fourier interpolation is extremely costly and prohibitive for practical applications, as it requires n x n y n z evaluations of complex exponentials for each Lagrangian point along the vortex lines. Affordable interpolation schemes are typically defined on physical space, and their accuracy depends on the order of the method and on the regularity of the field. In three dimensions, the most commonly used schemes are Nearest neighbours, where the value of the closest grid point is taken, and Tri-Linear or Tri-Cubic interpolation, performed in each direction by a line or a cubic polynomial. More recently, a new approach based on fourth-order B-splines has been proved to be specially well adapted to pseudo-spectral codes and found to be non-expensive and highly accurate [97]. In order to study the effect of different interpolation schemes, we study how a superfluid vortex ring moves in a static spatially dependent flow. We do not take into account yet the retroaction of the vortex filament on the normal fluid which we keep fixed. We consider a domain of size L x = L y = L z = 2π; for the normal fluid we choose the following superposition of ABC flows at different scales: v n (x, y, z) = n max n=1 (B cos ny + C sin nz, A sin x + C cos nz, A cos nx + B sin ny), (38) with A = 1, B = −1, C = 5 and n max = 10. We evaluate the field on a mesh with N = 128 and 256 collocation points in each direction. As the largest wavevector is √ 3 n max , the Fourier interpolation is exact (up to round-up errors).
We initialise a vortex ring of size R = 0.2387 (in non-dimensional units) and set the temperature at T = 1.95 K. We let the ring evolve with the static normal fluid in the background. As the ring evolves, it deforms in the presence of the highly non-homogeneous normal fluid. In order to obtain a quantitative comparison between different schemes, we measure the average radius of the vortex ring as a function of time R(t) = |s(ξ ) − s center | where s center = s and the average is performed over vortex points. The temporal evolution of R(t) is displayed in Fig. 2 for two different resolutions and different interpolation schemes. As expected, it is apparent that the error decreases when the number of collocation points of the normal fluid grid increases. It is also clear that tri-linear interpolation gives poor results.
To better quantify the accuracy of the interpolation methods, we compute the relative error of the radius with respect to the reference radius obtained by Fourier interpolation The time evolution of this relative error is displayed in Fig. 3. Remarkably, the B-spline interpolation reduces the interpolating errors considerably. The extra cost of the B-spline scheme is just one single fast Fourier transform independently of the number of points to be interpolated, and this is why we choose it for our local selfconsistent approach, FOUCAULT.

Friction distribution
We now discuss the back-reaction force of the vortex lines on the normal fluid. This force, seen from the normal fluid, is δ-supported on the N p Lagrangian points along the vortex lines {s i (t)} i=1,...,N p and needs to be distributed on the grid points where the normal fluid is computed. More precisely, we denote by (ζ, μ, χ), with ζ, μ, χ ∈ [0, 1), the neighbouring grid x y z points of s i (t). The force f ns (s i (t)) exerted by the vortex line on the normal fluid has to be distributed among neighbouring points using weights w ζ,μ,χ , as shown schematically in Fig. 4.
By definition, the weights satisfy where f i ns (t) := f ns (s i (t)), δ i is the length of the i-th vortex line element and the dependence on time of the friction force is explicitly showed, as it will play an important role in the subsequent part of the discussion. For the sake of simplicity, the external force has been omitted in Eq. (39) and a simple Riemann sum has been used to approximate the line integral. Note that a trapezoidal rule, which is a better approximation, may be numerically used by the simple replacement δ i → (δ i + δ i+1 )/2 and readjusting the indices of the sum. It is thus evident that our problem is formally identical to coupling discrete, point-like, active particles to the turbulent flow of a classical viscous fluid.
It is well known that in problems of active matter properties such as e.g. aggregation and condensation, mixing or/and growing of species may strongly depend on the choice of the force distribution method and the filtering scale [35,36]. The same issues appear in our self-consistent local model FOUCAULT. In order to illustrate this problem, we consider the simple case of a vortex ring moving in a initially quiescent normal fluid, similar to the setting studied in [30]. We consider a ring of radius R = 0.2387 and set the temperature at T = 1.95 • K. Because of friction, we expect the ring to shrink. We compare the temporal evolution of its radius employing different filtering methods. We distribute the force to the nearest neighbour grid points, i.e. w ζ,μ,χ = 0 for all grid points except the nearest one to s i . The resulting force field is then filtered by using either a moving average over N filter points or a Gaussian kernel of width N filter Δx, where Δx = L x /N x is the mesh size. Figure 5 displays the temporal evolution of the vortex ring radius for the different schemes. It is clear that the shrinking rate of the radius of the ring depends strongly on the filtering procedure which is employed.
This dependence is clearly spurious. A numerical method based on physics principles is hence needed. Since our context is turbulence, it makes sense to adopt the same rigorous regularisation approach which has been used to take into account the strongly localised response of point-like particles in classical turbulence [35,36]. The advantage of adopting this method is that the regularisation of the exchange of momentum between point-like particles and viscous flows is based on the physics of the generation of vorticity and its viscous diffusion at very small scales. In our case, the justification of this method arises from the very small Reynolds numbers of the flow based on the vortex core radius and the small velocity of the vortex lines with respect to the normal fluid.
We refer the reader to the original papers [35,36] for further details; in brief, this approach which we borrow from classical turbulence is based on the solution of the delta-forced, linear unsteady Stokes equation. Accordingly to the classical formulation, we introduce the timedelay R coinciding with the time interval during which the localised vorticity (generated by the relative motion between the vortex line and the normal fluid) is diffused to the relevant hydrodynamic scale of the flow, the spacing Δx of the computational grid which the normal fluid velocity field is calculated on. This finite time delay, R , regularises the delta-shaped nature of the friction force f i ns (t)δ(x − s i (t)) in a natural way via the fundamental solution of the diffusion equation: This solution, Eq. (40), is a Gaussian with standard deviation σ R = √ 2ν R . The resulting expression for the friction force exerted by the i-th vortex line element on the normal fluid at the time t at the point [35,36], yielding the following modified Navier-Stokes equation for the normal fluid velocity field: From the physical point of view, Eq. (41) implies that the strongly localised vorticity injected in the normal flow by the relative motion of the vortices is neglected until it has been diffused by viscosity to a characteristic length scale σ R = √ 2ν R . To be consistent, and in order to take into account the vortex induced disturbances as soon as the relevant hydrodynamic scales are affected, we choose the finite time delay R so that σ R /Δx = 1. Extensive tests performed in the original paper [35] ensure that σ R /Δx = 1 is a suitable choice.
In theory, for each i-th vortex line element, it is possible to compute the corresponding weight for each point of the numerical grid: it would be sufficient to integrate the Gaussian kernel g[x − s i (t − R ), R ] over the volume ΔV = ΔxΔyΔz centred on the grid point. For instance, as displayed in the schematic two-dimensional Fig. 6, the force generated by the i-th vortex element at s i contributes to the force computed on mesh point x, but also on mesh point z, the weights being the integral of g(x − s i , R ) over their respective cells, represented by black dashed squares in the Fig. 6. This procedure would have to be repeated on each grid point for each vortex element. Although exact, this approach would be extremely costly. In addition, referring again to Fig. 6, the weight corresponding to mesh point z is very small as we choose σ R /Δx = 1; the diffusion-based regularisation which we employ localises in fact the force in a sphere of radius σ R centred at s i . As a consequence, instead of distributing the force on each grid point by computing N x N y N z N p integrals g(x − s i , R ) all over the grid, we distribute the force only to neighbouring mesh points of s i , taking care of including the weights of far grid points. Proceeding in this fashion, the total force is preserved, as the integral , R ]dx = 1. As matter of example, in the sketch outlined in Fig. 6, the  Fig. 6, this corresponds to integrating g over the region coloured in green for point x and over the yellow region for point y. The generalisation to three dimension is straightforward. The advantage of this approach is that the space integrals of g may actually be computed analytically.
The computation of the integrals leads to the weights w ζ,μ,χ that will be used to distribute the force among the neighbouring grid points of s i , as illustrated in Fig. 4. Note that because the Gaussian kernel g is normalised to one, the requirement that 1 ζ,μ,χ=0 w ζ,μ,χ = 1 is satisfied by construction. To compute the weights w ζ,μ,χ , we first note that the integrals of g(x−s i , R ) can be factorised in each Cartesian direction: we will therefore only compute the contribution in the x direction, the computation of the contributions of the other directions is formally identical. The weight for the grid point of the left of s x i (indicated with s x i , corresponding to ζ = 0) results from the one-dimensional integral over (−∞, s x i + Δx 2 ), whereas the weight for the grid point on the right ( s x i + Δx, ζ = 1), stems from the integral over ( s x i + Δx 2 , ∞). The calculation of the one-dimensional weight is straightforward; we have: The integrals w μ [s

Numerical strategy and parallelisation
We have implemented the numerical integration of the fully-coupled local model taking advantage of modern parallel computing. The solvers of the VFM and the Navier-Stokes equations are of very different nature, but they need to interact only through the evaluation of the friction force. In this first version of the solver, we have opted for an hybrid OpenMP-MPI parallelisation scheme. The two solvers are handled by different MPI processes. Each MPI process contains many OpenMP threads, so that each solver is also independently parallelised by using this shared memory library. The evaluation of the friction force requires communication between the two solvers that do not have access to each other fields and variables. This communication is managed through the Message Passing Interface (MPI) library at the end of each time step. This scheme is naturally adapted to modern clusters that contain many nodes, with each node having a large number of CPUs with shared memory. In order to clarify this numerical strategy, in Appendix A we report a chart illustrating schematically the steps performed by each MPI process and the communication between the latter at the end of each time-step. Finally, we remark that the values of the time step for the Navier-Stokes equations and for the VFM solver need not be the same; typically, the VFM requires a smaller time step. In order to speed up the code, depending on the physical problem, each solver can perform sub-loops (in time) to ensure numerical stability and an efficient integration of the full model.

Applications
In this section, we show two physical applications of our fully-coupled local model FOU-CAULT : a superfluid vortex ring moving (i) in an initially stationary normal fluid and (ii) in a turbulent normal fluid. Our aim is not to introduce any new physical phenomena but simply validate our approach and show its computational capabilities.
(i) Vortex ring in initially quiescent normal fluid We first study the vortex configuration investigated in the pioneer work of Kivotides et al. [30]: a vortex ring moving in a initially quiescent normal fluid. It was observed that, due to the interaction between superfluid and normal fluid, two concentric vortex rings are created in the normal fluid, accompanying the superfluid vortex ring, forming a triple vortex structure. We integrate our fully coupled model using as initial condition a superfluid vortex ring of radius R = 0.2387 in a box of size 2π and set the temperature to T = 1.95 K. As the superfluid vortex ring travels, it sets in motion the surrounding normal fluid, thus losing energy and shrinking. A visualisation of the resulting flow is displayed in Fig. 7.
The superfluid vortex ring is displayed in green and the two normal fluid vortex rings which are generated are rendered in orange-reddish colours, forming the same triple vortex structure discovered by Kivotides et al. [30]. On the plane perpendicular to the vortex rings, we also show the normal fluid translational velocity in the direction of propagation of the three vortex rings, rendered in blue-white colours: a wake (or jet) behind the superfluid ring, recently identified [33], is apparent. In this case study, the number of vortex discretisation points is equal to 200 and the normal fluid resolution is 256 3 . The complete shrinking of the vortex ring took 2 hours of computational time. (ii) Vortex ring in turbulent normal fluid.
Our second application is the effect of turbulence in the normal fluid component on the dynamics of a superfluid vortex ring. As discussed in Sect. 3.2, we can easily add an external forcing F ext to the Navier-Stokes equations in order to sustain the normal fluid turbulence. We first produce a turbulent state using a volumetric external forcing at large scales (k sup = 1). We use a resolution of 256 3 points, which allows us to obtain a Reynolds number based on the Taylor microscale equal to Re λ = 127. Once the turbulent state is prepared, we study the evolution at temperature T = 1.95 K of a superfluid vortex ring of radius R = 35η, where η is the Kolmogorov length scale of the flow. The initial condition used in our self-consistent local model is displayed in Fig. 8 (top left panel). We clearly observe the turbulent fluctuations of the normal fluid displayed by the reddish rendering of the normal fluid enstrophy. As time evolves, the turbulent fluctuations destabilise the superfluid vortex ring, inducing large amplitude Kelvin waves (bottom left panel); After a couple of large-eddy-turnover times, the vortex ring self-reconnects, forming vortex loops which undergo further reconnections, creating a turbulent superfluid vortex tangle in equilibrium with the turbulent normal fluid. It is apparent from Fig. 8 that normal fluid fluctuations are responsible for an increase of the total vortex length. In Fig. 9, we display the temporal dependence of the total superfluid vortex length: the large amount of stretching is evident. The right panel shows that initially the vortex length does not increase exponentially, as predicted for the Donnelly-Glaberson instability [98] for a uniform normal flow, but with approximate power-law character. The figure also shows that, remarkably, there is a change of behaviour at times of the order of one T L , suggesting that normal fluid fluctuations may play an important role in quantum turbulence. This is an interesting question with important implications for current quantum turbulence experiments. It is natural to expect that the evolution of the vortex length should depend on temperature, but be independent [51] of the initial vortex configuration. Clearly, an accurate and physically sound solver of the fully coupled dynamics of normal fluid and superfluid like the solver that we have presented here is an important tool for the next studies of quantum turbulence. Concerning the details of the computation performed, the total computational time was 72 hours with a final number of vortex points equal to 131290. It is worth noting that in this circumstance of a very dense tangle, the heaviest calculation regards the computation of the evolution of the vortex tangle. Despite in fact the employment of the tree approximation [84] included in the code, the computational node advancing in time the vortex positions is absorbing most of the computational resources compared to the Navier-Stokes solver and the handling of the coupling.

Conclusions
We have presented a novel algorithm, named FOUCAULT, for the numerical simulations of quantum turbulence in helium II at nonzero temperatures. The main features of our approach are (i) the parallelised, efficient, pseudo-spectral code, capable of distributing the calculation amongst distinct computational cluster codes, and (ii) the forcing scheme employed for the normal fluid. These features allow the resolution a wider range of length scales compared to previous studies, from large quasi-classical scales to small quantum length scales smaller than the average inter-vortex distance. This is of fundamental importance in order to investigate the character of quantum turbulence shared with its classical counterpart, as well as the features which are specific to turbulent superfluid flows. From the point of view of the physical modelling, a novel feature of our approach is the local computation of the friction force, stemming from classical creeping flow analysis, which (unlike previous work) is independent of the numerical discretisation on the vortex lines. In addition, our approach implements an exact regularisation of the friction force based on the small-scale viscous diffusion of the disturbances generated by vortex lines moving with respect to the normal fluid.
After a detailed description of the distinct features of our algorithm, we have applied it to two problems: the motion of a superfluid vortex ring in an initially quiescent normal fluid and in a turbulent normal fluid. In the first problem, we have recovered the same qualitative features (the double normal ring structure and wake) observed in previous work; in the second we have demonstrated a new effect-the vortex stretching induced by turbulent fluctuations.
In conclusion, our novel algorithm paves the way for a new series of investigations of quantum turbulence in helium II at finite temperatures, particularly at length scales smaller than the average inter-vortex spacing in problems where the singular and quantised nature of the superfluid vortices is expected to play a fundamental role.

Data Availability Statement
This manuscript has associated data in a data repository. [Author's comment: This manuscript has associated data in a data repository. Additional metadata are available at: https://doi.org/ 10.25405/data.ncl.12587684.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A. Schematic illustration of numerical strategy and parallelisation
In this appendix, we include a chart (Fig. 10) describing schematically the parallelisation strategy. As reported in the main text, we have two distinct solvers (the VFM and the Navier-Stokes (NS)) each corresponding to an MPI process. Each MPI process has access to OMP_NUM_THREADS OpenMP threads that share memory. The two MPI processes only need to communicate when exchanging the normal fluid velocity field v n (x, t) (NS → VFM) and the mutual friction force F ns (x, t) (VFM → NS).