This paper investigates the applicability of smoothed particle hydrodynamics in simulations of superfluid helium4. We devise a new approach based on Hamiltonian mechanics suitable for simulating thermally driven and weakly compressible flows with free surfaces. The method is then tested in three cases, including a simulation of the fountain effect. We obtain remarkable agreement with referential and theoretical results. The simulations provide new physical insight, such as the pressure and temperature fields in a vessel experiencing the fountain effect.
I. INTRODUCTION
Quantum fluids are governed by laws of quantum mechanics rather than classical ones and represent a fairly recent addition to the family of various types of fluids with relevance both in fundamental science and in technological applications. Prime examples can be found in superfluid phases of helium isotopes or in Bose–Einstein condensates of rarefied gases. These fluids possess remarkable properties, such as the possibility of frictionless flow, the existence of numerous sound modes, the formation of superfluid films, quantized circulation around discrete vortex lines, or the ability to transport heat convectively with zero local mass flow. In addition to the fields of lowtemperature physics and fluid dynamics, this behavior makes quantum fluids appealing, via certain analogies,^{62} for researchers across multiple disciplines, including astrophysicists^{13,44} or cosmologists,^{22,65} for which they represent an accessible model system.
In this work, we deal specifically with the isotope helium4. The liquid phase, historically called He I, exists below 4.2 K at standard pressure. When cooled further by pumping on its vapors, helium4 undergoes a secondorder phase transition at 2.17 K and enters the superfluid phase, called He II. Most physical properties of helium are wellknown and tabulated.^{6}
He II exhibits all of the exotic phenomena mentioned above; for instance, it flows easily through narrow capillaries without friction and displays twofluid behavior (in the sense of Landau's model, see further below). The latter leads to the mechanocaloric effect, or the famous fountain effect,^{1} where a heater causes the motion of helium in the form of a fountainlike jet.^{2} Additionally, various forms of turbulent motion may exist in He II, including a state of twofluid “double turbulence;” for a phenomenological treatment of quantum turbulence, see Ref. 47.
Historically, experimental investigations of dynamics in He II may be divided into studies of thermally generated flows, such as the fountain effect, counterflow in a channel driven by a heater placed at its closed end,^{52} or mechanically driven flows, mostly by submerged oscillating structures.^{14,20,43} A superfluid wind tunnel experiment has been constructed,^{41} and a largescale Kármán flow experiment is now in operation.^{40} Methods for flow visualization using tracer particles^{25,26} or He excimers^{42} have been developed recently and applied to both thermally and mechanically driven flows.
The motion of superfluid helium4 can be simulated by solving the Gross–Pitaevskii equation,^{46} the Hall–Vinen–Bekarevich–Khalatnikov (HVBK) models,^{5,15} the onecomponent models,^{35} or by the vortexfilament method.^{10,19} However, many experiments are notoriously difficult to simulate numerically using meshbased techniques such as finite element or finite volume methods, as they contain either a free surface (fountain effect) or a moving boundary (oscillating bodies). Therefore, it is advantageous to use a method based on Smoothed Particle Hydrodynamics (SPH),^{31,33,58} which is meshfree and thus more suitable for these problems.
This paper presents a new numerical meshfree method for simulations of superfluid helium4. An SPH method for helium II has already been developed by Tsuzuki,^{54–56} who investigated vortices emerging in a rotating cylinder. Our approach is different because we discretize helium II using only one type of fluid particle with two densities as opposed to the work of Tsuzuki, where normal and super particles are distinguished. Our method comes at the cost of including additional convective terms, which consider that entropy, for instance, is advected by normal velocity and not by the collective or coflow velocity. However, it liberates us from assuming that superfluid behaves like a mixture and facilitates the implementation of heating bodies (a crucial feature in thermally driven flows, like the fountain effect). Additionally, we use a slightly different form of equations, a Hamiltonian variant of the twofluid model that goes beyond the Landau–Tisza approach. We explain this in detail in Appendix A. The Hamiltonian characteristic of the equations is advantageous also in the SPH numerical scheme.^{23}
The first macroscopic models of superfluid helium4 were proposed by Tisza et al.^{16,27,49,51} Their final model^{21} consists of four evolution equations, namely, the continuity equation, balance of momentum, evolution of superfluid velocity, and entropy balance. The Landau–Tisza model can also be derived from the Gross–Pitaevskii equation.^{53} However, the model has several limitations. First, it does not allow for nonzero superfluid vorticity (quantum vortices).
Second, the Landau–Tisza model is formulated in terms of five quantities (superfluid density ρ_{s}, normal density ρ_{n}, superfluid velocity $ v s$, normal velocity $ v n$, and entropy density s), despite having only four evolution equations. This inconsistency is overcome by setting a dependence of the ratio $ \rho n / \rho $ on temperature ( $ \rho = \rho s + \rho n$ being the total mass density). Although this setting closes the evolution equations, it goes against the nature of superfluid helium4, which is not a mixture of two fluids, but rather a single fluid with two motions, as expressed by Landau:^{27,29} “It must be particularly stressed that we have here no real division of the particles of the liquid into ‘superfluid’ and ‘normal’ ones….” The dependence $ \rho n / \rho ( T )$ actually defines how free energy depends on temperature, which has to be taken into account because otherwise compatibility with Hamiltonian mechanics would be violated.^{50} The HVBK model resolves the problem of four equations for five quantities by requiring both $ v n$ and $ v s$ to be divergencefree (incompressible).^{4,38} Another solution was proposed by Zilsel,^{64} who introduced a continuity equation for both ρ_{s} and ρ_{n}, despite that there are no real two components in the superfluid helium (as noted by Landau). Yet another solution is the socalled onefluid model,^{35} where entropy exhibits an extra motion instead of density.
Here, we follow a route that builds upon Hamiltonian structures of the evolution equations for superfluid helium4. The Hamiltonian structure of Landau–Tisza equations was derived from quantum commutators between the mass density and phase of the wave function.^{7,27} Moreover, Volovik and Dotsenko extended the Hamiltonian structure to a model related to HVBK dynamics by also considering Hamiltonian motion of quantum vortices.^{59–61} Alternatively, the Hamiltonian form of equations can be derived from the onefluid model.^{50} This is the structure we use in the current manuscript. Another Hamiltonian formulation of HVBK dynamics was obtained by Holm and Kuperschmidt,^{11,17,18} which contains additional state variables (vector and scalar potentials of the superfluid velocity).
II. SPH TWOFLUID MODEL
This section contains a twofluid model of superfluid helium4 that extends the Landau–Tisza model and its discretization within SPH.
A. Reversible part of the model
B. Irreversible forces
Although SPH can compute reversible flows, the stabilization of entropy helps to prevent the selfemergence of disorder. Numerical noise is a recognized problem in explicit SPH, which makes velocity converge to Boltzmann distribution and eventually leads to the unphysical gaslike behavior of simulated particles.^{23} We stabilize by adding a “small Laplacian” to the righthand side of the entropic balance. Particle shifting^{32} would be an alternative remedy.
C. Constitutive laws
.  Cavity .  Second sound .  Fountain effect . 

T_{0}  $ 1.9 \u2009 K$  $ 1.65 \u2009 K$  
ρ_{0}  1  $ 145.5 \u2009 kg / m 3$  $ 145.2 \u2009 kg / m 3$ 
s_{0}  $ 725.5 \u2009 J / ( kg \u2009 K )$  $ 335.0 \u2009 J / ( kg \u2009 K )$  
$ \chi n 0$  1  0.4195  0.1934 
$ \chi \u2032$  0  $ 5.697 \xd7 10 \u2212 4 \u2009 kg \u2009 K / J$  $ 5.851 \xd7 10 \u2212 4 \u2009 kg \u2009 K / J$ 
u_{1}  20  $ 40 \u2009 m / s$  $ 40 \u2009 m / s$ 
u_{2}  $ 18.83 \u2009 m / s$  $ 20.37 \u2009 m / s$  
C  $ 3902 \u2009 J / ( kg \u2009 K )$  $ 1861 \u2009 J / ( kg \u2009 K )$  
μ  $ 1 Re$  0  $ 10 \u2212 4 \u2009 kg / ( m \u2009 s )$ 
β  0  0  $ d \rho 0 s 0 2 200 u 2$ 
δr  $ 1 N$  $ L N$  $ d 40$ 
h  $ 3 \delta r$  $ 3 \delta r$  $ 2.8 \delta r$ 
δt  $ h 10 u 1$  $ h M u 1$  $ h 5 u 1$ 
$ t end$  40  $ 2 L u 1$  $ 0.3 \u2009 s$ 
H  $ 2 \xd7 10 \u2212 3 \u2009 m$ 
.  Cavity .  Second sound .  Fountain effect . 

T_{0}  $ 1.9 \u2009 K$  $ 1.65 \u2009 K$  
ρ_{0}  1  $ 145.5 \u2009 kg / m 3$  $ 145.2 \u2009 kg / m 3$ 
s_{0}  $ 725.5 \u2009 J / ( kg \u2009 K )$  $ 335.0 \u2009 J / ( kg \u2009 K )$  
$ \chi n 0$  1  0.4195  0.1934 
$ \chi \u2032$  0  $ 5.697 \xd7 10 \u2212 4 \u2009 kg \u2009 K / J$  $ 5.851 \xd7 10 \u2212 4 \u2009 kg \u2009 K / J$ 
u_{1}  20  $ 40 \u2009 m / s$  $ 40 \u2009 m / s$ 
u_{2}  $ 18.83 \u2009 m / s$  $ 20.37 \u2009 m / s$  
C  $ 3902 \u2009 J / ( kg \u2009 K )$  $ 1861 \u2009 J / ( kg \u2009 K )$  
μ  $ 1 Re$  0  $ 10 \u2212 4 \u2009 kg / ( m \u2009 s )$ 
β  0  0  $ d \rho 0 s 0 2 200 u 2$ 
δr  $ 1 N$  $ L N$  $ d 40$ 
h  $ 3 \delta r$  $ 3 \delta r$  $ 2.8 \delta r$ 
δt  $ h 10 u 1$  $ h M u 1$  $ h 5 u 1$ 
$ t end$  40  $ 2 L u 1$  $ 0.3 \u2009 s$ 
H  $ 2 \xd7 10 \u2212 3 \u2009 m$ 
D. Discrete space
E. Discrete time
Simple explicit integrators, like leapfrog, are commonly used in SPH codes. In this paper, we employ the following scheme, which is very similar to leapfrog, except we update entropy and density using midtime positions:

(initial step only) find rate of $ v$ and $ v s$

update $ v$ and $ v s$ by $ \delta t 2$ step

update x by $ \delta t 2$ step

recompute the cell list and find ρ

find the rate of s

update s by $ \delta t 2$ step

update x by $ \delta t 2$ step

recompute the cell list and find ρ

find the rate of $ v$ and $ v s$

update v and $ v s$ by $ \delta t 2$ step
In each step, rates are evaluated using the most recently computed values of $ s , \rho , v , v s , \u2009 and \u2009 x$.
III. NUMERICAL EXPERIMENTS
A. Liddriven cavity
Before we jump to examples involving actual superfluidity, we assess the validity of our model for classical fluid (Navier–Stokes equations), which can be understood as a special case of system (1) for $ \chi s = 0$. Clearly, our discrete equations (19) must work for classical fluids; otherwise, there would be no hope to use them on more complex problems.
Liddriven cavity is usually formulated as a dimensionless problem: Incompressible, viscous flow is confined in a box $ ( 0 , 1 ) \xd7 ( 0 , 1 )$. No slip boundary is prescribed at the left, bottom, and right walls. At the top, velocity is $ v = ( 1 , 0 ) T$, density is $ \rho 0 = 1$, and viscosity $ \mu = 1 Re$. Reynolds number $ Re$ has various values.
Walls are implemented using a layer of dummy particles, which are treated normally with the exception that their velocities are constantly zero. The lid is also implemented this way, but with velocity $ v a = ( 1 , 0 ) T$ appearing in the evaluation of viscous forces. The list of simulation parameters is displayed in Table I. We measure transverse velocities along x and y centerlines and compare them to a reference solution by Ghia et al.^{12} The results are shown in Figs. 1–3 and convergence curve in Fig. 4. Error grows approximately linearly with $ \delta r = 1 N$.
B. Second sound waves
C. Fountain effect (twodimensional)
Unlike the two previous examples, in superfluid fountain effect, coflow and counterflow are equally significant. Let us begin by describing our idealized setting. A cell (with shape indicated in Fig. 7) is submerged in liquid helium4 at temperature $ T 0 = 1.65 \u2009 K$. In its center, there is a point heat source of constant power $ W \u0307$. The cell is insulated by adiabatic walls and a thin superleak at the bottom.
Adiabatic walls of the cell are implemented using static dummy particles, which take part only in density and pressure computation. At the walls of the cryostat, let us prescribe noslip for $ v$, free condition for $ v s$, and Dirichlet condition S = S_{0} for entropy (which is almost the same as T = T_{0}). Consequently, the wall acts as a cooler, preventing the temperature from growing indefinitely. This condition is implemented using dummy particles, which we force to have constant entropy and zero coflow velocity, but we allow them to conduct heat.
The simulation never reaches a stationary state, but the jet speed is relatively stable after a brief transitional period. Figure 8 shows that the computed jet speed is very close to the value (31) but slightly lower on average. Figure 9 shows the shape of the fountain, velocity, and pressure. The simulation remains stable even during the violent impact of jetted helium. Figure 10 shows the temperature in the first millisecond of simulated time. After this short period, second sound waves fill the cell, and temperature gradients vanish, except for a steep jump at the superleak. The subsequent evolution is shown in Figs. 11 and 12. When the fountain connects with colder helium in the cryostat, a temporary drop in temperature occurs. Nonetheless, the effect of the second sound traveling upstream and entering the cell through the capillary was not observed. Table I lists all simulation parameters for reproduction purposes.
IV. CONCLUSION
We described a new energyconserving SPHbased numerical method for superfluid helium. We verified its validity in three essential test cases, demonstrating that the technique can accurately capture coflow and thermally driven counterflow. Our fountain effect simulation verifies the applicability of the formula for fountain jet speed by Landau and Lishitz.
The numerical scheme discretizes a Hamiltonian onefluid formulation generalizing the Landau–Tisza twofluid model. We also provide a closed form of the energy functional in terms of density, coflow velocity, entropy density, and superfluid velocity, which has been missing in the literature.
In the future, we would like to enrich our discrete system by implementing the dissipative effects of quantum vorticity. Particle simulation of Rollin films is another promising research direction, as well as simulation of helium droplets,^{24} or simulation of inertial particles in superfluid helium.^{45} Moreover, the meshfree nature of SPH is advantageous in multiphase flow simulations, including phase change.^{34,63} Consequently, we expect that this method could be extended to simulate solidification and evaporation in superfluid helium.
ACKNOWLEDGMENTS
O.K. and D.S. were supported by Project No. START/SCI/053 of Charles University Research program. M.P. was supported by Czech Science Foundation, Project No. 23‐05736S.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Ondřej Kincl: Investigation (equal); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (lead). David Schmoranzer: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Michal Pavelka: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available within the article.
APPENDIX A: DERIVATION OF THE MODEL
This appendix contains a derivation of model (1) from a Hamiltonian formulation of superfluidhelium dynamics.^{50,59} The field of superfluid density interacts with itself partly due to the convective term and partly due to the presence of quantum vortices. When the reversible (nondissipative) part of interactions mediated by quantum vortices is taken into account, we obtain a system of four Hamiltonian equations for four state variables (mass density ρ, momentum density $m$, volumetric entropy density $ s \xaf$, and superfluid velocity $ v s$) as follows:
In summary, the reversible part of our model [Eqs. (1)] together with the identification of derivatives of energy (9) and a form of the energy itself (determined up to the internal energy) are implied by Hamiltonian mechanics of superfluid helium4 [Eqs. (A1)] and by the requirement of Galilean invariance. The model is more precise than the standard Landau–Tisza model because it takes into account at least the reversible interaction between the field of superfluid velocity with itself due to the quantum vortices. This results, for instance, in the presence of convective derivative $ \u2202 t + v \xb7 \u2207$ in evolution equations for all the state variables, instead of having the superfluid velocity $ v s$ convected only by itself.^{53}
APPENDIX B: METHOD DERIVATION AND CONSERVATION LAWS
1. Reversible part
Let us imagine a classical solution to these equations, which lives in a timedependent domain Ω_{t} and such that the boundary $ \delta \Omega t$ is an adiabatic free surface, where we assume $ v n s \xb7 n = 0$ and p = 0. (Pressure is only relevant up to an additive constant. We fix this constant so that vapor pressure is zero.) The adiabatic free surface is our choice of donothing boundary condition. (Implementation of other boundary conditions, such as at the walls of a cryostat, necessitates specific treatment.)
2. Irreversible part
We will need two additional operators, both of which are well known and established.
In analogy to (B18), we define