Particle approximation of the two-fluid model for superfluid 4He using smoothed particle hydrodynamics

This paper presents a finite particle approximation of the two-fluid model for liquid 4He using smoothed particle hydrodynamics (SPH). In recent years, several studies have combined the vortex filament model (VFM), which describes quantized vortices in superfluid components, with the Navier–Stokes equations, which describe the motion of normal fluids. These studies led us to assume that coupling both components of the two-fluid model instead of using the VFM to describe the superfluid component enables us to approximate the system. In this study, we formulated a new SPH model that simultaneously solves both equations of motion of the two-fluid model. We then performed a numerical simulation of the rotating liquid 4He using our SPH. The results showed that the two major phenomena, the emergence of multiple independent vortices parallel to the circular axis and that of the so-called rigid-body rotation, can be reproduced by solving the two-fluid model using SPH. This finding is interesting because it was previously assumed that only a single vortex emerges when addressing similar problems without considering quantum mechanics. Our further analysis found that the emergence of multiple independent vortices can be realized by reformulating the viscosity term of the two-fluid model to conserve the angular momentum of the particles around their axes. Consequently, our model succeeded in reproducing the phenomena observed in quantum cases, even though we solve the phenomenological governing equations of liquid 4He.


Introduction
The bizarre behavior of liquid helium 4 (hereinafter, 4 He) has attracted the attention of condensed matter physicists for many years; however, the detailed dynamics of this fluid remain a mystery even now, eighty years after its discovery. Figure 1 presents an overview of the basic properties of 4 He in the form of a schematic phase diagram. In 1937, P. L. Kapitza found that the viscosity of liquid 4 He decreased upon cooling and almost reached zero below the λ-point (labeled point c in figure 1). Liquid 4 He was later confirmed to display viscosity in the hightemperature area above the λ-line, and to lose it in the low-temperature area below the λ-line. Once lack of viscosity is achieved in this latter temperature range, liquid 4 He displays many unexpected properties that characterize it as a 'superfluid.' By contrast, liquid 4 He within the high-temperature area above the λ-line is still classified as a 'normal fluid. ' The superfluidity condition is characterized by properties such as those in the following examples. Superfluid 4 He can penetrate a high-density capillary filter (e.g., gypsum), owing to the so-called superleak effect, whereas this is not possible for normal fluid 4 He. Liquid 4 He in a vessel equipped with a capillary filter at its bottom and an opening at its top, and dipped in a reservoir, gushes from the top of the vessel because the surrounding superfluid 4 He flows back into the container owing to the superleak effect. This process is known as the 'fountain effect.' In addition, the 'film flow effect' takes place when superfluid 4 He spontaneously comes out of its vessel as a result of the van der Waals forces between the helium atoms and walls being stronger than the forces between the 4 He atoms themselves.
Many physicists have made efforts to explain these interesting properties of superfluids at both the ontological and phenomenological levels. One year after the discovery by P. L. Kapitza, F.W. London explained the phase transition of liquid 4 He from a quantum mechanical standpoint, comparing it with a type of Bose-Einstein condensation (BEC) [1]. London's theory has been validated by many experiments. Apart from this, the 'two-fluid model' proposed by L. Tisza [2] and L. Landau [3], which assumes the coexistence of superfluidity and normal fluidity in liquid 4 He, has become the most acknowledged phenomenological description. A remarkable achievement of the two-fluid model is that it connects the macroscopic perspective of thermodynamics with quantum mechanics, such that the superfluid can be described within the framework of basic mechanics (for detailed descriptions of the two-fluid model, see section 2.1).
From the viewpoint of numerical analysis, a finite approximation of the phenomenological governing equations represents an essential step toward the direct numerical simulations of large-scale problems on a continuum mechanical scale, for example, the fountain effect or film-flow effect. However, to date, only a few studies have adopted such an approach; this is primarily because most previous studies took a quantum mechanical standpoint (e.g., numerical analyses of nonlinear Schrödinger equations [4,5], and quantum Monte Carlo calculations [6,7] with the aim of examining the static responses [8] or to explore quantum vortices and turbulence [9][10][11][12][13]). Conversely, previous studies [14][15][16][17] introduced direct numerical simulations of the twofluid model. Nevertheless, they focused on mesh-or grid-based approaches; some studies [14,15] involved the development of finite difference approximations, and others [16,17] led to finite element approximations. Notably, none of the previous studies reported a finite particle approximation of the two-fluid model.
Numerical simulations using the finite particle approximation are advantageous in many respects because of the discretization based on many-particle interacting systems. Each particle moves in space while interacting with every other particle, taking on either the state of a superfluid or that of a normal fluid. The mechanical picture, in this case, is much closer to the view of quantum statistical systems that follow Bose-Einstein condensates, compared with the case of grid-based methods. Using the same number of particles as real atoms in liquid 4 He in simulations remains challenging even today. Nevertheless, finite particle approximation has a high possibility of bridging the gulf between the phenomenological and microscopic concepts of liquid 4 He. For example, it should be possible to devise numerical models that replicate quantum effects such as vortex reconnections and incorporate them into the discretized two-fluid model. Furthermore, from a practical point of view, the particle approximation is capable of handling complex boundaries. It has great potential to perform simulations that are sufficiently realistic to reproduce a large-scale experimental setup. However, in addition to these potential advantages, the computational costs are much higher than those of other methods. Furthermore, because the non-uniformity of particle distributions often causes the numerical accuracy to deteriorate, it is necessary to select appropriate supportive techniques to stabilize the simulations.
The main purpose of this study is to present a scheme for the finite particle approximation of the two-fluid model. We propose to discretize the two-fluid model of superfluid 4 He using smoothed particle hydrodynamics (SPH) [18], a well-established Lagrangian particle approximation that is particularly popular in the field of astrophysics. By combining theoretical methodologies typical of two different branches of physics (lowtemperature physics and astrophysics), we aim to pioneer a new paradigm of Lagrangian particle mechanics targeting superfluids. According to quantum mechanics, in principle, it is not possible to separate the superfluid and normal fluid components in space. Let us discuss this point further. In recent years, several studies have combined the vortex filament model (VFM), which describes quantized vortices in superfluid components, with the Navier-Stokes equations, which describe the motion of normal fluids [19][20][21]. On the premise that this combined approach is acceptable, we can assume that coupling both components of the two-fluid model, instead of using the VFM for the superfluid component, would approximate the system. This approximation is allowed only when we use the two-fluid model to directly simulate large-scale problems on the continuum mechanical scale because technically, this breaks the microscopic laws of quantum mechanics. Scenarios such as these, which simultaneously solve for the two components of the two-fluid model, correspond to multi-phase flows in a classical fluid. Therefore, we allow the two components to be combined in space in the numerical tests in this study. We investigate the possibility of reproducing several quantum mechanical phenomena observed in liquid 4 He in this classical mechanical approximation.
The remainder of this paper is structured as follows. Section 2 describes the concepts of the two-fluid model and SPH, provides a summary of the wave equations in superfluid 4 He, and discusses the expression of thermodynamic parameters in the two-fluid model, in preparation for section 3. In section 3, we derive a finite particle approximation of liquid 4 He. In section 4, we discuss strategies for practical applications and introduce several improved numerical schemes for SPH. Specifically, we present a reformulation of the viscosity term in the two-fluid model to conserve the angular momentum of the fluid particles around their axes. In section 5, we demonstrate our model with three major numerical tests: the Rayleigh-Taylor instability, wave propagation analyses, and rotating cylinder simulations. Section 6 summarizes our results and concludes the paper.

