Our recent study suggested that a fully classical mechanical approximation of the two-fluid model of superfluid helium-4 based on smoothed-particle hydrodynamics (SPH) is equivalent to solving a many-body quantum mechanical equation under specific conditions. This study further verifies the existence of this equivalence. First, we derived the SPH form of the motion equation for the superfluid component of the two-fluid model, i.e., the motion equation driven by the chemical potential gradient obtained using the Gibbs–Duhem equation. We then derived the SPH form of the motion equation for condensates based on the Gross–Pitaevskii theory, i.e., the motion equation driven by the chemical potential gradient obtained from the Schrödinger equation of interacting bosons. Following this, we compared the two discretized equations. Consequently, we discovered that a condition maintaining zero internal energy for each fluid particle ensures the equivalence of the equations when the quantum pressure is negligible. Moreover, their equivalence holds even when the quantum pressure is non-negligible if the quantum pressure gradient force equals the mutual friction force. A zero internal energy indicates the thermodynamic ground state, which includes an elementary excitation state. Therefore, the condition can be sufficiently satisfied when the velocities of fluid particles do not exceed the Landau critical velocity, which is not a stringent condition for simulations with a characteristic velocity of a few $cm\xb7s\u22121$ in a laboratory system. Based on the above, we performed a simulation of rotating liquid helium-4 and succeeded in generating a vortex lattice with quantized circulation, known as a quantum lattice.

## I. INTRODUCTION

Since the beginning of the 20th century, a significant decrease in the viscosity of liquid helium-4 in cryogenic regions (approximately 2.1 K) has attracted attention of researchers in the field of low-temperature physics. The viscosity loss of liquid helium-4 as a continuum corresponds microscopically to the loss of molecular viscosity, which results from the van der Waals forces between helium atoms. This loss is even more microscopically attributed to the innate nature of atoms, i.e., the change in the energy states of helium atoms manifested in the cryogenic temperature range. Notably, the film flow is a phenomenon observed when liquid helium-4 creeps out of its container as it loses its viscosity and becomes a superfluid. The fountain effect, caused by a superleak of helium atoms passing through a porous medium, such as a plaster, is also well known. More precisely, the “fountain effect” refers to the thermodynamic effect observed in a vessel immersed in a reservoir, which is equipped with a heater, an opening at the top, and a capillary filter at the bottom. This arrangement allows the flow of a superfluid into the heater and causes the normal fluid generated by the heater to flow out through the opening. Because the viscosity loss properties of helium atoms in this temperature range enable them to penetrate such capillary filters, this penetration phenomenon is referred to as a “superleak.”

Understanding the dynamics governing the behavior of superfluid helium-4 is not only an outstanding academic achievement in low-temperature physics but also expected to clarify various related phenomena owing to the similarities with the dynamics of the phase transition of helium-4. Some related examples are liquid crystals^{1–4} and superconductivity^{5,6} in condensed matter physics, Hawking radiation^{7–10} and neutron star crusts^{11} in astrophysics, and turbulence.^{12–14} In addition to the value of such a scientific study, realizing robust and safe control of superfluid helium-4 in the bulk state by predicting its complex behavior under various conditions via large-scale fluid simulations is of paramount importance in the development of cryogenic cooling systems for x-ray observation satellites^{15,16} and space telescopes^{17,18} for high-energy astronomy; these simulations may be based on classical fluid dynamics while being capable of estimating the effects of macroscopic quantum phenomena. Therefore, a direct numerical analysis of superfluid helium-4 in the bulk state on scales ranging from centimeters to meters can be expected to facilitate the development of safe and robust cryogenic cooling systems and significantly enhance the engineering applicability of quantum liquids.

However, no method for simulating the behavior of such large-scale superfluid helium-4 has yet been established. Generally, methods that can simulate superfluid helium-4 can be categorized into four major types. The first involves treating superfluid helium as a quantum many-body system with a Bose–Einstein condensate^{19–22} and directly solving the Schrödinger equation. The second is the vortex filament model (VFM), which solves the dynamics of superfluid fields using a Lagrangian approach.^{23–27} The third method is the density functional theory (DFT), which is also a well-established method for obtaining the density states.^{28–30} Examples of the applications of the DFT to liquid helium-4 include studies on nanodroplets^{31,32} and electron bubble states.^{33,34} These three methods are quantum mechanical models that describe phenomena on the order of nanometers to micrometers. However, using them in the simulations of macroscopic cryogenic liquid helium-4 would require an excess of tens of billions of atomic-sized analytic particles, which is unrealistic even with the latest supercomputers. In contrast, the fourth method is a phenomenological model proposed by Landau^{35} and Tisza,^{36} known as the two-fluid model. It describes liquid helium-4 based on two components: a superfluid component that is an inviscid incompressible fluid and a normal fluid component that follows the incompressible Navier–Stokes equation; the model describes the system as a mixture or superposition of these components. However, this is a classical fluid approximation, and it cannot reproduce macroscopic quantum phenomena, such as a vortex lattice in rotating liquid helium-4.

Recently, an improved two-fluid model, including the conservation of angular momentum, was developed in our previous study.^{37,38} This model has attracted attention as a new approximation model that balances classical and quantum mechanical descriptions to address the drawbacks of existing approaches. Our previous study^{37} reformulated the motion equation of the normal fluid component of the two-fluid model to include a term accounting for the angular momentum conservation of particles around their axes. Let us refer to this rotational angular momentum as the “spin angular momentum” analogous to the corresponding term in quantum mechanics. As will be discussed later, our previous study^{37} treated the fluid forces of both components as a classical mechanical approximation and discretized the reformulated two-fluid model using smoothed-particle hydrodynamics (SPH), which is a well-established Lagrangian particle approximation model for flow problems originally developed in astrophysics.^{39} Interestingly, we observed the emergence of parallel spinning vortices presenting rigid-body rotation in the numerical simulation of rotating liquid helium-4. In addition, our subsequent study^{38} incorporated a vortex dynamics model into the reformulated two-fluid model, successfully capturing the phenomenon of vortex lattices via numerical simulations,^{38} also reported that the number of vortices generated in a system of rotating cylinders agreed with the theoretical solution based on Feynman's rule.^{40,41} However, it still included one free parameter that had to be optimized to determine the intensity of the angular velocity. A video of the simulation in Ref. 38 is available at https://arxiv.org/src/2105.03177v3/anc. These results indicate that a vortex lattice in rotating liquid helium-4 can be reproduced even by solving the two-fluid model based on a fully classical mechanical approximation that includes the fluid forces of both components. This is applicable under the condition that the viscosity is rederived to conserve the rotational angular momentum using the SPH and the vortex dynamics is incorporated into the system. In particular, a vortex lattice in rotating liquid helium-4 is deemed a purely quantum mechanical phenomenon. Thus, the results reported in Ref. 38 have gained significant scientific interest because they directly challenge the decade-long preconception that a fully classical fluid dynamics approach cannot reproduce quantum mechanical phenomena, such as vortex lattices in liquid helium-4.