Two-fluid model
The total mass density ρ is expressed as follows: where ρ n and ρ s are the mass densities of the normal fluid and of the superfluid, which respectively satisfy the laws of mass and entropy conservations as follows: Here, v n and v s are the velocities of the normal fluid and of the superfluid, respectively, while σ is the entropy density. By solving the Gross-Pitaevskii equation (a nonlinear Schrödinger equation for boson particles) and the Gibbs-Duhem equation simultaneously, we obtain the following equations of motion for liquid 4 He from a phenomenological perspective [22]: where D{ · }/Dt represents a material derivative, and v n and v s are the velocities of the normal fluid and superfluid, respectively. η n is the viscosity of the normal fluid. T is temperature, P is pressure, and σ is entropy density.
Hydrodynamic aspects of liquid 4 He have been intensively studied through experiments involving counterflow in a cylinder attached to a reservoir of the liquid at one point and a heater at another point. Under this condition, the liquid 4 He is in a separate phase of super and normal fluids, which flow in the opposite direction [23]. C. J. Gorter and J. H. Mellink found the mutual friction forces to be the reason for the discrepancy in equation (4) and equation (5) with respect to the counterflow experiments at high heat fluxes [24], which are given as Further studies by W.F. Vinen [25] through experiments with the second sound wave attenuations revealed a detailed quantitative expression for F sn as where ρ s is the mass density of a superfluid, κ is the quantum of circulation, v sn is the relative velocity of v s − v n , α is a friction coefficient, and L is the vortex line density which is time dependent. It is known that the decay of L can be described as where β v is a coefficient in the Vinen equation [26].

A brief overview of standard SPH
The fundamental concept of SPH is the approximation of the Dirac delta function δ in integral form by a distribution function W, which is called the smoothed kernel function, as follows [18]: where f is a physical value at the position r in a domain Ω and the parameter h is called a kernel radius. The smoothed kernel function W must have at least the following four characteristics; (a) W converges to the delta function δ as h approaches 0 ( d =  W lim h 0 ); (b) W satisfies the normalization condition (∫Wdr = 1); (c) W satisfies the symmetricity W( − r) = W(r); and (d) W is such that the value of W becomes 0 at a distance kh, which is called a compact support condition. Figure 2 represents a schematic view of a kernel function. An example of a smoothed kernel function W is the Gaussian, which is simply expressed as 2 2 where d is the dimension (d = 1, 2, or 3) and C is a normalization factor. The Gaussian kernel is frequently used in the field of astrophysics. In the field of incompressible SPH, in contrast, polynomial functions are often used [27,28].
Equation (10) is further discretized on the basis of the summation approximation using the position r j , the infinitesimal volume ΔV j , density ρ j , and mass m j of the jth particle as follows. where m j = ρ j ΔV j and N p is the number of fluid particles used to approximate the system. The gradient and the Laplacian of the physical value f(r i ) are axiomatically obtained from equation (12) using a vector analysis, as follows: Here, we have used the mathematical relationship ∇f = ρ[ ∇ ( f/ρ) + ( f/ρ 2 )∇ρ] to derive equation (13), so that the gradient operator has symmetricity of pressures between the ith and jth particles. For more details of the operators in standard SPH, refer to [29].

Wave equations of superfluids
According to the literature [3], the governing equations of equation (2) to equation (5) yield wave propagation equations regarding density ρ and entropy S, which produce multiple types of sound waves. Under a first approximation that neglects the viscosity term, the wave equations give two different sound waves, expressed as follows [30,31]: where c 1 and c 2 are called the first sound velocity and the second sound velocity, respectively. C v represents the specific heat at constant volume. Equation (16) can be rewritten using the number of 4 He atoms N and the mass of a 4 He atom M as where we have introduced a basic relation between the entropy density σ and the entropy S of σ = S/(NM).

Expression of thermodynamic variables in an elementary excitation model
L. Landau derived that the elementary excitation of liquid 4 He consists of two components, 'phonons' and 'rotons' [3]. The pressure P of liquid 4 He can be expressed as follows [32]: Here, let us denote the ratio of a circle's circumference as π, and the speed of sound as c. We also introduce constant values for μ, p 0 , and Δ as 1.72 × 10 −24 g, 2.1 × 10 −19 gcms −1 , and 8.9 K, respectively, according to the literature [33,34]. When T = 93 K, equation (18) can be approximated as [32,33]: where k B represents the Boltzmann constant and ÿ is the reduced Planck's constant. Note that equation (19) is similar to equation (5) in the literature [32]. The entropy S is obtained from the first-order partial derivatives of P with respect to temperature T as where we have ignored higher-order terms in the derivation of the second term in equation (21). Accordingly, the entropy density σ is obtained using the relation σ = S/(NM) as Here, both ζ and ξ are constant parameters because their components of N, M, π, c, μ, p 0 , Δ, k B , and ÿ are all constant and can be given as known parameters. Now, we have prepared all the parameters needed for the derivation of a particle approximation of liquid 4 He.

Discretization of the two-fluid model using standard SPH
First, equation (6), the motion equation for superfluid components, can be rewritten as: Using equation (13), we obtain the particle discretization of the right-hand side of equation (25), with respect to the ith particle, as follows: Here, the superscript of q i x represents the first letter of 'superfluid' or 'normal fluid'(x = s, n). The detailed evaluation of σ i is discussed in section 4.3.
The two-way coupling technique for SPH [35,36] that uses the weighted average method to compute relative velocity between two phases can be applied to the discretization of F sn . From equation (8), a mutual friction force ( ) F s i of the ith superfluid particle is calculated as follows: where C L is 2/3ρ s ακL and Ω i is the set of all the particles in the neighborhood of the ith normal fluid particles. The third term on the right-hand side in equation (25) is obtained by multiplying equation (29) by r -s 1 . Second, equation (7), the motion equation for normal fluid components, can be rewritten as: The first term on the right-hand side in equation (30) can be discretized similar to the case of equation (26). Meanwhile, the discretized expression of the second term with respect to the ith particle is obtained using equation (13) as follows.
where the breakdown of q i n is given as Meanwhile, a substitution of equation (14) into the third term directly leads to: where ν n is the kinetic viscosity, which is equal to η n /ρ n . Subsequently, the mutual friction force ( ) F n j of the jth normal fluid particle is obtained by summing up the contributions from each reaction force of the ith superfluid particle as follows: where Ω j is the set of all particles in the neighborhood of the jth normal fluid particle, while Ω i is that of the ith superfluid particle. Note that the summation inside the brackets in equation (35) is with respect not to j but to i in order for the two phases to conserve an equal amount of momentum exchange. In the end, the fourth term on the right-hand side in equation (30) is obtained by multiplying equation (35) by rn 1 . Now, the right-hand sides of equation (25) and equation (30) are appropriately discretized. The material derivative of velocities in the left-hand sides of equation (25) and equation (30) can be discretized using an explicit time-integrating scheme, e.g., the Verlet algorithm [37], in accordance with the usual manner in SPH simulations. Consequently, we have successfully derived a particle approximation of the governing equations of liquid 4 He.
4. Proposed strategies for practical application 4.1. Introduction of improved techniques to stabilize the simulations It is known that the standard SPH only permits a slight density difference between two phases when being applied to multi-phase flow problems because it does not assume a larger density gradient than that of the gradient of the smoothing kernel function in its formulation [38]. As a remedial measure, we introduce an improved SPH scheme that ensures the continuity of the pressure and pressure gradient at the interface between two phases, while satisfying mass and momentum conservation [39]. In this method, the density of the ith particle is computed as follows: Instead of equation (26), the pressure gradient term is then computed as whereP is the reduced pressure, given as˜≔ Additionally, several sophisticated techniques are appropriately introduced. First, we use a pair-wise particle collision technique [40] that gives a repulsive force when two particles are too close to each other. Specifically, the momentum conservation in the normal direction between two colliding particles is computed using this technique. Next, we introduce a one-dimensional Riemann solver [41] in each pair-wise particle, instead of using artificial viscosities. These two supportive techniques become valid only when the two particles are closer to each other than a distance equal to the initial particle distance d, as exceptional cases.
The non-uniformity of particle distributions leads to particle shortages in the local spaces in the simulation domain, which can cause the numerical accuracy to drastically deteriorate and cause unphysical pressure oscillations. Hence, the utilization of stabilization techniques to this end is indispensable. A new formulation of the pressure Poisson equation with a relaxation coefficient that contributes to obtaining smoothed pressure distributions in the implicit or semi-implicit SPH was proposed [42]. With respect to the use of explicit timeintegrating schemes in the SPH method, a weakly compressible SPH improved by adding artificial diffusive terms to the continuity equations was presented [43]. Several filtering techniques have also been reported. For example, a smoothing velocity involving interpolation among neighboring particles was proposed [29]. A density filtering technique known as the Shepard filter for the zero-th order correction among neighboring particles was developed [44]. This study adopts a pressure interpolant with a similar formulation as the Shepard filter, according to the literature [45].
In the case of classical fluids, we usually consider introducing surface tension models [46][47][48] to describe the mechanics between two phases accurately. However, it is safer to avoid using existing surface tension models in our case because these previous models build on the Young-Laplace relationship in continuum mechanics and assume the dynamics of classical turbulence. Above all, it is known that liquid 4 He exhibits quantum turbulence [49,50]. Although several recent studies report that quantum and classical turbulent systems both have common characteristics [51,52], for example, the Kolmogorov1941 energy spectrum [53], details of the correspondence between the two different turbulent systems are still in discussion. Therefore, we avoid using conventional surface tension models.
Meanwhile, we set the time increment, δt, to be smaller than the three conditions of numerical stability: the Courant-Friedrichs-Lewy (CFL) condition, the diffusion condition, and the body force condition [54]. We use an explicit time-integrating scheme, the velocity-Verlet algorithm [37,55], in accordance with the usual process for SPH simulations. It should be noted that there exists a possibility to violate stringent energy conservation when using an explicit time-integrating scheme. Regarding entropy conservation, we set constant values for the particles in the initial state and retained these values during simulations; this is further explained in section 4.3.