Figure 1 presents a schematic of the improved two-fluid model developed in our previous study.^{37,38} The involved method has two essential features. The first is the reformulation of the viscosity term in the normal fluid component of the two-fluid model. Although the ordinary two-fluid model can describe the macroscopic dynamics of cryogenic liquid helium-4, angular momentum conservation, which is the essence of quantum mechanics, has not been explicitly formulated. In this regard, in a study on numerical simulations of mesoscale flows considering microfluidics, the Navier–Stokes equation with angular momentum conservation for the rotational motion of molecules constituting a polar fluid was formulated by Condiff and Dahler.^{42} Moreover, the conservation of angular velocities of fluid particles was explicitly rederived in the Lagrangian form. In addition, in recent studies, a method^{43} has been developed to discretize this rederived Navier–Stokes equation using smoothed-dissipative particle dynamics^{44} for the simulation of mesoscale flows, such as a blood flow, with heat dissipation.^{45,46} These studies led us to believe that applying the same technique to the discretization of the motion equation for the normal fluid component may enable angular momentum conservation in the two-fluid model, thereby making it feasible to reproduce the macroscopic quantum phenomena of liquid helium-4. The second feature is the finite particle approximation using Lagrangian particle mechanics, which is well suited for these types of scenarios, i.e., scenarios regarding fluid particles as coarse-grained helium atoms. It also facilitates spin angular momentum conservation of the fluid particles by maintaining the angular velocity constant for each spherical fluid particle, as illustrated in the lower left side of Fig. 1. In our previous studies,^{37,38} the SPH was employed to this end.

In addition, greater emphasis should be placed on our classical approximation detailed in Refs. 37 and 38, which includes the fluid forces of both components and solves their motion equations in a fully coupled manner, such as a multiphase flow in classical fluid dynamics. We will further describe this point. The original two-fluid model proposed by Landau treats the aforementioned two components independently; however, the existence of mutual friction forces was proposed by Gorter and Mellink.^{47} Following this, the coupling of the two components was gradually confirmed. In a recent related study on liquid helium-4, the superfluid component was solved using the VFM and coupled with the Navier–Stokes equation.^{23,25–27} In our previous studies, we directly coupled the Navier–Stokes equation not with the VFM but with the motion equation of the superfluid component in the classical two-fluid model, i.e., we employed the motion equation of an inviscid incompressible flow. For quantum mechanical corrections, we conserved the spin angular momentum to ensure close resemblance between the fluid system and a quantum mechanical system. In a broader sense, among these existing coupling approaches, our approach can be classified as one with a coarse-grained model.

Importantly, we observed the emergence of multiple spinning vortices forming a rigid-body lattice, despite solving the two-fluid model based on a fully classical approximation incorporating the fluid forces of both components. Therefore, regardless of the above, the following hypothesis can be suggested: discretizing the two-fluid model using particles in the SPH formalism can extract the nature of a multiparticle interacting system, and a fully classical mechanical approximation may be equivalent to solving a many-body quantum mechanical equation under specific conditions. This hypothesis is valid only for large-scale problems because the fully classical approximation disregards several laws of quantum mechanics. In our previous studies, the SPH discretization of the two-fluid model with angular momentum conservation and its application to large-scale problems of rotating liquid helium-4 made it possible to observe the phenomena of vortex lattices. However, no study has provided theoretical evidence supporting the equivalence of the microscopic equation of motion of a many-body quantum system and the phenomenological motion equation of the superfluid component of the two-fluid model in the SPH formalism.

To address this gap, this study aims to demonstrate the existence of this equivalence. First, we derive the SPH form of the motion equation of the superfluid component of the two-fluid model, i.e., the motion equation driven by the gradient of the chemical potential obtained using the Gibbs–Duhem equation. Then, we derive the SPH form of the motion equation for condensates from the Gross–Pitaevskii (GP) theory, i.e., the motion equation driven by the gradient of the chemical potential obtained from the Schrödinger equation of interacting bosons. We then compare these two discretized motion equations in their SPH forms derived independently from microscopic and macroscopic perspectives to identify a condition wherein both become equivalent. Notably, we report that maintaining the internal energy zero for each fluid particle ensures the above equivalence. We also discuss the implications of this finding in terms of real-world cases and detail its incorporation into liquid helium-4 simulations as a stepping stone toward accurate future numerical reproduction of macroscopic quantum phenomena, such as film flows and fountain phenomena.

## II. BRIEF REVIEW OF EXISTING THEORIES

### A. Gross–Pitaevskii theory

The time-dependent many-body Schrödinger equation and its quantum Hamiltonian *H* are as follows:^{48}

where $\psi (r1,r2,\u2026,rN,t)$ represents a many-body wavefunction, *N* denotes the number of particles, *m _{q}* denotes the mass of a particle, $\u210f$ denotes the reduced Planck's constant, $U(ri,t)$ represents the external potential, and $V(ri\u2212rj)$ represents the mutual interaction potential between the

*i*th and

*j*th particles.

For *N* identical interacting bosons, the wavefunction $\psi (r1r2,\u2026,rN,t)$ is symmetric with respect to the exchange between the *i*th and *j*th particles, as follows:

The Bose–Einstein condensate^{19} assumes that all particles occupy the same single-particle state. Thus, the many-body wavefunction $\psi (r1,r2,\u2026,rN,t)$ can be decomposed into a product of single-particle wavefunctions $\varphi (ri,t)$, as follows:^{49–51}

where $\varphi (ri,t)$ satisfies the normalization condition as follows:

The condensate wavefunction is described as follows:

Based on the definition given in Eq. (6) and the normalization condition given in Eq. (5), we can obtain $\u222b|\Psi |2dr=N$. Therefore, the condensate wavefunction $\Psi $ can be expressed as a function of the condensate density $n(r,t)$, as follows:^{52}

where *i* denotes an imaginary unit and $\theta (r,t)$ represents the phase of the condensate wavefunction.

Substituting the Lagrangian function of the Hamiltonian given in Eq. (2) into the Euler–Lagrangian equation to satisfy the least-action principle yields the following time-dependent GP equation for the condensed matter wavefunction:^{53,54}

By substituting Eq. (7) into Eq. (8) and comparing the real and imaginary parts of both sides of Eq. (8), we can obtain the following series of equations for the condensate density $n(r,t)$:^{52,54}

where we can deduce the relationship between the velocity of the superfluid component and the phase of the condensate wavefunction as follows:

### B. Smoothed-particle hydrodynamics

In SPH, a discrete physical quantity $\phi $ is represented as a continuous quantity using the Dirac delta function, *δ*, which is further approximated using a distribution function *W*, known as the kernel function,^{55} as follows:

where kernel function *W* is selected to satisfy the following conditions:

A straightforward example of *W* satisfying the given conditions in Eqs. (14)–(16) is the Gaussian kernel^{39} expressed as $W(r\u2212r\xb4)=C/hd\u2009exp[\u2212|r\u2212r\xb4|2/h2]$, where *C* denotes a normalization constant, *h* denotes the kernel radius, and *d* represents the dimension. However, in the simulations of incompressible flows, a polynomial function^{56,57} that converges to zero at a distance of *kh* is selected as the kernel function because it satisfies the normalization condition given in Eq. (16) by simply integrating over a distance *kh* without having to integrate over the entire domain. In this study, a cubic spline kernel function^{58} is used in the numerical experiments described in Sec. IV.

The gradient of $\phi $ is expressed as follows:

Here, Eq. (18) is derived from Eq. (17) using Gauss's divergence theorem, after representing Eq. (17) in the form of integration by parts.^{55}