A reformulation of the viscosity term in the two-fluid model
The viscosity term in equation (34) also requires improvement because this study focuses on the problem of liquid 4 He rotated by a cylindrical vessel along a single direction at a constant angular velocity. In this case, the fluid particles' angular momentum around their axes (often called 'spin angular momentum' as an analogy for the corresponding term in quantum mechanics) could contribute to the dynamics of the entire system. However, many SPH series, including those reported in the literature [39], do not consider the spin angular momentum of fluid particles. This is mainly because these previous works focus on governing equations that omit the rotational contribution by the spin angular momentum of particles. Although neglecting rotational terms does not cause a problem in many cases, it can result in severe issues in certain cases, particularly for a fluid containing insolvable matters such as vesicle suspensions.
In 1964, D. W. Condiff derived the Navier-Stokes equations in the Lagrangian form with spin angular momentum conservation from Cauchy's linear momentum principle, aiming to reproduce the tiny effect of rotating molecules [56]. The resulting equation is expressed as Here,h,h r , andx are the shear viscosity, rotational viscosity, and bulk viscosity, respectively.h r determines the magnitude of rotational forces and must satisfy the relationship¯h h   0 r [57]. The third term on the righthand side becomes 0 in the case of an incompressible SPH because of the condition ∇ · v = 0.
It can be understood that equation (39) is an extension of the ordinary Navier-Stokes equations; the terms on the right-hand side, except for the first term, converge toh v 2 as the parameterh r becomes 0 when ∇ · v = 0. In addition, K. Müller discretized equation (39) using standard SPH in his study on formulating angular momentum conservation in smoothed dissipative particle dynamics [57]. Given these previous studies, we propose introducing the same extension as equation (39). Namely, we extend η n ∇ 2 v in equation (34) to the sum of (¯¯) h h +  v r 2 and¯w h 2 r , which implies thath h = n . We then compute (¯¯) h h +  v r 2 using equation (14) and¯w h 2 r for the ith particle using K. Müller's model [57], as follows: where ω i and ω j are the angular velocities of the ith and jth particles, respectively. We set the ω of each particle to in the initial state, where d indicates the initial particle distance, v max represents the estimated maximum velocity, and C ω is a constant parameter determining the scale ratio to v max . The value of C ω lies between 0 and 1.0, depending on the target problems in the simulation tests.