Let us consider a case wherein the integral domain Ω is divided into *N _{p}* small volumes $\Delta Vi$ ($i=1,2,\u2026,Np$). Based on the summation approximation, Eqs. (13) and (18) can be expressed in discrete forms using mass

*m*, density

_{i}*ρ*, and small volume $\Delta Vi$ and their relationship $\Delta Vi=mi/\rho i$, as follows:

_{i}where $Wij=W(|ri\u2212rj|,h)$ and $\u2207Wij$ denote the gradient of *W _{ij}*. The remainder of our analysis only requires the expression of the gradient of $\phi (ri)$ in Eq. (20). Additional details on improved SPH operators for the gradient, basic and improved SPH operators for rotation and divergence, and several numerical techniques necessary for stable SPH simulations can be found elsewhere.

^{38,43,55,59}

Before proceeding, we will attempt to clarify our concept of the SPH discretization of liquid helium-4 using Fig. 2 and explain some of the terms referred in the subsequent sections. First, as already mentioned, we apply the classical Lagrangian picture, as depicted in Fig. 2(b), to the helium-4 continuum, which is presented in Fig. 2(a). Specifically, we assume a scenario wherein two types of classical fluid particles are mixed: normal fluid particles following the Navier–Stokes equation and inviscid fluid particles following the momentum equation of a classical inviscid flow. Importantly, when discretizing a continuum using SPH, fluid particles are assumed to correspond to the discretization points of the continuum. Figure 2(c) illustrates an example of the application of the Voronoi tessellation, which is commonly used for the spatial discretization of classical incompressible fluids, to a domain of the helium-4 continuum. The gray dashed lines represent the boundaries of the Voronoi diagram. The black circle represents the neighboring area of the *i*th particle or the discretization point in the continuum. The bottom part of Fig. 2(c) depicts the particle distribution within the corresponding kernel function of the *i*th particle in the SPH calculation.

All fluid particles in this calculation are classical fluid particles [in fact, if these fluid particles are to obey quantum mechanics, they cannot be unambiguously labeled, as shown in Figs. 2(b) and 2(c), because identical particles cannot be distinguished in quantum mechanics]. Thus, in this paper, “normal fluid particles” or “superfluid particles” refer to classical fluid particles unless otherwise stated. It should be emphasized that superfluid particles refer to classical fluid particles that obey the momentum equation of an inviscid flow. Moreover, note that only in Sec. II A, the term “particle” implies an atom; otherwise, it implies a classical fluid particle. As mentioned, the two types of classical fluid particles (normal and inviscid) denote the discretization points of the continuum, as depicted in Fig. 2(c). Specifically, they are virtual particles. This concept is based on several recent mathematical studies demonstrating that SPH is a generalized particle method, wherein each fluid particle does not serve as a physically significant particle but simply acts as a discretization point in a Voronoi cell finite element.^{60–62} In our studies, we adopted this concept and postulated that these two types of virtual particles have no physical significance other than as discretization points, or fragments, of a continuous two-phase flow. Nevertheless, future research may provide some quantum physical interpretation of these virtual particles.

Because the two types of fluid particles maintain their volumes, both contribute to the macroscopic fluid phenomena uniquely observed as a whole. In summary, our model adopts a type of one-fluid system, wherein the fluid phenomena of the entire system are focused upon, whereas the motions of individual virtual particles have no physical significance from a microscopic viewpoint. The concept of a one-fluid system for liquid helium-4 is not very prominent; however, it has a long history.^{63} Previous studies reviewed in Ref. 63 described the dynamics of liquid helium-4 from a microscopic viewpoint, primarily based on quantum mechanics. Our studies differ from these studies in that we attempted to describe identical dynamics in terms of the finite particle approximation of a continuum. Importantly, the concept of our virtual fluid particles is different from the two-component concept of the ordinary two-fluid model. This discussion highlights that the motions of individual virtual fluid particles do not represent the behavior of a normal fluid or superfluid component of the ordinary two-fluid model. We reemphasize the following. Our particle approximation is based on the concept of a generalized particle method wherein a fluid particle only serves as a discretization point in a continuum. Specifically, all fluid particles are classical fluid particles that have no physical significance from a microscopic viewpoint.

The idea indicating that superfluid particles retain their volumes may be innovative. However, experiments have demonstrated that liquid helium-4 condensates produce only approximately 13% of the total atoms even at approximately absolute zero.^{64} Therefore, we do not observe a massive loss of liquid helium volume in the cryogenic range. Our model assuming a constant volume of inviscid fluid particles is consistent with these facts and has a certain validity from a phenomenological standpoint. Therefore, only the fluid phenomena of the system as a whole have physical significance. Because all particles are classical fluid particles, we refer to our approach as a “fully classical mechanical approximation.” By contrast, in the conventional two-fluid model, the macroscopic spatial domain is occupied by the normal fluid components, which microscopically interact with quantum vortices.^{24,27,65} Therefore, the computational results of fluid simulations based on our model, which employs a one-fluid system, must be compared with those of normal fluid component simulations of the two-fluid model based on conventional methods, which employ a two-fluid system. This is important when comparing the simulation results of this study with those of previous studies cited in Sec. IV.

In addition, the mutual friction forces are added to the system of equations, as depicted in Fig. 1, in our actual calculations. However, in the simulations based on our model, the magnitudes of the mutual friction forces are sufficiently smaller compared to those of the fluid forces, such as pressure gradient forces, when using the microscopic models of mutual friction forces reported in Ref. 66. Nevertheless, we cannot deny the possibility that these terms may affect the system as disturbances in certain scenarios, even in our scheme, which is based on classical fluid mechanics. Furthermore, considering the future extension of our model to microscopic phenomena in the regime governed by quantum mechanics, retaining these terms at this stage may be useful for formal consistency with the general two-fluid model, wherein a mutual friction term is introduced. Therefore, we retained these mutual friction force terms in the SPH model in our series of studies.

Regardless of the above, in the numerical simulations conducted in our previous study using the SPH model, we observed the following phenomena consistent with the phenomenological definition of a vortex lattice: (1) the emergence of multiple vortices spinning in the same direction and (2) their collective rotation around the center of a cylindrical container at a constant speed, i.e., rigid-body rotation. Nevertheless, to prove that the vortex lattice observed in our previous study corresponds to a real-world vortex lattice, it is necessary to theoretically demonstrate the validity of coupling the two motion equations of the two-fluid model, which was experimentally used in our previous studies.^{37,38} In the subsequent sections, we present that the motion equation of an inviscid flow and the nonlinear Schrödinger equation of a bosonic system in SPH forms are equivalent if certain conditions are satisfied. Coupling the motion equation of an inviscid flow and the Navier–Stokes equation is seemingly inadequate to capture the dynamics of quantum fluid phenomena because both these equations are classical equations of motion. However, if we prove that the motion equation of an inviscid flow is equivalent to the GP equation in the SPH form under appropriate conditions, we can claim that our coupling method is equivalent to the method that couples the GP equation with the Navier–Stokes equation in such cases. Therefore, demonstrating the equivalence of the nonlinear Schrödinger equation and the motion equation of an inviscid flow in their SPH forms strengthens the validity of coupling the two equations of the two-fluid model from a theoretical perspective.

## III. EXPRESSIONS OF THE GP EQUATION AND TWO-FLUID MODEL IN SPH FORMS