Entropy density estimation
It can be said from equation (27) and equation (31) that each of q i s and q i n functions as a parameter that determines the degree of the effect of temperature gradient forces. The entropy density σ i is dominant in the behavior of q i s because the density ρ i fluctuates typically less than one percent compared with its averaged density in incompressible SPH [58,59]. Additionally, the behavior of σ corresponds to that of entropy S because the value of NM is constant; thus, q i s is proportional to the entropy S. Figure 3(a) shows the dependence of entropy S on temperature T, indicating that the effect of the temperature gradient force controlled by q i s on the superfluid component diminishes as the temperature T decreases, and it eventually reaches zero at absolute zero. Based on the result, we can say that the parameter q i s reflects well the characteristic that superfluids have almost zero or a very small temperature gradient, corresponding to extremely high heat conductivity.
Regarding the parameter q i n , the second sound velocity c 2 of equation (33) should be semi-empirically determined according to the literature, exemplified by [60,61]. In experiments, it is known that c 2 emerges when the system reaches temperatures lower than the critical temperature of approximately 2.17 K, increases as the cooling of the system increases, and finally reaches approximately 35 ms −1 after experiencing a sluggishness at approximately 20ms −1 [61]. Let us assume in this discussion that c 2 is constant. In this ideal case, the parameter χ becomes constant; thus, q i n becomes proportional to ( ) s -T i i 1 . Also, the behavior of σ corresponds to that of entropy S, similar to the case of q i s . After all, the behavior of q i n is determined by that of (TS) −1 . Figure 3(b) shows the dependence of (TS) −1 on the temperature T with a semi-logarithmic scale. The effect of q i n on fluid 4 He is confirmed to greatly increase as the temperature T i decreases. Remember that the parameter q i n is the discretized formula of equation (17), expressing the ratio of superfluid density to normal density, which generates the second sound wave. Figure 3(b) suggests that the effect of ρ s /ρ n increases as the temperature T decreases, on the premise that c 2 is constant; the effect is more amplified in actual cases because the velocity c 2 itself increases, as mentioned above [61].
As discussed above, both parameters q i s and q i n are confirmed to reflect the phenomenological behavior of superfluid 4 He qualitatively. However, it is also found that q i n is too sensitive to be used for practical applications. As a remedy, we give the parameters of ρ s and ρ n as constant initial parameters; we compute equation (31) from equation (27) multiplied by − ρ s /ρ n . We also compute the entropy σ i in the parameter q i s from equation Here, the reason we introduce C e is because equation (21) does not consider the effect of interactions among particles in its derivation. In that sense, equation (21) can be said to express a relationship between entropy and temperature in a quantum ideal gas. Hence, it is necessary to consider the volume ratio of gas to a liquid. In this article, we estimate the value of C e as the volume ratio of gas to liquid, each of which has the same density under a certain pressure (C e ≈ V liquid /V air ), as the first approximation level.