Let us consider a typical case wherein the mutual interaction potential *V* is given by a scalar multiple of the delta function *δ* as follows:

where *g* denotes the coupling constant, which represents the interaction strength.

Equation (21) represents a contact interaction model that assumes weak interactions between atoms. Therefore, it is essentially suitable for reproducing the behavior of dilute gases. Conversely, it is not suitable to reproduce the behavior of superfluid helium-4 based on strong interactions. Therefore, fluid equations derived from the GP equation are often used only to discuss their qualitative behavior. However, in SPH, the delta function is replaced by a kernel function, as in Eq. (13). Hence, Eq. (21) acts as a strong interaction model that determines the contribution of particles in the neighboring area with an effective radius. Thus, for SPH calculations, Eq. (21) can be expected to provide a relatively accurate approximation of superfluid helium-4 behavior.

By substituting Eq. (21) into Eq. (10), which represents the third term on the right-hand side of Eq. (10) in a discrete form using Eq. (12), and then operating $\u2207$ from the left on both sides of Eq. (10), we obtain the following:

Here, we use Eq. (11) to derive the left-hand side and first term on the right-hand side of Eq. (22). The expression provided in the square brackets in the fourth term on the right-hand side is known as the “quantum pressure,” hereinafter denoted as *P _{q}*. If the spatial variation in the condensate density profile, $n(r,t)$, is small, the quantum pressure can be neglected.

^{52}The chemical potential $\mu (r,t)$ is then expressed as follows:

^{51}

The set of second to fourth terms on the right-hand side of Eq. (22) corresponds to $\u2212\u2207\mu /mq$. Using the relationship $\u2207(vs2/2)=vs\xd7(\u2207\xd7vs)+(vs\xb7\u2207)vs$ and the condition of irrotational flow of the superfluid component $\u2207\xd7vs=0$, we can rewrite Eq. (22) as follows:

Using Eqs. (20) and (23), we can further rewrite Eq. (24) in a discretized form of SPH. The resulting equation for the *i*th fluid particle is as follows:

where the right-hand side can be obtained as follows:

Here, *U _{j}* and

*n*are the abbreviations of $U(rj,t)$ and $n(rj,t)$, respectively. $Pq(i)$ denotes the discretized form of the quantum pressure

_{j}*P*in Eq. (23), and it can be represented using a series of SPH operators. We denote the quantum pressure gradient force $\u2212\u2207Pq(i)/mq$ as $fq(ri,t)$ for subsequent discussions. In this manner, by approximating the delta function

_{q}*δ*with the kernel function

*W*, we derive the SPH form of the motion equation of condensates from the GP theory. Specifically, this represents the motion equation driven by the gradient of the chemical potential obtained from the Schrödinger equation of interacting bosons. In summary, we arrange Eq. (26) as follows:

Here, the subscript QM on the left-hand side symbolically indicates that the right-hand side is an equation derived from the microscopic equation of motion.

Let us recall that the phenomenological equation of motion of the superfluid component of the two-fluid model can be obtained by substituting the Gibbs–Duhem equation, i.e., $\u2207\mu =(V/N)\u2207P\u2212(S/N)\u2207T$, relation $\rho =(Nmq)/V$ for density, and relation $\sigma =S/(Nmq)$ for entropy density into Eq. (24). In these expressions, *V* denotes the volume of the entire system, *m _{q}* is the mass of a particle (an atom),

*N*is the number of atoms,

*S*is the entropy, and

*T*is the temperature. The resulting equation is as follows:

Using Eq. (20), we obtain another SPH form of the gradient of the chemical potential, *μ*, derived from the Gibbs–Duhem equation, as follows:

where *P _{j}* and

*T*denote the pressure and temperature of the

_{j}*j*th fluid particle, respectively. In summary, we obtain the following:

Here, the subscript TM on the left-hand side symbolically indicates that the right-hand side is an equation derived from a thermodynamic viewpoint. The upper part in Fig. 3 presents a schematic of the derivation of the two different SPH forms in Eqs. (27) and (30).

In the following, we compare the discretized motion equations given in Eqs. (27) and (30), which are separately derived from a microscopic and macroscopic perspective, to identify a condition under which they become equivalent. Let us focus on a case wherein the spatial variation in the condensate density, $n(r,t)$, is sufficiently small; therefore, the quantum pressure, $Pq(r,t)$, can be negligible. Because $fq(ri,t)$ in Eq. (27) vanishes owing to $Pq(r,t)=0$, we can derive a condition for the equivalence of Eqs. (27) and (30) by comparing each term, as follows:

Using the total density ($\rho =(Nmq)/V$), entropy density ($\sigma =S/(Nmq)$), and chemical potential ($\mu j=Uj+gnj+Pq(j)$) for the *j*th fluid particle, Eq. (31) can be represented as follows:

$Pq(j)$ of *μ _{j}* is zero in this case. Equation (32) represents a specific case of the thermodynamic Euler equation,

^{67},$ST\u2212VP+N\mu =U\u0303$, for the

*j*th fluid particle when its internal energy is zero, i.e., $U\u0303=0$. Here, the ̃ symbol is added to the internal energy to avoid confusion with the potential energy,

*U*. Equation (32) indicates that the equivalence of Eqs. (27) and (30) holds when the internal energy of each fluid particle is zero. In other words, this equivalence holds if the fluid system is in the ground state. Because this is a thermodynamic condition, the “ground state” includes an elementary excitation state [refer to Refs. 35 and 68–70 for a detailed explanation of the elementary excitation]. Herein, the ground state in quantum mechanics is referred to as the “real ground state.” We refer to the ground state derived explicitly from classical thermodynamics as the “thermodynamic ground state.” If we deduce that a system is in the thermodynamic ground state based on classical thermodynamics, we cannot conclude that the system is in the real ground state. This is because we cannot deny the possibility that the system is in an elementary excitation state. This information on elementary excitation can only be obtained by considering quantum statistical mechanics. In this regard, Eq. (32) represents a specific case of the thermodynamic Euler equation of each fluid particle when the internal energy is zero. This equation represents the classical thermodynamic relationship. Hence, Eq. (32) suggests that the system is in the thermodynamic ground state. Therefore, there still exists a possibility that the system is in the elementary excitation state. Let us recall that our purpose is to satisfy the relationship in Eq. (32) to establish the equivalence of the two different SPH forms. Particularly, satisfying Eq. (32) indicates that the motion equation of an inviscid flow and the nonlinear Schrödinger equation of a bosonic system are equivalent in their SPH forms when the density variation is moderate. Once this equivalence is validated, the motion equation of an inviscid flow can be coupled with the Navier–Stokes equation with theoretical validation, unlike in our previous studies, which experimentally adopted the coupling method. As long as the fluid speed does not exceed the critical speed for an elementary excitation, we can assert that the system is in a state close to the real ground state, i.e., a state that can be sufficiently regarded as the ground state from a thermodynamics perspective. This can be achieved when the speeds of all fluid particles do not exceed the critical speed because the total fluid speed is equal to or lower than the speed of an individual fluid particle in SPH. In this case, we can conclude that the system satisfies Eq. (32) for all fluid particles, and hence, the equivalence holds.

In addition, because *V*, *N*, and *S* are constant parameters, the total differential of Eq. (32) can be obtained as follows:

This indicates that the Gibbs–Duhem relationship holds for each fluid particle when Eq. (32) is satisfied. These findings can be summarized as follows: If the spatial variation in the condensate density is sufficiently small such that the quantum pressure is negligible and if the internal energy of each fluid particle remains zero, the motion equation of the superfluid component becomes equivalent to the motion equation of the condensates derived from the GP theory. Here, in the former, fluid particles are driven by the gradient of the chemical potential obtained using the Gibbs–Duhem equation, and in the latter, fluid particles are driven by the gradient of the chemical potential obtained from the Schrödinger equation of interacting bosons. The middle part of Fig. 3 presents a schematic summary of these results.

If the quantum pressure, $Pq(r,t)$, is non-negligible, the quantum pressure gradient force, $fq(ri,t)$, in Eq. (27) denotes the mutual friction force for the superfluid component of the two-fluid model in a counterflow. This implication is based on the following consideration. From a review of previous studies, we can conclude that the velocity distribution in counterflow experiments exhibits a flat profile along the moving direction as the amount of heat input to the system, *W _{in}*, increases.

^{71–73}This observation is different from the parabolic profile that follows the Hagen–Poiseuille equation of a laminar flow.

^{74,75}Specifically, the temperature gradient becomes proportional to

*W*cubed as

_{in}*W*increases. Gorter and Mellink introduced a pair of mutual frictional forces $Fsn(r,t)$ into the two components, as shown in the upper part of Fig. 1, to explain the discrepancy between the results of Landau's two-fluid model and experiments.

_{in}^{47}Hall and Vinen corroborated that the mutual friction force can be attributed to the interactions between the normal fluid component and quantum vortices.

^{66,76–78}In this manner, herein, the mutual friction forces are introduced as quantum mechanical corrections. In our formalism, $\u2212fsn(r,t)=\u2212(1/\rho s)Fsn(r,t)$ is added to the right-hand side of Eq. (28). Hence, the discretized form of the mutual friction force, $\u2212fsn(ri,t)$, for the

*i*th fluid particle is added to the right-hand side of Eq. (30).

Equations (c) and (d) in the lower part of Fig. 3 represent Eqs. (27) and (30), respectively, with the addition of $\u2212fsn(ri,t)$. The following important findings are obtained based on these two SPH forms. If we denote the chemical potential at *P _{q}* = 0 as $\mu \xaf$, based on Eq. (23), relation $\mu =\mu \xaf+Pq$ holds, indicating that $\mu \xaf$ is the classical thermodynamic chemical potential. In addition, from Eq. (26) and relation $\mu =\mu \xaf+Pq$, the first term in Eq. (c) represents the discretized form of the gradient of the classical chemical potential, $\mu \xaf$, which naturally corresponds to the first term in Eq. (d), owing to the Gibbs–Duhem relation. In contrast, $fq(ri,t)$ denotes the gradient of the quantum pressure scaled by $\u2212mq$, and $fsn(ri,t)$ denotes the quantum mechanical correction. In brief, the first terms in Eqs. (c) and (d) represent the discretized equation of the classical chemical potential, and the second terms in Eqs. (c) and (d) represent the discretized equation of a quantum mechanical correction given by either $fq(ri,t)$ or $fsn(ri,t)$. Therefore, it is reasonable to hypothesize that the equivalence of (c) and (d) is established when each term is equal to the corresponding term. The condition for the equivalence of the first terms in Eqs. (c) and (d) can be derived by comparing each term in the sum, which satisfies the relation $STj\u2212VPj+N\mu \xafj=0$, where $\mu \xafj$ is the classical chemical potential of the

*j*th fluid particle. From the Euler equation, we can observe that the equivalence of the first terms holds if the internal energy of each fluid particle is zero. Thus, Eqs. (c) and (d) are equivalent if their second terms are equal. These findings can be summarized as follows. If the quantum pressure is non-negligible and if the quantum pressure gradient force equals the mutual friction force, the equivalence of the two different SPH motion equations separately derived from a microscopic or macroscopic viewpoint is established. This is valid under the condition that the internal energy of each fluid particle is zero.

Thus, we determined the condition for the equivalence of the motion equations when the quantum pressure becomes non-negligible owing to the large spatial variation in the condensate density profile. However, this scenario rarely occurs for most incompressible flows based on SPH because the density fluctuation is typically below 1% of its averaged density in SPH calculations.^{79–81} Therefore, as long as Eq. (32) is satisfied for each fluid particle, the equivalence of the two different SPH forms is almost always guaranteed. Nevertheless, as an exceptional case, the quantum pressure may be non-negligible at the interface of two different types of fluid particles owing to the differences in their densities. However, the detailed dynamics at such interfaces are yet to be revealed. In the following, we present our perspectives on the equivalence of the two SPH forms at the interface of two different fluid particles.

In our opinion, the equivalence of the quantum pressure gradient and the mutual friction forces at the interface of two different fluid particles can be maintained, particularly in the case of weakly compressible flows included in explicit SPH simulations of incompressible flows. A thought experiment on a counterflow can explain the foregoing. As the heat input increases, the relative velocities of both fluid components increase, and the mutual friction forces must be included. Concurrently, in the SPH simulations with an explicit time-integrating scheme, the spatial change in the velocities of fluids generates a density variation in time according to the equation of continuity, owing to the weak incompressibility. This variation is minimal, below 1% for most parts of the flow, as mentioned above; however, it becomes evident at the interface owing to the density difference and yields a steeper density gradient with time evolution. Consequently, a quantum pressure gradient force can be induced with the increasing heat input, similar to the case of mutual friction forces. In this manner, the scenarios in which both forces coexist are similar, and we can expect their equivalence to hold to some extent. Interestingly, the weak compressibility in explicit SPH is often regarded as a drawback for classic incompressible fluids; however, it is a significant characteristic that can reproduce the effect of quantum mechanical corrections in quantum fluids. Specifically, the density variation in explicit SPH can create quantum many-particle-interacting systems; this could be one reason for the successful reproduction of the quantum mechanical phenomenon of vortex lattice formation in our previous study.

## IV. APPLICATION TO PRACTICAL PROBLEMS

We can assert that the objective of this study has already been accomplished based on the discussion in Sec. III and the following results. If the spatial variation in the condensate density is sufficiently small such that the quantum pressure can be negligible, the two discretized motion equations in their SPH forms, which are derived separately from the GP equation and the two-fluid model, are equivalent, provided that the internal energy of each fluid particle is zero. This equivalence holds even when the quantum pressure is non-negligible if the quantum pressure gradient force equals the mutual friction force. These results theoretically justify our method of coupling the motion equations of the two-fluid model, which was adopted experimentally in our previous studies. Coupling the motion equation of an inviscid flow and the Navier–Stokes equation is seemingly inadequate to capture the dynamics of quantum fluid phenomena because both these equations are classical equations of motion. However, at least in the SPH form, because the motion equation of an inviscid flow is equivalent to the GP equation under appropriate conditions, our coupling method is equivalent to the method that couples the GP equation with the Navier–Stokes equation in such cases. The latter method is seemingly more consistent with well-accepted methods that couple the quantum mechanical equations with the Navier–Stokes equation than the former, which couples the motion equation of an inviscid flow with the Navier–Stokes equation. However, our analysis indicates that they are equivalent, at least in their SPH forms, as long as the internal energy of each fluid particle remains zero when the quantum pressure is negligible. In summary, our method of coupling the motion equation of an inviscid flow with the Navier–Stokes equation is valid for simulating liquid helium-4 if the aforementioned condition is satisfied. Therefore, in this section, we present the use of this method in simulating the dynamics of liquid helium-4, as in our previous studies. As a developmental step in our study, we discuss the implications of these findings for solving real-world cases and their incorporation into liquid helium-4 simulations using our two-fluid model to accurately reproduce actual phenomena.