The variables of the system
The variables for the ith particle are given as a set of ( ) i . P i is calculated from ρ i by solving a one-to-one thermodynamic correspondence between P i and ρ i given by the Tait equation of state [62]. Meanwhile, ρ i and σ i are calculated using equation (36) and equation (41), respectively. m i and T i keep their initial values in this paper. ( ) v s i and ( ) v n i are not necessarily distinguished in implementations because each particle only takes either a superfluid or normal fluid state at a time. Thus, the variable of the ith particle can be reduced to the set of (r The parameter C L includes one time-dependent parameter of the vortex line density L and three other constant values: the superfluid density ρ s , a friction coefficient α, and the quantum of circulation κ, where L is calculated from equation (9) at every time step. In this article, the parameter β v of equation (9) is calculated according to table 1 and table 3 in the literature [63]. To summarize, the constant parameters characterizing the system are given as a set of (¯¯)   right; the small secondary vortices at the bottom of the simulation domain emerge when Re = 100 and become larger when Re = 1,000, similar to the benchmark solutions [64].
The reason we examined the two cases of Re = 100 and Re = 1,000 is as follows. In section 5.3, we present the numerical experiments we carried out by rotating 4 He; this was done by setting the outer diameter of the circular vessel to 0.2 cm, the temperature T to 1.6 K, and by setting the angular velocity of rotating around the axis of the cylinder to 5 rad/s. It is difficult to estimate the Reynolds number of the total fluid system because only the normal fluid components have viscosity. Hence, we considered ideal cases in which the simulation domain is only filled with normal fluid. In this case, the Reynolds number is approximately 176.6 when setting the abovementioned diameter, temperature, angular velocity, and kinetic viscosity ν n (=μ n /ρ n ) at T = 1.6 K to the parameters v c , L c , and ν c . We also estimated Re to be approximately 1064.9 when evaluating ν c by the parameter ν, which is equal to μ n /ρ. We then take up the closest cases of Re = 100 and Re = 1,000 to these estimations, among those that can be compared with the benchmark solutions.

Numerical analyses
In this section, we demonstrate our method by performing three major numerical tests: the Rayleigh-Taylor instability, entropy wave propagation analyses, and rotating cylinder simulations.
In section 5.1, we carry out a Rayleigh-Taylor instability analysis to obtain a fluid-mechanical description of the interactions between two components when simultaneously solving both equations of motion of the twofluid model. Recently, the linear growth of the mushroom structure similar to that in classical fluids has been reported to be observed in the direct simulation of the Gross-Pitaevskii equation in the case of two-component  [65]. Because of the similarity of the systems, we expect to observe a similar type of instability in our case. We examine whether we can find the emergence of the linear growth of the mushroom structure on the condition of a separate state of superfluid and normal fluid. Meanwhile, simulating wave propagation is challenging for the SPH, particularly in the case of an explicit time-integrating scheme. We discuss the applicability of SPH to entropy wave propagation problems by examining the effect of density fluctuations on the error in calculations with a hypothetical case in section 5.2.
Furthermore, we perform rotating cylinder simulations, as discussed in section 5.3. In this test, the superfluid is rotated by a cylindrical vessel in a single direction at a constant angular velocity. A similar setup has been reported in previous experimental studies [3,66]; these studies noted interesting phenomena that were different from those in classical fluid cases. Specifically, multiple quantum vortices parallel to the cylinder axis emerged and rotated around their respective axes in the same direction as the cylindrical vessel. These vortices were spontaneously aligned, forming a lattice, the so-called 'quantum lattice,' and revolving around the cylindrical axis while maintaining constant relative positions with each other, similar to rigid body rotation.
Here we investigate the possibility of reproducing several quantum mechanical phenomena observed in liquid 4 He even when solving the phenomenological governing equations of liquid 4 He using our SPH method.

Rayleigh-Taylor instability analysis
We set the simulation domain of (L x , L y ) to be (100 μm, 200 μm) and place the superfluid particles at y > 0.5L y and the normal fluid particles at y < 0.5L y . We set the resolution of (N px , N py ) to be (240, 480), where N px and N py are the number of particles arranged in each direction. The interface between the two phases is disturbed by the function ( ( ) ) p = y L x L 1.0 0.15 sin 2.0 x x , and the value of the gravity constant of 9.8 m/s 2 is set in the y-direction, in a manner similar to the case in a classical fluid [67]. We arrange multiple layers of frozen particles as walls, alongside the boundary of the simulation domain. No-slip conditions are imposed on the wall particles as boundary conditions. This is acceptable because both normal and superfluid components have frictional forces against the normal fluid particles forming the wall.
The reference temperature T 0 is set to be 1.6 K. The parameters (ρ s , ρ n , α, β v , κ, η n ) are determined by reference to the values at T 0 = 1.6 K in table 1 and table 2 in the literature [63]. The resulting density ratio of ρ s /ρ n becomes 5.024. Meanwhile, we set the vortex line density L 0 to be 1 × 10 6 cm −1 as an initial state, in reference to the literature [26]. The parameters (h,h r , C ω ) are set to be (η n , 0, 0) in this test. We estimate the dimensionless parameter C e to be 1.428 × 10 −3 considering the volume ratio of V liquid /V air at 1.013 25 bar [68]. The value of NM in equation (22) is set equivalent to r L n x 2 because only the normal fluid components carry entropy. Figure 5 shows snapshots of the Rayleigh-Taylor instability simulation. The red part and blue part in figure 5(a) respectively represent the superfluid particles and normal particles. As a result, we confirmed the emergence of the characteristic mushroom structure in the linear growth regime. The mushroom structure shows its growth in the slanting direction because of wall boundary conditions; this is the same as the case for a classical fluid [69]. The time from the start to the final snapshot is approximately 7.49 ms in physical time, which corresponds to approximately 2.34 in the time scaled by the factor ( ) g L x . Figure 5(b) represents the magnitude of the viscosity term calculated using equation (34); the magnitude of viscosity is confirmed to be zero in normal fluid components. Figure 5(c) shows the magnitude of the velocity of particles for reference. Consequently, we obtained the description of the interactions between two components when simultaneously solving both equations of motion of the two-fluid model; it was observed to show the linear growth of the mushroom structure unique to the Rayleigh-Taylor instability problems, similar to ordinary fluids.

Effect of density fluctuations on the numerical analysis of the entropy wave equation
Let us examine the effect of density fluctuations on the calculation error of wave propagation problems with a hypothetical case. The wave equation is given as where c is the speed of sound. We impose closed boundary conditions: x In this case, the analytical solution of equation (42) is given as follows [70]: Here, C m,n and C m n , * correspond to the coefficients of the sine terms and cosine terms in the double Fourier expansion of equation (45), respectively. In this benchmark test, we discretize the left-hand side of equation (42) using the second-order central difference method [71] while replacing the right-hand side with the Laplacian operators of SPH or the finite-difference method (FDM); we can focus on the problem of numerical error in the spatial direction.
We set the simulation domain of (L x , L y ) to be (1.0, 1.0) and set the resolution of (N px , N py ) to be (400, 400). We take two cases, setting the speed of sound c to be 135 ms −1 and 20 ms −1 , assuming two representative cases of the second sound waves of 4 He in the low-temperature region, as mentioned in section 4.3. We set the pair of (m, n) to be (4, 3) and set the initial velocity of each particle to be 0 everywhere in the simulation domain. In this case, the values of C m,n and C m n , * become 0 and 1. Figure 6 shows snapshots of the analytical solution of the targeted problem during the time from 0 s to the half-period of time.
In this test, we set the same mass for all the particles. In this case, only the non-uniformity of particle arrangements contributes to the fluid density fluctuation because of the weighted calculation according to the positions of particles. Let us introduce the initial distance between two particles, l, and a distance parameter dl. Consider the case where we isotropically arrange particles with an interval of l in each direction and shift the particles by dl in an arbitrary direction. In this case, the averaged small volume ΔV 0 is l 2 , and the maximum one ΔV is (l + dl) 2 . We define a volume fluctuation rate δD as ((ΔV ) − ΔV 0 )/ΔV 0 × 100. Because of the uniformity of mass and the relationship ΔV = m/ρ, it is possible for us to regard δD as a criterion for density fluctuation. Additionally, the physical time is set to 0.03 s, roughly assuming the time from the emergence of the second wave  to the decay of the vortex line density [26]. In addition, we introduce a stabilization technique presented in the MPS (moving particle semi-implicit) [72] to suppress the increase in variance caused by repeated calculations (see appendix). We then compare the simulation results of SPH with that of FDM and the analytical solutions given by equation (45). Figure 7 shows the variation of S at a certain point during the simulation with C s = 135, in the respective cases. The solid line indicates the analytical solution, and the solid line with the cross symbol shows the simulation result of FDM. The remaining seven lines with symbols indicate the simulation results of SPH with different values of δD between 0 and 3 per cent. All results were confirmed to show excellent agreement at this viewing scale. Figure 8(a) shows further analyses of the results in figure 7, demonstrating the dependence of the relative error on the parameter δD in the respective cases. Here, the relative error E r is defined as the average of |S i − S a |/|S a + A| × 100 over all the particles, where S i is the value of S of the ith particle and S a is the analytical value of the point on the particle. A is the amplitude of the oscillation, which is equivalent to 1 in this case and is added to the denominator to avoid division by zero. The solid curve indicates the fitting curves using the least squares method. The relative error δD in the case of SPH was found to be a minimum at approximately δD = 0.5; subsequently, it increased as δD increased by a power of two. Meanwhile, figure 8(b) shows the simulation results in the case where C s = 20. The dependence of E r on the parameter δD shows the same increase as in the case where C s = 135.
The results in figure 8 demonstrate the applicability of SPH to wave propagation problems of liquid 4 He, even though the relative error E r depends on the parameter δD by a power of two. This occurs because the relative error is small enough to be less than 0.25 per cent, even in the case where δD = 2. Here, δD is known to be approximately 0.5, in usual cases. Thus, the relative error is approximately 0.012 per cent for C s = 20, and 0.018 per cent for C s = 135. Through the discussion in this section, it is shown that we can perform wave propagation analyses within the range of these permissible errors.  The reason why δD is estimated to be 0.5 is explained as follows. The Mach number, the ratio of the maximum velocity to the speed of sound, is denoted as M a . Given that the temperature is constant, the relative error of density is proportional to the square of M a , which is described as ( ) r r r -= M 0.5 a 0 0 2 as in [73]. Because of the uniformity of mass, we can obtain the relationship ( ) where ΔV 0 = m/ρ 0 and ΔV = m/ρ. Because ΔV 0 is l 2 and ΔV is (l + dl) 2 , a simple calculation leads to the following relationship between δD and M a : Figure 9 shows the relationship between M a and δD described by equation (48). In the simulation of incompressible SPH, M a is given as an input parameter, which is typically set to 0.1. The value of δD at this time is approximately 0.5, as shown in figure 9.

Rotating cylinder simulations
The outer diameter of the circular vessel was set to be 0.2 cm. We set the angular velocity of rotating around the cylinder axis to 5 rad/s in the counterclockwise direction. The resolution of (N px , N py ) was set to be (500, 500). We set the temperature T 0 to 1.6 K and determined the parameters (ρ s , ρ n , α, β v , κ, η n ) by referring to the values at T 0 = 1.6 K listed in a previous report [63]. Meanwhile, the parameters (h,h r , C ω ) are set to be (η n , η n , 0.01) in this simulation test. The wall particles, made of normal fluid particles, are arranged along the inner circumference of the vessel. No-slip conditions are imposed on the wall particles as boundary conditions. This is acceptable because both normal and superfluid components have frictional forces against the normal fluid particles forming the wall. In the initial state, superfluid or normal fluid particles are randomly generated in proportion to the density ratio and are distributed inside the vessel. Figure 10 shows an image of the rotating cylinder simulation at t = 0.5 s. In this figure, (a) represents the superfluid particles in red and normal fluid particles in blue; (b) shows a view colored by the magnitude of velocities of the particles; and (c) shows an enlarged view of the center of (b) colored by the magnitude of rotational forces given by equation (40). A movie of the complete simulation is provided as supplemental material, and its online version is available at the link in the caption of figure 10. This movie has a viewing time of 36 s, which corresponds to 5 s in terms of physical time.
Interesting observations were noted. Moments after commencing the simulation, normal fluid particles gathered and formed clusters, which changed to vortices revolving around their individual axes, with their pores filled with normal components. These vortices changed their locations by interacting with each other, while occasionally connecting with smaller vortices or absorbing the normal particles that failed to form vortices. The vortices located near walls are affected by the forced rotation of the outer vessel. Because the effective range of forced rotation continuously increases in space as physical time increases, the vortices gradually rotate around the circular axis, independent from their own rotations. Figure 11 shows the distributions of the velocity of all particles at (a) t = 0.5 s and (b) t = 5 s. The horizontal line indicates the distance between a particle and the center of the circular axis, and the vertical line represents particle velocity. From figure 11(a), it is confirmed that the velocities of particles located in the vicinity of the wall are considerably greater than those of particles in the inner area because of the effect of the forced rotation. The velocities of normal fluid particles in the inner area at L < 0.06 cm are confirmed to fluctuate more widely than those of the superfluid particles, because the normal particles form the vortices and rotate around the axes of the vortices. As time passes, the vessel's forced rotation gradually affects the inner part of the vessel. Figure 11(b) Figure 9. Relationship between the mach number Ma and the parameter δD.
shows the velocity distribution for t = 5 s. The white dotted line represents the fitting line obtained by the least squares method. It is understood from (b) that all the vortices revolve around the circular vessel'axis with almost the same angular velocity of |ω| = 5 rad/s and that their averaged velocities follow a relation of L|ω|. This result suggests that the vortices exhibit rigid-body rotation.
As mentioned above, previous experiments have reported that multiple quantum vortices parallel to the cylinder axis emerged and rotated around their respective axes in the same direction as the cylindrical vessel. These vortices are arranged such that they form a lattice, called the 'quantum lattice.' They revolve around the circular axis while maintaining constant relative positions with each other, similar to rigid body rotation. Through comparisons of our simulation results and the experimental results, it can be said that our current model succeeds in reproducing the first and third phenomena, i.e., the emergence of multiple independent vortices parallel to the circular axis and that of 'rigid-body rotation. ' Unfortunately, our current model needs improvements in terms of reproducing quantum lattices. It would be necessary to generate uniform vortices to ensure that the repulsive forces among them balance each other, thereby forming a lattice. Figure 12 shows a histogram of the pore areas in the vortices for t = 5 s, which is obtained from the resulting image using ImageJ [74][75][76]. pixel 2 gives the unit of area. The histogram visually clarifies the discrepancy between the target phenomena and our simulations; the histogram should tend to a delta function as the vortices become uniform. A promising strategy to realize such convergence is to incorporate the 'quantization of circulation' in SPH formalism to ensure that each vortex has the same size and energy; this needs to be studied in future works.
Additionally, we performed the same simulation test for different values of C ω . It was confirmed that the magnitude of C ω determines the balance between the strongness of vortices and that of forced rotation. At Figure 10. An image of the rotating cylinder simulation at t = 0.5 s: (a) a view representing the superfluid particles in red and normal fluid particles in blue; (b) a view colored by the magnitude of velocities of the particles; and (c) an enlarged view of the center of (b) colored by the magnitude of rotational forces given by equation (40). A movie of the complete simulation is provided as supplemental material, and its online version is available at https://www.satoritsuzuki.org/video-preprintarxiv-200411052. approximately C ω > 0.015, the vortices swell by connecting with or absorbing each other. Consequently, several enlarged vortices continue moving inside the simulation area, almost regardless of the vessel's forced rotation. In contrast, the vortices hardly rotate around their own axes at approximately C ω < 0.01.
It has been assumed that only a single vortex emerges when addressing similar problems without considering quantum mechanics. However, our simulation results show that the two major phenomena of rotating liquid 4 He, i.e., the emergence of multiple independent vortices parallel to the circular axis and that of the so-called 'rigid-body rotation,' can be reproduced by solving the two-fluid model, which is the phenomenological governing equation of liquid 4 He, by using SPH. The rotational cylinder problem is one of the few problems exhibiting different phenomena depending on basic mechanics or quantum mechanics. It is meaningful that our model succeeded in reproducing the phenomena observed in quantum cases, even though we solve the phenomenological governing equations of liquid 4 He.

Conclusions
In this study, we theoretically derived a finite particle approximation of the two-fluid model of superfluid 4 He using smoothed particle hydrodynamics. To solve the governing equations, we introduced an auxiliary equation representing the microscopic relationship between entropy and temperature, which was derived from quantum statistical mechanics in elementary particle physics. We then formulated a new SPH model that simultaneously solves both equations of motion of the two-fluid model. In particular, we presented a reformulation of the viscosity term in the two-fluid model to conserve the angular momentum of the fluid particles around their axes. Through numerical tests, several sophisticated techniques have been introduced. Specifically, we utilized a Riemann solver among neighboring particles to stabilize the simulation as an alternative to artificial viscosity. We demonstrated a Rayleigh-Taylor instability simulation, confirming the emergence of the characteristic mushroom structure in the linear growth regime. We investigated the applicability of SPH to entropy wave propagation analyses by examining the effect of density fluctuations on the errors in calculations with a benchmark case. The simulation results show good agreement with those of the FDM and the analytical solution.
Furthermore, we performed a simulation of rotating liquid 4 He. It was found that two major phenomena of rotating liquid 4 He-the emergence of multiple independent vortices parallel to the circular axis and that of the so-called 'rigid-body rotation'-can be reproduced by solving the two-fluid model, which is the phenomenological governing equation, by using SPH. This finding is interesting because it was previously assumed that only a single vortex emerges when addressing similar problems without considering quantum mechanics. Our further analysis found that the emergence of multiple independent vortices can be realized by the aforementioned reformulation of the viscosity term to conserve the angular momentum of the particles around their axes. It is meaningful that our model succeeded in reproducing the phenomena observed in quantum cases, even though we solve the phenomenological governing equations of liquid 4 He.
Consequently, our simulation model has been confirmed to reproduce the phenomenological characteristics of superfluid 4 He. We have, thus, pioneered a new paradigm of Lagrangian particle mechanics, which will hopefully stimulate further numerical studies on superfluids in low-temperature physics.