As previously mentioned, the condition requiring that the internal energy of a fluid particle be zero implies that the fluid particle is in the thermodynamic ground state. Because this is a thermodynamic condition, the ground state includes an elementary excitation state. One approach to sufficiently satisfy this condition is to ensure that the velocity of each fluid particle does not exceed the Landau critical velocity. Under the assumption of a weakly interacting Bose–Einstein condensate, it is generally known that the magnitude of the Landau critical velocity is given by the speed of sound.^{82} For helium-4, the critical velocity becomes lower than the speed of sound owing to the presence of rotons.^{83,84} Although this depends on the problem and measurement method, experiments involving highly sensitive flow measurement systems based on diaphragm displacement with a DC detection circuit have revealed that the critical velocity generally ranges from 7.7 to 11.6 $m\xb7s\u22121$. Moreover, it ranges from approximately 0.16 to 0.22 $m\xb7s\u22121$ in the slowest case observed during certain irregular transition states.^{85} Let us now focus on flow problems wherein the characteristic velocity of the flow is given by a few $cm\xb7s\u22121$, as a typical case of laboratory systems. Here, the fluid speed is always lower than the critical velocity when considering the simulations of liquid helium-4 at the laboratory scale. Therefore, the condition for equivalence can naturally be satisfied.

In addition, in quantum hydrodynamics, the circulation, *C*, of a vortex is given by an integer multiple of the quantum of circulation, *κ*, as expressed by $C=q\kappa $. Such a topological defect enables vortices to maintain dynamic stability. Here, *q* is typically one, and a situation wherein the circulation becomes more than twice *κ* is rarely encountered because the higher energy states are unstable. In the rotation problem, a vortex lattice is formed owing to the balance between the repulsive forces acting among multiple parallel vortices and the Magnus force. The latter acts as a centripetal force owing to the forced rotation of the outer vessel. Accordingly, in the simulations of rotating liquid helium-4, it is necessary to ensure that the magnitude of the circulation in each vortex is approximately equal to *κ* to accurately reproduce the effect of quantization of circulation. In our previous study,^{38} we monitored the circulation of clusters formed by accumulation of fluid particles at every time step. Therein, if the circulation of a cluster was greater than the quantum of circulation, *κ*, it was identified as a vortex. In this case, we obtained the repulsive forces based on the vortex dynamics theory^{38,51} between the vortex and other vortices or clusters to prevent further vortex growth. Consequently, we succeeded in reproducing the phenomenon of vortex lattice formation; however, we obtained a histogram wherein the circulation of each vortex was continuously distributed for circulation values greater than 1–3 times the *κ* value. However, we failed to quantize the circulation; therefore, we could not regard the formed vortex lattice as a “quantum lattice.” As a remedial measure, in this study, we apply a stronger constraint on the system; under this constraint, repulsive forces are applied to the interactions among clusters in advance to prevent them from merging when the sum of the circulations of two approaching clusters is estimated to be greater than *κ* to reproduce topological defects, and the constraint is expressed as follows:

where $Fint(k,l)$ denotes the resulting interaction force between the *k*th and *l*th clusters or vortices, and $fint(k,l)$ denotes the repulsive force of the vortex dynamics model. We use the same $fint(k,l)$ as in Sec. III D in Ref. 38. $\Gamma x\u2009(x=k,l)$ represents the circulation of the *x*th cluster or the vortex. *w _{c}* and

*γ*are model parameters that confine the circulation of a vortex at approximately the quantum of circulation,

_{c}*κ*, and they are set to approximately 0.9 and 1.1, respectively. Equation (34) can prevent vortices from accelerating and exceeding the critical velocity because of the unbalanced interactions among them, which result from the emergence of unrealistic vortices with circulations much greater than

*κ*.

We performed a numerical simulation of rotating liquid helium-4 using SPH under similar computational conditions as in Ref. 38, except for the conditional judgment represented in Eq. (34). Normal and superfluid particles were randomly distributed according to their density ratios of $\rho s/\rho $ and $\rho n/\rho $ in a cylindrical container with an outer diameter of 0.2 cm; this container began to rotate at a speed of $5\u2009rad\u2009s\u22121$ after the simulation was initiated. The temperature was maintained at 1.6 K throughout the simulations. Practical simulations using SPH require the introduction of several established techniques to ensure numerical stability. Additional details on these techniques can be found elsewhere Sec. II C in Ref. 38. For the boundary conditions, we first considered the concept of the conventional two-fluid model in the target system and then approximately reproduced the mechanics of our model. Specifically, we assumed a situation wherein the normal fluid component of the conventional two-fluid model is arranged alongside the outer vessel. In this scenario, the wall can interact with the normal fluid component via the viscosity force, and it can also interact with the superfluid component via the mutual friction force. Therefore, in general terms, the walls are in nonslip states against the fluid component. Based on this consideration, we applied a nonslip boundary condition on the system by arranging viscid particles alongside the outer vessel. When the rotation of the outer vessel began, we imposed the velocity of the vessel as a moving boundary condition to the wall particles. The results are presented in the upper row of Figs. 4(a)–4(c). These illustrate the two-component view, the spatial velocity distribution depicting the local velocity fields of the spinning vortices scattered inside a global constant-velocity field, and an enlarged view of multiple spinning vortices represented by rotational forces, respectively. Similar to Ref. 38, it is confirmed that multiple spinning vortices form a rigid-body lattice rotating under a constant-velocity field comprising the local velocity fields of spinning vortices that are indifferent to the motion of the overall velocity field. In more detail, similar to our previous study, we observed small clusters close to the periphery of the vortex lattice. These clusters are bundles of fluid particles that do not grow into vortices; we cannot regard them as vortices because they are just clusters that do not spin. The physical significance of these clusters before they develop into vortices is still under discussion. However, we also observed small clusters of normal fluid particles that were initially attached to the wall and began to detach from the wall owing to the centripetal force resulting from the forced rotation of the outer cylinder.

Furthermore, we obtained a histogram of circulations, which is presented in Fig. 4(d), where the horizontal axis indicates the circulation scaled by the quantum of circulation, *κ*, and the vertical axis indicates the number of vortices. This figure reveals that all circulations of the vortices are distributed within a range of $0.9\u2009\kappa $–$1.4\kappa $, which is slightly wider but almost similar to the experimental result.^{86} In this study, the circulations appear to be quantized at approximately *κ*. We note the existence of a small amount of measured data in the region where $C/\kappa $ is smaller than one in the experimental results, which is not reproduced in the present simulations. Notably, an accurate reproduction of the histogram of circulation depends on Eq. (34). Our previous study^{37} has already revealed that the simple SPH discretization of the equations of the two-fluid model with angular momentum conservation can reproduce the rigid-body rotation of multiple vortices spinning in the same direction, even without using a vortex dynamics model. In a subsequent study,^{38} we formulated the vortex dynamics in SPH, describing the Magnus force generated by the rotation of each vortex and the repulsive force resulting from the interaction between vortices spinning in the same direction. The problem was to identify a basis for the determination of conversion of a gathering of particles into a vortex. Because vortex circulation is generally quantized to approximately $q=\kappa $, we used the aforementioned vortex dynamics model when the circulation of a gathering of particles was $q=\kappa $. Consequently, we generated a vortex lattice in our previous study.^{38} However, the quantization of circulation failed because numerous vortices with larger circulations than those observed in the experiments were generated. In a continuous model like the SPH, even if the vortex dynamics forces are applied to clustering particles after the system recognizes them as a vortex, the clustering process proceeds before the applied repulsive force effectively acts on them. Based on this, Eq. (34) improves the conditional judgment in Ref. 38 to be more sensitive to changes in the circulations of vortices. We can state that the conditional judgment in Eq. (34) is still tentative and needs further sophistication, primarily from theoretical aspects. Nevertheless, it is still significant that we succeeded in quantizing the circulations of vortices based on our continuum mechanical approach. The experimental values are for 1.3 K. At this stage of our study, the aim is to qualitatively compare the behavior at similar low temperatures. The fact that the circulation is quantized at approximately $C/\kappa =1$ still remains valid, and we believe that this degree of temperature difference does not affect this discussion within the scope of this study. Nevertheless, a more quantitative comparison based on experiments at the same temperature is planned for the future.

We previously demonstrated that in our system, by adjusting the parameter controlling the magnitude of the angular velocity around the axis of each fluid particle (*C _{w}* in Ref. 38) the number of generated vortices could match the theoretical value estimated by Feynman's rule. However, the sizes of the vortex lattice and its vortices were different from the experimental values.

^{87,88}Our current view on this point can be detailed as follows. Let us recall that the balance between the Magnus force generated by the spinning of the vortices and the repulsive force between the vortices determines the vortex arrangement. As the angular velocity of fluid particles become small, the magnitudes of the Magnus force, repulsive force, and circulations of the fluid particles forming each vortex reduce because they depend on the angular velocity, and hence, the system forms a smaller vortex lattice. At this time, if the angular velocity of each fluid particle becomes lower than a certain critical angular velocity, the rotational energy of the vortex becomes insufficient, and the vortex is unable to spin; consequently, it is no longer a vortex. By contrast, if the fluid particles are large, i.e., if the resolution of the system is insufficient, the rotational energy of each vortex in the calculation is larger than that of a real vortex. Nevertheless, we can maintain sufficient rotational energy to observe vortex spinning by setting a higher angular velocity for each fluid particle. In summary, the scale of a vortex lattice can be controlled by the radius of a fluid particle, i.e., the resolution applied to the system, and the magnitude of the angular velocity of a fluid particle. This suggests that our simulation results display a vortex lattice magnified by the amplified angular velocities of the fluid particles. In all cases, it is necessary to ensure that the velocity of a fluid particle does not exceed the critical speed. To ensure the foregoing, the angular velocity considered is the angular velocity of classical fluid particles, with no quantum mechanical significance. Their interpretation in terms of quantum mechanics should be defined in future studies. Once such quantum mechanical explanations are given, we may move on to the next stage of the study, such as a detailed justification for using internal spin degrees of freedom of fluid particles and repulsive forces, explicitly formulating the local thermal equilibrium in our SPH model, and a detailed description of the internal dynamics of vortices.

In these numerical experiments, the density ratio of the two different virtual fluid particles causes the low-density normal fluid particles to aggregate and form vortices. As explained in Sec. II B, all fluid particles considered in our calculation are classical fluid particles that only act as discretization points of the helium-4 continuum. The two types of virtual particles (superfluid and inviscid) have no physical significance other than serving as discretization points, or fragments, of a continuous two-phase flow. Particularly, they are virtual particles. This study adopted a one-fluid system, wherein we focus on the fluid phenomena of the system as a whole; the motions of individual virtual particles have no physical significance from a microscopic viewpoint. It is possible to display these two types of virtual fluid particles separately; nevertheless, our focus is only on the entire fluid field. Therefore, a vortex cannot be decomposed into physically significant quantum components in our current model. Nevertheless, as Figs. 4(b) and 4(c) display the velocity and force of fluid particles, physical quantities that are defined even for fragments of continuous two-phase flow, this definition suffices to confirm the phenomena of vortex lattices. In future studies, we can possibly represent the internal structure of a vortex by introducing subscale particles inside it and by modeling their interactions with microscopic Lagrangian dynamics, such as molecular dynamics, or by combining with the DFT. In particular, the SPH formulation is well suited to be combined with other Lagrangian approaches for this multiscale physics.

Furthermore, we simulated two-dimensional counterflow problems. We focused on a simplified scenario wherein a periodic condition was applied to a rectangular domain. Specifically, we set the simulation domain of (*L _{x}*,

*L*) as (1, 0.5 cm) and set periodic boundary conditions along the x direction. The characteristic velocity, which is the maximum velocity of the fluid particles calculated by multiplying the Mach number with the speed of sound, was $2.4\u2009cm\xb7s\u22121$. We set the average temperature as $1.6\u2009\u2009K$ and the temperature difference, $\Delta t$, as $0.01$ or $0.1\u2009K$ in the x direction. Fluid particles were driven by a constant-temperature gradient calculated from the input temperature difference. The other input parameters were the same as those in the rotation problem. To satisfy Landau's criterion, we measured the velocity profile along the x direction when the velocities of the fluid particles were lower than the characteristic velocity. Figure 4(e) presents a snapshot of the simulation. As shown in Fig. 4(f), the velocity profile in the horizontal direction is confirmed to be parabolic or flat when $\Delta t=0.01$ or $\Delta t=0.1\u2009\u2009K$, respectively. The velocity profile appears as a center-flattened profile in the x direction when the amount of heat input to the system increases. This profile characteristic is qualitatively consistent with the results of experiments

_{y}^{72}and numerical simulations

^{25,89}of the normal fluid component often observed in a deformation state.

^{25}Meanwhile, a tail-flattened profile is observed in the transition regime from the laminar to turbulent flow.

^{90}Thus, our SPH model can be instrumental in exploring the detailed mechanism and conditions under which these different profiles are observed. In particular, reproducing a tail-flattened profile is one future task; however, improving the boundary conditions may be necessary in such cases. A straightforward explanation for the simulation results of the counterflow can be provided from the viewpoint of classical fluid dynamics as follows. In a fully classical approximation including the fluid forces of both components and their interactions, the normal fluid component acts as a drag force against the superfluid component. In this test case, when we increased the temperature difference from $\Delta t=0.01$ to $\Delta t=0.1\u2009\u2009K$, the drag force increased and prevented the entire fluid from growing to create a parabolic profile before the fluid reached the critical velocity. In summary, our numerical scheme based on SPH was demonstrated by two representative flow problems: rotation and counterflow problems of cryogenic liquid helium-4.

The resolutions of each simulation were approximately 200 000 particles for the rotation problem and 52 400 particles for the counterflow problem. All simulations were performed on a GPU NVIDIA Geforce RTX 2080 Ti system. However, it is likely that larger numerical simulations with higher resolutions may capture more precise phenomena. In particular, in this study, we did not examine the internal dynamics within each vortex. This was because the stability of vortices had to be discussed after the development of a theoretical scheme for the inner structure of vortices in our SPH model. Nevertheless, these detailed dynamics should be explored in numerical simulations with a sufficiently high resolution to capture the internal dynamics of a vortex. Another approach could be to combine our approach with other methods, e.g., the DFT, which is suited to describe the underlying physics in such a subscale regime, instead of simply applying SPH to the two-fluid model. As mentioned in Sec. II B, we can interpret fluid particles in SPH not necessarily as physically meaningful particles but as analytic points. Based on this, the differential operators in SPH can also be used in a DFT model. Nevertheless, it still appears promising to discretize the two-fluid model, conserving the angular momentum with ultra-scale resolutions using billions of fine-grained analytic particles to capture the dynamics of vortex nucleation or the inner structure of vortices. By comparing the vortex profiles from the high-resolution SPH with those from the DFT calculations, we can also expect to obtain the quantum physical meaning and optimal values of the model parameters of the SPH model, such as *C _{w}*, which controls the magnitude of the angular velocity around the axis of each fluid particle. We acknowledge that a few billion analytical particles can be used in SPH calculations within a reasonable computational time on recent supercomputers;

^{91,92}thus, we can state that SPH has a significant potential for diverse approaches. In addition, three-dimensional numerical simulations are required to explore the dynamics of vortex reconnection

^{93}or vortex rings,

^{94}or the dynamics of the rotating droplets in a spherical configuration,

^{95}and droplet beams.

^{96}To achieve this, it is essential to develop multi-GPU simulation codes that can efficiently realize computationally expensive large-scale simulations involving billions of particles; however, it is necessary to approach this problem from all angles in science and engineering.

## V. CONCLUSION

Apart from our own previous study, no other study has attempted to reproduce the dynamics associated with macroscopic quantum fluid phenomena using a fully classical mechanical approximation model. Our study was motivated by a purely academic interest aimed at directly reproducing macroscopic quantum phenomena, such as a film flow and the fountain effects. Our previous studies revealed that a vortex lattice in rotating liquid helium-4 can be reproduced even when the two-fluid model is solved in a fully classical mechanical approximation that includes the fluid forces of both components. This is valid under the condition that the viscosity is rederived to conserve the rotational angular momentum using SPH and the vortex dynamics is incorporated into the system. Furthermore, a fully classical mechanical approximation of the two-fluid model using SPH may be equivalent to solving a many-body quantum mechanical equation under specific conditions. However, no study has provided theoretical evidence supporting the existence of this equivalence of the microscopic motion equation of a quantum many-body system and the phenomenological motion equation of the superfluid component of the two-fluid model in an SPH formalism.

This study demonstrated the existence of this equivalence in the following manner. We first derived the SPH form of the motion equation of the superfluid component of the two-fluid model, i.e., the motion equation driven by the gradient of the chemical potential obtained using the Gibbs–Duhem equation. Furthermore, we derived the SPH form of the motion equation of condensates from the GP theory, i.e., the motion equation driven by the gradient of the chemical potential obtained from the Schrödinger equation of interacting bosons. We then compared these two discretized motion equations in their SPH forms, which were separately derived from a microscopic or macroscopic perspective, to identify a condition of their equivalence. Consequently, we discovered that the thermodynamic condition requiring the internal energy of each fluid particle to be zero ensures this equivalence when the quantum pressure is negligible. Our results also indicated that this equivalence holds even when the quantum pressure is non-negligible if the quantum pressure gradient force equals the mutual friction force.

These results theoretically justify our method of coupling the motion equations of the two-fluid model, which was adopted experimentally in our previous studies. Coupling the motion equation of an inviscid flow with the Navier–Stokes equation is seemingly inadequate to capture the dynamics of quantum fluid phenomena because both these equations are classical motion equations. However, at least in their SPH forms, because the motion equation of an inviscid flow is equivalent to the GP equation under appropriate conditions, our coupling method is equivalent to the method of coupling the GP equation with the Navier–Stokes equation in such cases. The latter method is seemingly more consistent with well-accepted methods that couple quantum mechanical equations to the Navier–Stokes equation than the former, which couples the motion equation of an inviscid flow with the Navier–Stokes equations. However, our analysis indicated that they are equivalent, at least in their SPH forms, as long as the internal energy of each fluid particle remained zero when the quantum pressure was negligible. In summary, our method of coupling the motion equation of an inviscid flow with the Navier–Stokes equation was demonstrated to be valid for simulating liquid helium-4 under a scenario wherein the aforementioned condition was satisfied. Therefore, it is still acceptable to adopt the method of coupling the motion equation of an inviscid flow with the Navier–Stokes equation in simulating the dynamics of liquid helium-4, as in our previous study.

In addition to the case of negligible quantum pressure, we also discussed the condition for the equivalence of the equations when the quantum pressure was non-negligible owing to the large spatial variation in the condensate density profile. However, this scenario is rarely encountered for most incompressible flows based on SPH because the density fluctuation is typically maintained below 1% compared to the averaged density in SPH. Thus, as long as the internal energy of each fluid particle is zero, the equivalence of the two differential SPH forms is almost always guaranteed. However, in exceptional cases, the quantum pressure may be non-negligible at the interface of two different types of fluid particles owing to the difference in their densities. However, the detailed dynamics at such an interface are still under discussion and require further study.

The condition requiring each fluid particle to have a zero internal energy can naturally be satisfied in several numerical simulations with a characteristic velocity of a few $cm\xb7s\u22121$ for a laboratory system. This condition indicates that each fluid particle is in the thermodynamic ground state. Because this is a thermodynamic condition, the ground state includes an elementary excitation state. One approach to sufficiently satisfy this condition is to ensure that the velocity of each fluid particle does not exceed the Landau critical velocity. For liquid helium-4, several experiments have indicated that the critical velocity generally ranges from 7.7 to 11.6 $m\xb7s\u22121$. Moreover, it ranges from approximately 0.16 to 0.22 $m\xb7s\u22121$ in the slowest case. Let us focus on flow problems wherein the characteristic velocity of the flow is given by a few $cm\xb7s\u22121$, as a typical case of laboratory systems. Here, when considering the simulations of liquid helium-4 at the laboratory scale, the fluid speed is always lower than the critical velocity. Therefore, we performed a simulation of rotating liquid helium-4 with a sophisticated constraint such that the velocities of the fluid particles did not exceed the Landau critical velocity. Consequently, we generated a vortex lattice with quantization of the circulation, known as a quantum lattice.

Our method provides a perspective on the dynamics of liquid helium-4 based on a fully classical mechanical approximation that is different from the descriptions obtained using conventional methods. Our theory and simulation results demonstrate that our proposed scheme can be used to appropriately describe the macroscopic dynamics of liquid helium-4. Particularly, a vortex lattice with quantized circulation, i.e., a quantum lattice, can be reproduced by this method based on classical fluid mechanics, although this has been previously considered as a purely quantum mechanical phenomenon. We expect that our approach provides a new methodology for describing cryogenic liquid helium-4.

## ACKNOWLEDGMENTS

This study was supported by JSPS KAKENHI Grant No. 22K14177. The author would like to thank Editage (www.editage.jp) for English language editing. The author in particular, acknowledges Professor Katsuhiro Nishinari and the administrative staff at the Nishinari Laboratory and the RCAST, University of Tokyo. The author is also grateful to his family for their warm encouragement.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Satori Tsuzuki:** Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## References

_{4}