A hydrodynamic approach to reproduce multiple spinning vortices in horizontally rotating three-dimensional liquid helium-4

Liquid helium-4 has a long history of use in various engineering and scientific fields. In particular, it has been used as a refrigerant for superconducting magnets in magnetic resonance imaging (MRI) and has been studied for cooling superconducting motors for rapid long-distance mobility. In general, two types of cooling methods are used in superconducting systems. The first is low-temperature superconductivity, an example of which is a niobium–titanium alloy cooled by liquid helium to near absolute zero (0 to about 4.2 K).^1,2 The second type is high-temperature superconductivity, in which a special metallic compound is cooled to a relatively high temperature of approximately 100 K using liquid nitrogen to achieve a superconducting state. Low-temperature superconductivity has the advantage of high critical current density. However, according to quantum mechanics, liquid helium loses its viscosity when cooled to near absolute zero owing to Bose–Einstein condensation (BEC) to the ground state of energy, and this quantum effect often causes a severe leakage phenomenon (super leakage). In addition, it only enters the solid state under high pressure and thus must be used as a coolant in the liquid state. In other words, the stable control of liquid helium in a cooling system is very difficult in low-temperature superconductivity. Stable cooling with liquid nitrogen enables high-temperature superconductivity. However, the critical current density is small, and the current that can maintain the superconducting state is weak; therefore, it cannot withstand practical use. Thus, neither cooling method is one step away from practical application owing to the trade-off between the difficulty of cooling and the critical current density; this is a factor that hinders the full industrialization of superconducting motors.

If it were possible to perform large-scale simulations on recent supercomputers to reproduce the dynamics of bulk liquid helium-4, the following would be expected. First, it would be possible to precisely predict and control the occurrence of super leakage due to quantum effects in the dynamics of superfluid helium; this will enable the stable control and safety of superfluid helium in large-scale cooling systems at a low development cost. Consequently, the practicality of superfluid helium as a cold agent in low-temperature superconductivity is expected to be dramatically enhanced. In addition, it is expected to significantly reduce the consumption of liquid helium during cryogenic experiments. It has been reported that the carbon dioxide emitted from the production and transportation of liquid helium can reach 587 kg-CO₂/l.³ If the dynamics of bulk liquid helium can be simulated by large-scale simulation technology using a supercomputer, the number of experiments using liquid helium can be reduced by half, leading to a reduction in carbon dioxide emissions.

This study developed a new simulation technique that can reproduce the three-dimensional (3D) fluid behavior of superfluid helium-4 in a cryogenic cooling system on a real scale. As of yet, most simulation methods for liquid helium have been based on quantum mechanical scale simulations that aim to elucidate the principles, and there has been little progress in the development of numerical methods for the so-called “bulk state” quantum fluids. Next, we discuss the properties of helium-4 at cryogenic temperatures. In principle, the mechanism of liquid helium can be described using quantum mechanics, which is the governing law of the nanoscale. Helium-4 is a so-called Bose particle consisting of two neutrons and two protons. Bose particles have quantized angular momentum times an integer. In this paper, “helium” indicates helium-4 unless otherwise noted. Liquid helium-4 can be described as a many-particle interaction system of Bose particles. However, unlike ordinary materials that solidify in low-temperature regions, helium-4 exists as a liquid even near absolute zero. Therefore, liquid helium has the strange property of locally obeying quantum mechanical laws while behaving as a fluid, a continuum, as a whole. In a many-particle interacting system of Bosonic particles below a critical temperature near absolute zero, helium atoms in a particular region undergo Bose–Einstein condensation (BEC), in which they fall into the ground energy state; this minimizes the interaction between helium atoms. (Helium is a monatomic molecule, but we use the term “molecular viscosity” in accordance with the chemical engineering doctrine.^4,5) In other words, liquid helium-4 loses its viscosity at cryogenic temperatures below the critical temperature, that is, it becomes an inviscid fluid.

The fact that liquid helium becomes an inviscid fluid causes peculiar fluid phenomena, such as slipping through the microstructure of the porous medium (gypsum), which is not observed when liquid helium is composed only of viscous normal flow components, and film flow caused by the viscosity between the wall surface and helium atoms being larger than the viscosity between helium atoms. These phenomena are often described as “macroscopic quantum phenomena” because the disappearance of microscopic interactions, i.e., of molecular viscosity, affects the macroscopic fluid behavior of bulk liquid helium. In addition to the phenomena that can be observed with the naked eye, vortex lattice phenomena observed on the nanometer scale in rotating liquid helium are also macroscopic quantum phenomena that can be observed between the hydrodynamic and quantum mechanical regimes. To the best of our knowledge, no study has successfully simulated the vortex lattice phenomena of rotating 3D liquid helium-4 observed in the boundary regime between these two different dynamical regimes using a method based on classical fluid mechanics.

It is widely recognized that the macroscopic behavior of liquid helium-4 can be phenomenologically described by the two-fluid model proposed by Landau, which consists of an inviscid fluid and incompressible Navier–Stokes equations. However, it should be noted that this assumption has only been verified by the issue of countercurrent flow, in which shear viscosity is dominant. Macroscopic quantum phenomena observed in the rotation problem, such as vortex lattices, have yet to be numerically reproduced using two-fluid models. The authors hypothesized that the main reason for this is that the ordinary two-fluid model does not reflect the quantization of the spin angular momentum of individual helium-4 atoms when liquid helium-4 is considered a microscopic many-particle interacting system. As stated previously, the spin angular momentum of the Bose particles was quantized. If fluid particles, which are virtual constituent particles in fluid mechanics, can be interpreted as a type of coarse-grained model of such microscopic particles, then the angular momentum related to the spin of the constituent atoms should also be conserved in the fluid particles.

Considering this hypothesis, the author recently developed a two-fluid model with spin-angular momentum conservation.^6–8 The following overview provides a general introduction to this subject. The Navier–Stokes equation with conservation of rotational angular momentum around the axes of fluid particles was originally derived from the constitutive law of continuum mechanics by Condiff and Dahler in 1964.⁹ Half a century later, Müller et al.¹⁰ developed a particle approximation of the Navier–Stokes equation with spin-angular momentum conservation using smoothed-particle hydrodynamics (SPH)^11,12 for mesoscale thermal flow. Considering these studies, the author adopted the Navier–Stokes equations with spin-angular momentum conservation for normal fluids and employed SPH to discretize the equation systems using particles. In other words, we describe liquid helium-4 as a system of numerous particle interactions composed of two distinct types of classical fluid particles that maintain their rotational angular velocity around their axes under the assumption of rigid body rotations; this is referred to as a two-fluid model with spin-angular momentum conservation. More importantly, our two-fluid model differs from conventional models in terms of the physical picture. The conventional two-fluid model was derived by combining the nonlinear Schrödinger equation for bosons with the thermodynamic Gibbs–Duhem equation. Conversely, our model considers the bulk of liquid helium-4 at low temperatures as an inviscid fluid following continuum mechanics, which contains residual viscosity at finite temperatures. In summary, the proposed model is based on a two-phase flow. Despite the fully classical mechanical picture, the resulting system equations were consistent with an ordinary model [see Eqs. (1) and (2)]. In particular, we consider the bulk liquid helium-4 to be an inviscid fluid, assuming that the viscous fluid component remains at finite temperatures. As the temperature decreased, the amount of the viscous fluid component decreased, ultimately becoming a fully inviscid fluid at absolute zero. From the Lagrangian perspective, our two-fluid model is a mixture of fluid particles from inviscid and viscous fluid components. Both fluid particles are classical fluid particles; they are virtual particles. Thus, our model is based on the mixing of two classical fluids and is similar to a two-phase flow model. However, as described below, we allow for a certain “nature of particles” of both components as a correction. In this context, the physical picture in which the viscous fluid component remains bead-like in an inviscid fluid is a more accurate representation of the system than a two-phase flow model. This is also close to the image of a “nanofluid,” or nanoparticle dispersion system, where the nanoparticles maintain fluidity. Details of this technique are described in Refs. 6–8. This will be outlined later in this paper. In our previous simulations of rotating two-dimensional (2D) liquid helium-4, we observed the formation of multiple spinning vortices that occurred independently of the overall fluid motion. After incorporating the vortex dynamics model into our model, we found that multiple spinning vortices formed a lattice with a regular arrangement. The number of vortices generated is determined by the strength of the rotational angular velocity. For appropriate values, the number of vortices agrees well with the theoretical solution of the number of vortices calculated from Feynman's law. It is also noteworthy that the spinning of each vortex moves independently without affecting the behavior of the entire fluid. In other words, the vortices do not dissipate. In a quantum vortex, the vortices do not dissipate because of the circulation quantization. This important property can also be reproduced. Consequently, the vortex lattice phenomenon in the rotating 2D liquid helium-4 was simulated by solving the two-fluid model with spin angular momentum conservation.

This study presents the observation of independently spinning three-dimensional vortices in 3D liquid helium-4 in a cylindrical container rotated horizontally using the aforementioned two-fluid model with spin angular momentum conservation. In the context of quantum hydrodynamics, the dynamics of the connections or reconnections of 3D vortex lines differ from those in classical dynamics owing to quantized circulations. Consequently, the topological variations in the 2D vortex lines are not analogous to those in the 3D vortex lines. For example, two vortex lines can occupy skew positions relative to one another in space, resulting in only two distinct patterns of reconnection between vortex line-line interactions and vortex line-wall interactions. In addition, from the perspective of polyhedral geometry, angle management in 3D space requires the use of quaternions to circumvent singularities. The dynamics of quantum fluids in 3D space are complex and, in fact, not an extension of those in 2D space; the cross section of a vortex line in real cases is projected onto a plate in 2D space. The geometric complexity in 3D space indicates that computational methods in 3D space must be verified and validated independently of the results in 2D cases. Therefore, it is unclear whether the authors' method, which was effective in solving the 2D rotating problem, can numerically reproduce the dynamics of 3D liquid helium-4. For clarity, we excluded the vortex dynamics model effective for 2D liquid helium-4 designed in our previous studies⁷ when applying our two-fluid model to 3D cases. We only focused on discretizing our two-fluid model in 3D space. Artificial manipulations such as those proposed by Schwarz^13,14 were not considered in this study. Nevertheless, we succeeded in reproducing several 3D vortices spinning in the horizontal direction during the horizontal rotation of liquid helium-4. The scientific significance of our numerical results is discussed in Sec. IV. The remainder of this paper is organized as follows. In Sec. II, we provide an overview of the proposed two-fluid model and the 3D calculations on a parallel computing platform. In Sec. III, we present the numerical results of our simulations of 3D liquid helium-4 rotating horizontally. In Sec. IV, we discuss the scientific significance of the results obtained. Finally, Sec. V concludes the study.

Our two-fluid model differs from the ordinary two-fluid model in its derivation and interpretation of the equation. Here, the ordinary two-fluid model is derived from the following microscopic relations. More specifically, the Navier–Stokes equation for an inviscid fluid with a temperature gradient term was derived by solving the Gross–Pitaevskii equation, which is a nonlinear Schrödinger equation, and the Gibbs–Duhem thermodynamic relation in conjunction. This equation describes the dynamics of the superfluid component. In short, the usual two-fluid model is obtained by making “macro-corrections” to the microscopic relational equations. The governing equation for the normal flow component was obtained by subtracting the equation for the superfluid component from the momentum balance equation of the total fluid.^15, Figure 1(II) schematically shows the relationship among the Gross–Pitaevskii equation, Gibbs–Duhem thermodynamic relation, and the ordinary two-fluid model. In a domain where quantum mechanics is dominant, the superfluid and normal flow components exist independently, and the two components together are considered to conserve mass and momentum (although entropy is said to be carried only by the normal flow component). This microscopic picture was proposed by Tisza¹⁶ and Landau¹⁷ and was subsequently verified by countercurrent flow experiments. Originally, the two components were considered completely independent; however, experiments have shown that when the heat input is large, the two components exert a mutual frictional force proportional to the velocity difference between them.¹⁸ 

By contrast, the concept and derivation of our two-fluid model are quite the opposite. Notably, we obtain a system of two equations of the same form as in the two-fluid model described above by making particle corrections to the fluid equations, which are macroscopic expressions of relations in the domain dominated by continuum mechanics. Figure 1(I) presents a schematic of our methodology. The details are as follows: First, bulk liquid helium at cryogenic temperatures is considered an inviscid fluid. Subsequently, we assume that the viscous fluid component remains at a finite temperature, as illustrated in Fig. 1(a). The amount of the viscous fluid component decreased as the temperature decreased; at absolute zero, it became a completely inviscid fluid. In other words, our two-fluid model assumes the mixing of the fluid components of inviscid and viscous fluids. Thus, the “fluid particles” in our model are virtual particles in classical fluid mechanics. Liquid helium, on the other hand, is usually considered an incompressible fluid. However, helium atoms are known to cause Bose–Einstein condensations (BECs) when cooled to the critical temperature. (Interestingly, although liquid helium is commonly described as a multiparticle system of Bose particles, only 13% of it actually causes BEC even at absolute zero.¹⁹) In any case, the condensate components were no longer fluid. These state changes to and from the condensates are described using quantum statistical mechanics. In other words, from a microscopic viewpoint, the composition ratio of the condensate component to the normal fluid component fluctuates around the average value. Therefore, it can be assumed that the fluid volume fluctuated around the mean value of the equilibrium state. From a macroscopic perspective, the bulk of liquid helium-4 was compressed by a few percent to approximately 13% in the low-temperature regime, and the volume changed around the average. Accordingly, this physical picture is compatible with the fluid computational model of the explicit SPH method. The explicit SPH method uses a finite-particle approximation of the incompressible Navier–Stokes equations. However, instead of implicitly solving the Poisson equation for pressure, as is typically done for incompressible fluids, this method solves the equation of state of a fluid, for example, the Tait equation of state.²⁰ This approach essentially reproduces density fluctuations.^21–23 An important point is that we assumed weak compressibility for liquid helium. The ability to solve the Navier–Stokes equations for an incompressible fluid while reproducing its weakly compressible properties is a feature of explicit SPH. In summary, in our two-fluid model, liquid helium is considered an inviscid fluid, and it is assumed that the viscous fluid component remains a function of temperature. In other words, it is a mixture of inviscid and viscous fluids. In addition, weak compressibility was observed. Weak compressibility was reproduced using the explicit SPH method.

In our two-fluid model, we introduce the properties of multiparticle systems as a microscopic correction as follows: First, we assume that the ratio of the number of virtual fluid particles in the inviscid and viscous components is equal to that of the quantum mechanical superfluid and normal fluid components calculated using the BEC theory at this temperature. In other words, the ratio of the number of virtual fluid particles in the inviscid and viscous components can be described as $Ns/N$ and $Nn/N$⁠, where N_s and N_n are the numbers of helium atoms in the superfluid and normal flow components, respectively, and N is the total number of atoms that satisfy $Ns+Nn=N$⁠. The respective density ratios were assumed to follow BEC theory. In other words, the densities of the inviscid and viscous fluids normalized by the total density ρ are $ρs/ρ$ and $ρn/ρ$⁠, respectively, and they satisfy the relations $ρs/ρ=Ns/N$ and $ρn/ρ=Nn/N$⁠. Thus, $ρs+ρn=ρ$⁠. In our model, the virtual fluid particles have the same volume V_p, regardless of their particle type. Therefore, the density and mass of the fluid particles are proportional and can be expressed as $ρs=ms/Vp$ or $ρn=mn/Vp$⁠, where m_s and m_n are the masses of the fluid particles in the inviscid and viscous fluid components, respectively. Under this physical interpretation, the total density of liquid helium remains at ρ because the total mass M of the fluid remains constant at any temperature and the total volume V is also constant. Thus, it is reasonable to solve the system of equations for an incompressible fluid (although the use of explicit SPH also reproduces weak compressibility, as described above). The key points of our model are as follows. Let the mass of the total fluid be M, the density ρ, and the volume V; $ρ=M/V$⁠. Let ρ_s be the density of the inviscid component and ρ_n be the density of the viscous component. Let M_s be the total mass of the inviscid component and M_n be the total mass of the viscous component of the entire fluid. The relationship $ρs=Ms/V$ or $ρn=Mn/V$ holds for the entire fluid. From the relation $ρ=ρs+ρn=(Ms+Mn)/V=M/V$⁠, $Ms+Mn=M$⁠. Accordingly, the mass M is conserved throughout the system, and its component ratio is determined by the temperature T. The volume of the fluid is also always V since the total density ρ is always constant. Even if we correspond the composition ratio of the fluid particles to that of the actual microparticles, the continuum nature of the fluid is not lost.

Although the two components may exhibit particle properties as described above, in principle, they are “fluid components,” and therefore, an interface between the two components is considered to exist under certain conditions. Even if the temperature of the fluid is constant, the density difference causes a temperature gradient between the two fluid particles near the interface, which acts as interfacial tension. In summary, an interfacial tension proportional to the temperature gradient is generated in the opposite direction between the two fluid components, as shown in Fig. 1(b). In summary, the inviscid component follows the Euler equation with a temperature gradient term, whereas the viscous component follows the incompressible Navier–Stokes equation with a temperature gradient term in the opposite direction [Figs. 1(a) and 1(b)]. At this stage, only the temperature gradient force contributes to the interfacial tension between the two components. Temperature is a macroscopic state quantity. In other words, only the macroscopic physical properties contribute to interfacial tension. However, as the temperature approaches absolute zero, more local interactions between the two components are expected to contribute to interfacial tension. In other words, local properties have emerged as multiparticle systems. Therefore, in addition to the macroscopic temperature gradient force, the local interaction force $Fsn$ between two fluid particles is explicitly described as part of the interfacial tension. In classical fluid mechanics, each phase is subjected to an action proportional to the local velocity difference between the different phases at the interface of a two-phase flow. In the Lagrangian picture, the action force between two different fluid particles is proportional to the local velocity difference between the two components. More specifically, as shown in Fig. 1(c), the ith fluid particle is subjected to a force proportional to the velocity difference between its own velocity and the average of the velocities of different types of fluid particles, and its reaction force is reflected in the different types of fluid particles. Because we use the explicit SPH method described above, the interface model of the two-layer flow of the explicit SPH method can be used as is. We examined the conditions under which the interface between the two phases should be considered. First, the “interface” between the two phases is a fundamental concept in hydrodynamics. Therefore, the interface effect becomes non-negligible when the two components exhibit hydrodynamic rather than many-body properties. Liquid helium-4 is expected to behave as a uniform inviscid fluid when no external force is applied. However, the assumption of a single-phase flow of an inviscid fluid is broken by a strong external shock or heat, and an interface effect between the two components appears. As this situation simultaneously reveals the multi-particle nature of the fluid, the interfacial effect operates only in the direction of the applied impact force, and the force does not propagate in other directions as in a continuum. In the horizontal rotation problem, liquid helium-4 was forced to rotate horizontally; therefore, fluidity with weak compressibility may manifest strongly in the horizontal direction. Accordingly, the driving force due to the temperature gradient force was assumed to act in the horizontal direction. Interface effects in the vertical direction were not considered in our model because no external forces acted in this direction.

In this manner, we adopted the following assumptions. (i) It is a mixture of a viscous and an inviscid fluid at a finite temperature. (ii) The ratio of helium atoms in the ground state to the total number of atoms follows quantum statistical mechanics for bosons and thus fluctuates around its statistical average, enabling us to assume that the volume of the fluid fluctuates; in other words, weak compressibility can be assumed. (iii) Local interaction forces are exerted between two different types of fluid particles. Moreover, in another recent two-fluid model, the vortex filament model^24–27 was solved instead of the inviscid equation, in which the rotation of the vortices can be considered. In contrast, (iv) we assume spin angular momentum conservation of fluid particles around their axes; thus, the viscosity term is decomposed into three terms in Eq. (2). The aforementioned (i) corresponds to Figs. 1(a) and 1(b), (ii) corresponds to Fig. 1(e), (iii) corresponds to Fig. 1(c), and (iv) corresponds to Fig. 1(d). It should be emphasized that this paper is the first to clearly describe that cryogenic liquid helium-4 is a system in which viscous fluid components “remain” in the inviscid fluid at finite temperatures and that the residual components decrease with decreasing temperature, becoming a “uniform inviscid fluid” at absolute zero.

Let us summarize the physical significance of each parameter in Eqs. (1) and (2) with reference to the literature.^7, ρ_n and ρ_s are the mass densities of the normal and inviscid fluid components, respectively, which satisfy the relationship $ρ=ρn+ρs$⁠, where ρ is the mass density of the total fluid. $D{·}/Dt$ denotes the material derivative. $vn$ and $vs$ are the velocities of normal and inviscid fluid components, respectively. P, T, and σ denote the pressure, temperature, and entropy density, respectively. $η¯, ξ¯$⁠, and $ηr¯$ indicate the shear, bulk, and rotational viscosities, respectively, where $ηr¯$ determines the impact of the rotational forces.¹⁰ The rotational force here indicates the rotational force around the axis of each fluid particle. The vector $ω0$ in the fifth term on the right-hand side of Eq. (2) represents the angular velocity of each fluid particle along its axis. While the vorticity is defined on a continuum, the vector $ω0$ is a variable of the microscopic constituent particles, which in its original formulation represents molecular particles and later approximates fluid particles. Hence, the vorticity and $ω0$ are different variables. This is explained further. Originally, Condiff and Dahler⁹ reformulated the Navier–Stokes equations from Cauchy's equation of motion and the conservation law of angular momentum to include the internal freedom of molecular particles to account for the effect of their spins on the entire polar fluid, resulting in a Navier–Stokes equation with spin-angular momentum conservation. Because molecules are polarized in polar fluids, the influence of the rotational motion of individual molecules on the bulk fluid needs to be considered. A schematic of the derivation of this equation is presented in the  Appendix (Fig. 9). It is worth noting that Condiff and Dahler divided the total angular momentum (M) per unit mass of the fluid into the angular momentum of the bulk fluid (⁠ $r×u$⁠) and the contribution due to the internal degrees of freedom of the molecules, that is, the rotational angular momentum (I) [see Fig. 9(d) in the  Appendix]. The vectors r and u represent the coordinates and velocity of the local fluid fragment, respectively. I can be further decomposed into $I·ω0$⁠, where $I$ is a tensor field, which is a scalar multiple of the unit dyadic if we postulate a uniform and isotropic rotational field (and thus a rigid-body rotation) for each constituent particle. In this context, $ω0$ denotes the angular velocity of a molecule's rotation around its axis. Thus, $ω0$ is in principle a quantity defined for a molecule. However, $ω0$ is often redefined for each fluid particle comprising a collection of molecules, assuming that each fluid particle is a representative point among the subordinate molecules based on a coarse-grained approximation. Therefore, the term “fluid particle” here is closer to the definition in the context of molecular fluid dynamics and thus may be not the same as a fluid particle as defined in classical fluid mechanics. When considering mesoscale flow problems, such as molecular fluids, fluid particles are often interpreted in both the classical and microscopic senses; we still refer to them as classical to compare the quantum mechanical description of the ordinary two-fluid model in this paper.

In the Lagrangian picture, the difference between the vorticity $∇×u$ owing to the rotational motion of the bulk fluid and the vorticity owing to the rotation of the molecules can be assumed to cause strain. Specifically, we assume that the deviation of the vorticity from the rigid-body rotation yields strain because the molecules rotate rigidly. When the vorticity $∇×u$ at a point is a rigid-body rotation, the local rotation of the fluid is synchronized with the spin of the molecule; thus, $∇×u−2ω0=0$⁠. In this case, no viscous stress owing to local rotation occurs. Otherwise, the difference in $∇×u−2ω0$ induces strain in the fluid, making it deviate from rigid-body rotation, and contributes to the asymmetric part of the stress tensor [see Figs. 9(F) and 9(I) in the  Appendix]. In the third through fifth terms on the right-hand side of Eq. (2), the viscosity term is divided into three contributions: shear, volume, and rotational viscosities. Thus, the spin-angular momentum-conserving Navier–Stokes equations can be rephrased as Navier–Stokes equations with reformulated viscosity terms. This model can easily connect microscopic and fluid mechanics in a continuum mechanical regime. To date, it has been applied to mesoscale flows in which the rotation of molecules or fluid particles dominates, such as in polar fluids and suspension flows. Previously, Müller adopted SPH to discretize the spin-angular momentum-conserving Navier–Stokes equations to simulate two immiscible flow;¹⁰ their model was further applied to red blood cell flow.²⁸ 

Mathematically, SPH is a finite particle approximation of a continuum field that can decompose the quantities defined in the field into neighboring particle contributions. The intrinsic part of SPH is a representation of a physical quantity in an integral form using the Dirac delta function δ and its approximation by the finite distribution function W, called the smoothed kernel function. We then discretize the approximation form based on the concept of the sum approximation for numerical simulations. This paper does not focus on the derivation of the SPH models because they have been discussed in our previous studies;^6–8 however, the basic concepts of SPH and its discretized models used in this study are illustrated in Fig. 10 in the  Appendix. In Refs. 6–8, the viscosity term of the two-fluid model for cryogenic liquid helium-4 was reformulated to conserve the spin-angular momentum. In other words, the viscosity term was reformulated so that the two-fluid model can be applied not only to shear viscosity-dominated problems such as countercurrent flow but also to rotation-dominated problems. We then discretized the equations using the SPH. Consequently, we observed the formation of multiple spinning vortices that occurred independently of the overall fluid motion in our SPH simulations of the rotating 2D liquid helium-4.⁶ After incorporating the vortex dynamics model into our model, we found that multiple spinning vortices formed a lattice with a regular arrangement.^7,8 The number of vortices generated was determined based on the strength of the rotational angular velocity. For appropriate values, the number of vortices agreed well with the theoretical solution of the number of vortices calculated using Feynman's law.⁷ Compared with our previous studies, this study is the first to address a 3D problem. The system does not assume gravity in any direction, and considering numerical stability, the validation in the 2D cases is effective for the 3D cases as well. As we focus on the disk-shaped problem and consider only the horizontal rotation around the vertical axis (z) of the disk, the developed rotational viscosity term for the two-dimensional model can be used directly. However, because of the topological differences in vortices between the 2D and 3D cases, it needs to be studied whether the previous method, which is effective in solving the 2D rotating problem, can numerically reproduce the vortices in 3D liquid helium-4; large-scale simulations for 3D cases will be reported and discussed later.

The fourth term on the right-hand side of Eq. (2) becomes 0 owing to the incompressibility condition $∇·v=0$⁠. In this respect, the current model always imposes the incompressibility condition $∇·v=0$⁠, and the simulation is slightly weakly compressible by solving using the explicit SPH method. Therefore, at present, the fourth term is meaningless; however, in the future, it may be solved directly under a compressibility condition. Therefore, we describe the system of equations in a generic form as in Eq. (2). The viscosity term in Eq. (2) is an extended version of the usual expression: the set of the third to fifth terms converges to $η¯∇2vn$ because $ηr¯$ converges to 0 when $∇·v=0$⁠. This implies that $η¯$ corresponds to the shear viscosity η_n of the normal fluid. Meanwhile, the local interaction force $Fsn$ is expressed as follows: As mentioned earlier, $Fsn$ is proportional to the local velocity difference between the two components, and can be expressed as $Fsn=Clvsn$⁠, where C_l is a coefficient. This form exactly corresponds to the mutual friction forces between the two components in quantum mechanics; in this paper, we estimate C_l assuming that $Fsn$ corresponds to the mutual friction forces $Fsnq$ of the ordinary two-fluid model, where $Fsnq$ can be expressed as $Fsnq=2/3ρsακLvsn$⁠, where $vsn$ is the relative velocity $vs−vn$⁠, κ is the quantum of circulation, α is the friction coefficient, and L is the time-dependent vortex line density. Unlike in our previous studies, we use the initial value of L before its decay, $L−1(0)$⁠, where L satisfies Vinen's equation: $L−1(t)=L−1(0)+βvt,$ where β_v is a coefficient. For further details, refer to.²⁹ 

Let us explain the physical significance of each parameter in Eqs. (3) and (4) with reference to Refs. 6 and 7. In Eq. (3), M represents the mass of a helium-4 atom, N indicates the total number of helium-4 atoms in the system, π represents the ratio of the circumference of a circle to its diameter, and c is the speed of sound. $ℏ$ indicates the reduced Planck's constant and k_B represents the Boltzmann constant. The constant values μ, p₀, and Δ are given by $1.72×10−24 g, 2.1×10−19 gcms−1$⁠, and $8.9 K$⁠, respectively.^31,32 Equation (3) is always true within the scope of this study because it is valid if $T≪93 K$⁠.³¹ In the elementary excitation model, the first term in the parentheses on the right side of Eq. (3) represents the “photon” contribution, and the second term in the same parentheses represents the “roton” contribution; after finding the total pressure from the summation of the photon and roton contributions, the entropy density in Eq. (3) is obtained by the thermodynamic relation, such that the partial derivative of pressure with respect to temperature.⁶ Equation (3) represents the relationship between entropy density σ and temperature T in a quantum “ideal gas.” To obtain the value of liquid helium-4, we decompose σ as $σ=Ceσ0$⁠, where C_e indicates the volume ratio of liquid to gaseous helium-4 under an atmospheric pressure of 1.013 25 bar, which is $1.428×10−3$⁠.³³ Under constant temperature, we calculated σ₀ using Eq. (3) and maintained a constant value for entropy conservation. Meanwhile, in Eq. (4), ρ₀ and p₀ represent the fluid density and pressure in the initial state, respectively. α_P is the specific heat ratio that affects system incompressibility. β_P denotes the reference pressure, which typically corresponds to p₀. Equations (1)–(4) represent the governing equations of the system.

Notably, for the dynamics of cryogenic liquid helium-4, both the derivation of the governing equations from the mechanical interpretation in Figs. 1(a)–1(e) and quantum mechanics resulted in similar mathematical forms. The first model assumes a two-phase flow model in which the two components are separated and the fluid particles have a volume exclusion effect on each other. Particle correction has been introduced as an exception. Conversely, the second or ordinary two-fluid model assumes that both components coexist in space and that neither type of particle exerts volume exclusion effects on each other, and introduces mutual friction as an exception. In summary, the governing equation can be identical regardless of whether the system comprises separate or coexisting (overlapping) particles. However, separate and overlapping constituent particles lead to different physical meanings, and the former and latter models provide the classical and quantum-mechanical descriptions of the same problem. In particular, the physical meaning of “particles” is different; the former indicates classical fluid particles, and the latter represents quantum minute particles. This inspired us to reconceptualize the system in terms of the degree of volume exclusion between the two components, hereafter, referred to as the mutual volume exclusion (MVE) effect. Here, the MVE effect can be rephrased as the effect of classical fluid forces, such as surface forces, body forces, or interface tensions. This eliminates the need to distinguish the physical meanings of the particles between the two models, considering that the MVE is the primary factor causing the particles to separate or overlap.

Ultimately, the two models can be unified into one model in terms of the MVE effect as a variable. This model has the following advantages. New insights into the connection between classical and quantum mechanics can be obtained by using the MVE effect as a unified metric of the system. This may also provide a new perspective on the relationship between the classical two-phase flow and the ordinary two-fluid model of unified hydrodynamics. In engineering and industry, it may become possible to capture classical and quantum fluid phenomena with only one model by adjusting the degree of MVE of the constituent particles depending on the location and condition. Practically, there is no need to switch computational codes depending on the conditions, and a single computational code can be used to analyze the entire large-scale system.

In such a unified model, the transition from classical to quantum particles can, in principle, be automatically performed using high-resolution computational calculations. This is because as the computational resolution increases, the size of the fluid particles in the classical domain gradually decreases, and the excluded volume effect of individual fluid particles also decreases. In other words, a low computational resolution in a local region implies that the size of the analytic particles at that location is large, and a large analytic particle size implies that the excluded volume effect of the individual particles is large, which is suitable for describing the location where the classical fluid is dominant. Conversely, a high computational resolution in a local region implies that the size of the analytical particles at that location is small, which implies that the excluded volume effect of individual particles is small, which is appropriate for a local region where quantum hydrodynamics dominate. In summary, in the unified model, the convergence of the computational resolution is expected to correspond to the convergence of the physical phenomena. In this case, the introduction of multiresolution particles makes it possible to reproduce the behavior of the fluid in the entire system in a single model.

To this end, the immediate challenge of our approach is to provide evidence that increasing the resolution brings the simulation results closer to the phenomena observed in the quantum-mechanical regime, even when using our classical mechanical two-fluid model with spin-angular momentum conservation.

In high-performance computing, an effective parallel computing algorithm for the simulation of many-particle interacting systems should be chosen based on the characteristics of the target problem, computational scheme, and properties of the parallel computing environment. In general, these problems include gravitational many-body problems,^34,35 molecular dynamics (MD) problems,^36–38 Lagrangian fluid analysis using SPH,^39,40 deformation analysis of structures using the mesh-free method,^41–45 and powder problems using the distinct element method (DEM).^46–49 The interaction domains and connectivity between particles depend on the target problem. From a computational perspective, gravitational many-body problems and MD problems result in N-body calculations owing to the large interaction range. Consequently, the computational cost of the interaction calculations is $O(N2)$⁠. Conversely, DEM calculations for powders require interaction calculations of the repulsion and friction between the contacting particles. Therefore, the memory reference cost is more significant than the computational cost when appropriate techniques such as neighbor particle listing and particle data sorting^50–52 are implemented to reduce the interaction calculation cost to O(N). In contrast, fluid calculations using particle methods, such as SPH, stand at a middle ground between the N-body and DEM calculations in terms of both memory reference and computational cost, even if we use a technique similar to DEM. In a 3D calculation, each particle interacts with approximately a few hundred particles within the interaction domain, which is called the kernel radius. Moreover, the numerical accuracy of the spatial discretization in SPH is typically equivalent to that of the central difference method in the finite difference method (FDM). In addition, the time discretization accuracy was between the first- and second-order accuracies with respect to the resolution, depending on the time-integration schemes. Accordingly, it is imperative to capture the underlying real physics with the highest possible accuracy using as many analytical particles as possible. This is one of the reasons why the field of large-scale simulation of particle methods frequently appears in the ACM Gordon Bell Prize.⁵³ 

In recent years, graphics processing units (GPUs) have been installed in world-class supercomputers as accelerators to improve the computational performance, memory utilization, power efficiency, and other metrics. Multi-GPU computing can speed up simulations and ensure that sufficient memory is available. The development of applications that operate efficiently on GPU supercomputers requires CUDA programming.⁵⁴ This also requires making good use of hierarchical memory structures, such as the limited amount of device memory on the GPU board and high-speed shared memory on the chip. The GPU supercomputer consists of many nodes with distributed memory connected by high-speed networks. GPUs are installed on each compute node as computational accelerators. For particle simulations on a memory-distributed system such as a supercomputer, it is efficient to divide the entire computational domain in space and assign a node with accelerators to each divided subdomains. Among particle problems, calculations performed in gravitational many-body and MD are based on long-range interactions. Therefore, the computational cost is dominated by floating-point operations rather than memory accesses, which allows for high execution performance, even on GPU supercomputers with low byte/flops. However, in SPH simulations that compute short-range interactions, most computational time is spent on memory accesses, making it impossible for a low byte/flop supercomputer to achieve sufficient performance.

The computational cost per GPU for particle simulations using SPH depends on the number of particles in the subdomains to which each GPU is assigned. In a parallel computing platform such as a supercomputer, the distribution of particles among nodes changes in both time and space, resulting in a non-uniform distribution of particles per node. At this point, the computational load in regions that are densely populated with particles increases significantly, making large-scale computations impractical. Therefore, it is essential to introduce dynamic load balancing between nodes to achieve sufficient execution performance in simulations of short-range interaction-based particle methods, such as SPH. This involves repartitioning the computational domain according to particle distribution to ensure that the computational load on each node is uniform. However, achieving dynamic load balancing in simulations of short-range interaction-based particle methods on GPU supercomputers is challenging because of hierarchical memory structures. The implementation is further complicated by the need for specialized technologies specific to large-scale parallel computing, such as CUDA programming and the MPI library,⁵⁵ for the development of the computational code. In addition, when the particle distribution undergoes significant changes and moves in a random direction across nodes, it is essential to efficiently manage the particle data on the GPUs to minimize the communication cost between GPUs across nodes via host CPUs and between GPUs and CPUs within each node. It is critical to consider efficient computation methods for individual GPUs considering their different memory structures, such as the device memory on the GPU and the high-speed shared memory on the chip.

It is imperative to develop dynamic load balancing among the GPUs to exploit the high computational efficiency of multiple GPUs. Various load balancing techniques have been proposed in related fields, including the slice-grid method, which moves the domain boundary to a spatially partitioned region; graph theory-based partitioning (such as k-means), as represented by ParMETIS^56,57 and Zoltan;^58,59 and the hierarchically structured tree method, which uses a recursive structure of cells to partition a domain, as exemplified by kd-tree decomposition,⁶⁰ orthogonal recursive bisection (ORB),^61–63 and domain decomposition using a hierarchically structured grid with space-filling curves.⁶⁴ Recently, we conducted a study on optimal dynamic domain decomposition for particle simulations of short-range interactions compatible with low byte/flops supercomputers, such as GPU supercomputers. Subsequently, we developed a multi-GPU computing framework that enables large-scale particle simulations using SPH. The SPH computational framework⁵² developed by the author was used in this study. Based on the nature of the target problem, the 2D slice grid method was used as the domain decomposition method. Figure 2 illustrates a snapshot of the domain decomposition in the case where the 19.6 × 10⁶ fluid particles were divided into 32 subdomains, each of which was assigned to a single GPU to speed up the simulations.

Unlike the free surface problem, the particle distribution in the rotating liquid helium-4 simulation was not excessively biased. Therefore, to prevent the computational load from being excessively biased, we divided the domain at regular intervals, such as every several hundred steps, so that the overhead of domain re-segmentation was negligible, and the load was balanced. A workflow diagram of the calculations is presented in Algorithm 1. The SPH calculation process for the 3D simulation of liquid helium-4 using our two-fluid model was, in principle, the same as that for the 2D case reported in our previous study. The velocity Verlet method^65,66 was used as the time integration method, and a nonslip condition was imposed on the lateral wall boundaries. As this is a three-dimensional calculation, a slip condition is imposed on the top and bottom surfaces, assuming that the inviscid fluid continues. After the particle density was calculated, the particle pressure was determined using the Tait's equation of state. After applying the pressure filter,⁶⁷ the updated particle density and particle pressure were used to calculate the pressure and temperature gradient forces. The viscosity term was separately calculated for the rotational and shear viscosities of the particles. The collision model⁶⁸ computed in line 17 is a correction model that preserves the volume effect of the fluid particles. The local interaction forces between the viscous and inviscid fluid particles were then computed alternately, the boundary conditions were calculated, and the physical quantities of the particles were updated according to the velocity Verlet algorithm. For a detailed description of the computational model using SPH, refer to Refs. 6 and 7.

Given the significant increase in the computational cost associated with 3D calculations, it is imperative to run simulations on a multi-GPU platform. As previously mentioned, dividing the computational domain into subdomains and assigning a GPU to each subdomain is a highly efficient approach. To ensure the connectivity of the domains between subdomains, it is necessary to copy the updated particle data to the neighboring subdomains each time the particles are in the halo region, which is between the boundaries with the neighboring domain and the interior of the domain by the interaction length (in our implementation, twice the kernel radius), update their properties. As shown in Algorithm 1, five communications of particle data within the halo domain were required. In addition, the particle positions are updated after all computations in a single step were completed. Each GPU then transfers the out-of-domain particles that have left the subdomain to the target neighboring subdomain. Consequently, communication between neighboring GPUs was required six times. Frequent repetition of communication fragments the particle array and reduces the performance; therefore, periodic reordering of particle data is performed. This process enables efficient multi-GPU computation using SPH.

Let us set the outer diameter of a cylindrical container in the horizontal direction to 0.2 cm with a thickness of 0.04 cm in the vertical direction (z-direction). The resolutions of $(Nx,Ny,Nz)$ are given by (500, 500, 100), where N_a $(a=x,y,z)$ indicates the number of particles in the direction in a rectangular domain that covers the cylindrical container. We set the rotational angular velocity around the cylinder axis (z-axis) to 5 rad s⁻¹ counterclockwise. Similar to the two-dimensional cases in our previous work, we set the temperature T to 1.6 K and determined the values of $(ρs/ρ,ρn/ρ,κ)$ with those at T = 1.6 K listed in Ref. 69. In addition, η_n was obtained from the relationship of 0.566  $κρn$ with reference to Ref. 69. The shear viscosity $η¯$ and rotational viscosity $ηr¯$ have the value of η_n. We set the parameter $Cω$⁠, which determines the magnitude of the angular velocity $ω0$ and thereby determines the resulting number of vortices, to be 0.016 to maintain a similar degree of magnitude of the rotational forces as those in 2D cases reported in our previous studies. Normal fluid particles were placed along the inner circumference of the container as walls, and no-slip conditions were imposed on the wall particles as boundary conditions. Conversely, we set slip conditions at the top and bottom of the vessel. We applied a well-established wall and dummy particle method;^70,71 we update the density and pressure of the inner wall particles similarly to fluid particles, and retain the initial pressure and density values of the outer wall particles. Figure 3 illustrates the geometrical setup of the target system. In the initial state, superfluid or normal fluid particles were randomly generated in proportion to the density ratio and distributed in the vessel.

The computational conditions in this study were adopted from the parameter values determined in our previous study.⁶ In addition to satisfying the Courant–Friedrichs–Lewy (CFL) condition, several benchmark problems were solved to validate the stability and numerical accuracy of the calculations. Specifically, we simulated the propagation of the wave equation using the proposed SPH method and compared its numerical solution with the numerical and theoretical solutions obtained using the difference method. In addition, we compared the simulation results with the exact solution for the case in which the Reynolds number R of the normal flow component was 100, and for the case in which R was 1000 in 2D driven flow problems. We obtained very good agreement between the simulation results and the exact solutions obtained by the high-order FDM.⁷² Furthermore, assuming an ideal case in which the two components behave as perfect classical fluids, we solved the Rayleigh–Taylor instability problem under conditions where the density difference $(ρs/ρn)$ was five times, which is the same as that for the problem covered in this study, and confirmed that the characteristic mushroom structure could be reproduced. These results were confirmed for 2D problems. Although we targeted a 3D problem, we did not need to impose gravity on the system. In addition, the target problem was quasi-2D rather than symmetric 3D, and dynamic changes in the vertical direction were not dominant. Accordingly, we used the values used in previous simulations. Note that we set δt to be sufficiently small to satisfy the CFL condition for all resolution cases. We set δt to $5.0×10−6$ s and computed 180 000 time steps, which is 0.9 s in physical time. In the simulation, under the aforementioned conditions, the total number of normal and inviscid fluid particles, including the wall particles, was 19 636 400. We boosted our simulation using 32 GPUs installed on the Wisteria/BDEC-01 supercomputer at The University of Tokyo, Japan. We applied a two-dimensional slice-grid technique to the dynamic domain decomposition in the horizontal direction and decomposed the computational domain into 32 subdomains (four in the horizontal direction and eight in the depth direction), and a single GPU was assigned to each subdomain. Consequently, approximately 5.5 days were required to compute all 180 000 time steps of the simulation. For comparison, we performed simulations with a lower resolution using $(Nx,Ny,Nz)$ = (250, 250, 50) to investigate the dependence of the observed phenomena on the computational resolution.

Figure 4 shows snapshots of the rotation simulation with a lower resolution using 2 454 000 particles from the beginning to 0.79 s in physical time. The fluid particles of the inviscid and viscous components were visualized in color according to the strength of the rotational forces calculated by the fifth term of Eq. (2) using the volume rendering technique. The wall particles were removed for visualization. A video of the simulation is presented in Fig. 4 as integral multimedia in Fig. 5. A preprint version of the video is available at https://www.satoritsuzuki.org/video-horizontally-rotated-3d4he-2p4m. A single forced vortex that spans the container is generated in the rotation problem of a single-phase flow in classical fluid mechanics. However, it can be confirmed that the behavior of the fluid is similar to that of a quantum fluid. In other words, we observed that the low-density and viscous components coalesced to form multiple spinning vortices. By contrast, in this simulation, the vortices repeat arbitrary coupling and separation while continuously changing. This is characteristic of classical fluids. It is reasonable that the dynamics of the vortex should be similar to those of classical flow because it does not currently contain artificial rules for vortex recombination, as observed in quantum fluids. In the future, if, for example, Schwarz's rule^13,14 is introduced, it will be possible to get closer to the second image of the quantum vortex. Snapshots from the beginning of the simulation to 0.79 s in physical time in the high-resolution case with approximately 19.6 × 10⁶ particles are shown in Fig. 6. A video of the simulation is presented in Fig. 6 is provided as integral multimedia in Fig. 7. A preprint version of the video is available at https://www.satoritsuzuki.org/video-horizontally-rotated-3d4he-19m. Low-density and viscous components agglomerate to form multiple spinning vortices, which is common in the low-resolution case. This is also the same as in the low-resolution case in that the vortices experience arbitrary connection and separation with continuous change. However, unlike in the low-resolution case, the individual vortices are thinner and the number of vortices is larger, even though the computational conditions are the same.

A description of the dynamics of liquid helium-4 at cryogenic temperatures based on classical fluid dynamics yields Eqs. (1) and (2), which are similar to the conventional two-fluid model derived by combining the nonlinear Schrödinger equation for bosons, that is, the Gross–Pitaevskii equation, with the thermodynamic relation Gibbs–Duhem equation. In addition, our two-fluid model includes a conservation term for angular momentum related to the rotation of classical fluid particles. The fact that we could reproduce multiple spinning vortices, such as quantum vortices, using a method based on classical fluid mechanics would be groundbreaking because multiple spinning vortices and their lattice phenomena in liquid helium-4 were thought to originate purely from quantum mechanics. The simulation results showed that the individual vortices were thinner and that there were more vortices in the high-resolution simulation than in the low-resolution simulation. This is reasonable because if we focus on the cross section at z = 0, the number of vortices observed on the cross section should be a few hundred according to Feynman's law in two-dimensional cases.⁷ In addition, the vortex centers of quantum vortices are essentially of the order of angstroms (although quantum vortices are thought to interact even at large distances). Therefore, the results of the high-resolution simulation were closer to the theoretical solution.

As mentioned previously, the usual two-fluid model and our fluid model have the same form of equations, but their interpretations differ. Usually, the two components–superfluid and normal flow–overlap, and the interaction between them is interpreted as occurring only through the mutual friction term. In other words, these two components coexist in a single region. This is a quantum mechanical picture. However, our model considers liquid helium as an inviscid fluid in a continuum, with a viscous component remaining at finite temperatures. The two components interact as a multiphase flow under certain conditions; however, the nature of a multiparticle system appears locally. In this case, the temperature gradient along the horizontal direction was considered. In addition, a local interaction force proportional to the velocity difference, which is often used in particle–fluid coupling, acts between the two-fluid particles (we also call this the mutual friction force because it has the same form as the quantum mechanical mutual friction force). Thus, multiphase flow and particle-fluid coupling are often modeled according to classical fluid mechanics. Because the fluid particles are only fragments of continuum mechanics, they do not overlap and exhibit an excluded volume effect on each other. In summary, unlike the typical two-fluid model, our two-fluid model, which is based on classical fluid mechanics, is a one-fluid system. In this sense, it makes sense that high-resolution simulations approach the quantum mechanical picture more closely than low-resolution simulations. In multiparticle simulations, the computational resolution is determined by the diameter of the fluid particles. As the computation increases in resolution, the fluid particles decrease in size. Simultaneously, the excluded volume effect of the fluid particles also decreased; therefore, the particles were more likely to slip through and overlap with each other. Consequently, our method describes a one-fluid system at a low resolution, in which the two components are separated to describe a single fluid. However, as the resolution increases, the system is expected to become a two-fluid system. In other words, we can expect our two-fluid model to be smoothly connected and equal to the two-fluid model. That is, in this problem, the convergence of the calculations coincides with the convergence of the physical phenomena.

In our model, the fluid particles were subjected to pressure gradient forces, temperature gradient forces, shear viscosity, rotational viscosity, and mutual friction forces. The shear and rotational viscous forces act only between the components of the viscous fluid. In a typical case, the distribution of these forces acting on the fluid particles was approximately 72.9% pressure gradient force, 8.28% temperature gradient force, 18.07% shear viscous force, and approximately 0.71% rotational viscous force. The mutual friction force was approximately 0.01%–0.04% that of the rest. Although the rotational viscous force is less than 1% of the force applied to the fluid particles, it is understandable that it causes an independently spinning vortex. The mutual frictional force is negligible but may act as a disturbance, and which force induces which physical phenomena should be investigated in future studies.

Several experimental studies^17,73 on rotated cryogenic liquid helium-4 have reported peculiar phenomena compared to those observed in rotation problems in classical hydrodynamics. Specifically, multiple quantum vortices emerge and spin around their axes parallel to the vertical center axis of the cylindrical vessel, eventually forming a regular lattice. The lattices rotate around the cylindrical axis while maintaining a constant relative position, which is called rigid-body rotation. Of particular interest is the emergence of multiple independent spinning nondissipative vortices; this phenomenon in a 3D problem was numerically reproduced for the first time in this study. To clarify the implications of the findings of this study and for future tasks, we investigated the dependence of the number and average size of vortices on the numerical resolution, that is, the number of fluid particles, at two different elapsed times. The number of vortices was measured in the following manner. Because vortices exhibit vertical structures, we measured the number of vortices that crossed the z = 0 plane. We identified and labeled the vortices from the intensity distribution of the rotational forces exerted by the respective vortices on the z = 0 plane using the ImageJ software.^74–76  Figures 8(a) and 8(b) present this process; Fig. 8(a) is an example snapshot of the volume rendering of the simulation, showing the intensity of the rotational force at time of 0.32 s in the lower resolution simulation with 2.4 × 10⁶ particles; this subfigure corresponds to Fig. 4(c); Fig. 8(b) shows the scalar distribution of the intensity of the rotational force in the z = 0 plane in Fig. 8(a). Figures 8(a) and 8(b) clearly show that the edges of each vortex can be identified because the intensity of the rotational force is maximized at the edges of the vortices. Finally, Fig. 8(c) shows the correlation between the numbers of observed vortices and particles, and Fig. 8(d) shows the relationship between the average vortex size and the number of particles at 0.32 and 0.79 s.

The significant observations are summarized as follows: First, vortices tend to move freely or separately, absorbing neighboring smaller vortices and changing into a smaller number of larger vortices; the total number of vortices decreases with time. In brief, local interactions among vortices are classical. In practice, fewer vortices were observed, as shown in Fig. 8(c) for the triangular symbol (at 0.79 s) compared to that for the cross symbol (at 0.32 s). However, it should be emphasized that independent vortices rarely disappear unlike in ordinary classical fluids. In summary, multiple spinning, independent, and non-dissipative vortices emerged for horizontally rotated 3D cryogenic helium-4 at a lower resolution; however, local interactions among the vortices were classical, except for their non-dissipative nature.

However, this phenomenon was observed less frequently as the resolution increased. In Fig. 8(d), as the resolution increases from low (2.4 × 10⁶ particles) to medium (8.62 × 10⁶ particles) to high (19.6 × 10⁶ particles), the average size of the vortices decreases and no significant difference in size is observed between two different cases of elapsed time. Specifically, the average size of the vortices can be considered to reach a steady state over time. First, considering that the proposed method does not introduce constraints on the topology changes in the vortices, the interactions between the vortices inevitably exhibit classical fluidity. Nevertheless, as the resolution increases, the vortices become smaller. This can be described as follows. In our simulations, the diameter of the fluid particles decreased as the resolution increased, but the value of the parameter $Cω$⁠, which determines the magnitude of the angular velocity, was set to be inversely proportional to the size of the particle such that the rotational force at the edge of the individual fluid particles, expressed as their product as $r|ω0|$⁠, remained the same. This ensures that the torque (rotational force) on each fluid particle is constant regardless of the resolution. Consequently, although the number of particles required to form a vortex remained constant, the average size of the vortices decreased because the size of the constituent particles decreased. Then, the smaller sized vortices with the same magnitude of rotational force are distributed over the entire domain. Because the average size of the vortices decreases, the frequency of vortex–vortex interactions decreases; a higher number of multiple rotating vortices with small sizes can exist compared to the low-resolution case.

As we will see later, the average size of the vortices at 19.6 × 10⁶ particles is comparable to that of the average vortex–vortex distance derived from Feynman's rule in quantum hydrodynamics; from a local approximation point of view, it can be considered approximately equivalent to the velocity field generated by each quantum vortex. From a phenomenological perspective, the fact that a number of multiple independent vortices can exist with smaller sizes comparable to the velocity field of quantum vortices as the resolution increases suggests that vortices become more quantum mechanical as the resolution increases. In other words, our two-fluid model might directly describe the microscopic behavior as the resolution increases. There are two possible explanations for this observation. The first possibility is that this claim is incorrect, and that artificial control of the vortex topology is still required. In quantum mechanics, circulation is quantized. Therefore, we may need to introduce a constraint such as the Schwarz equation to maintain the circulation or each vortex constant, leading to a change in the system into a quantum fluid system. The second possibility is that such a constraint can be introduced automatically. According to the Young–Laplace equation, the surface tension of a classical fluid is proportional to its curvature. As the vortex becomes smaller, its surface tension increases, which may solidify the vortex and act as a type of topological control. To validate this, simulations with the same physical conditions and at higher resolutions must be performed.

Based on the results of this study, the vortices have more opportunities to interact with each other over longer periods. Consequently, the vortices gradually integrate and decrease to the same degree as the ideal number of vortices calculated by Feynman's rule, which is equal to 284.7 in the corresponding 2D case (refer to Ref. 7 for specific calculations) with 19.6 × 10⁶ particles in Fig. 8(c). However, even in this case, the number of vortices decreased with time. Thus, the system did not reach a steady state with respect to the number of vortices in the time direction. Among the two possibilities, in the former, a comparison with the theoretical value of Feynman's rule must be discussed in detail after the corresponding 3D physical model of circulation control is introduced, and the number of vortices becomes stable in the time direction. By contrast, in the latter, higher resolution simulations are required. In Fig. 8(c), the curve with the triangular symbol is exponential relative to the computational resolution, whereas the curve with the cross symbol is linear with respect to the computational resolution. Therefore, we can expect to reach a stationary state in the time direction when using tens to 100 × 10⁶ particles for the entire domain without introducing such a topology control model. In the former case, developing a 3D vortex dynamics model based on classical fluid mechanics as an alternative to Schwarz's rule^13,14 is another approach. Our results for the 2D problem have been presented in previous studies,^7,8 where we proposed a model that balances the Magnus force generated by each rotating vortex with the repulsive force as a vortex–vortex interaction, which prevents the vortices from moving closer than necessary. Specifically, by calculating the vortex circulation value and introducing a numerical model in which the vortex dynamics model becomes effective when the circulation exceeds a certain criterion, we successfully controlled the equal arrangement and circulation of vortices.^7,8 Hence, in principle, it is possible to develop similar dynamic vortex models for 3D problems. However, several issues need to be resolved, such as the physical correspondence between the classical vortex dynamics model and Schwarz's rule,^13,14 including whether they are in one-to-one correspondence. Moreover, in future research, our 2D vortex dynamics model should be extended to changes in 3D topology.

For 19.6 × 10⁶ particles, 500 fluid particles were arranged in each horizontal direction, as mentioned previously. Because the diameter of the disk was 0.2 cm (⁠ $=2.0×10−3m$⁠), the diameter of a particle is $4×10−6m$ or 4 μm. If we assume ten particles per vortex in each horizontal direction, the size of a single vortex would be approximately 40 μm. In contrast, the scale of a vortex core is on the order of angstroms. However, because the effective range of a vortex is determined by the Biot–Savart law, which describes the velocity field around the vortex, the area of influence of the vortex is considered to be its effective size. In our case, the average vortex distance was calculated using Feynman's rule, which yielded approximate $10−4m$⁠.⁷ The corresponding number of vortices is shown by the dotted line in Fig. 8(c). It is reasonable to understand that the vortex observed does not correspond to an actual quantum vortex core, but to the effective velocity field created by the quantum vortex.

In summary, the next task is to control the circulation of vortices by ensuring an appropriate vortex–vortex interaction to stabilize the number of vortices in the time direction. As described above, this may be achieved by introducing an artificial circulation control model by referring to Schwarz's rule, or by simply increasing the resolution, expecting that an increase in the strength of the surface tension would solidify the vortex. We regard the latter case as promising for the following reasons. In Fig. 8(c), the curve with the triangular symbol is exponential relative to the computational resolution, whereas that with the cross symbol is linear with respect to the computational resolution. Extrapolating these results enables us to expect that the vortices may reach an approximately stationary state in the time direction when tens to 100 × 10⁶ particles are used in the entire domain. Another approach is to develop a 3D vortex dynamics model based on classical hydrodynamics, as an alternative to Schwarz's rule. We achieved some success in the 2D calculations using models representing this effect in the Magnus force and vortex repulsion models,⁷ as described above. However, several issues need to be resolved, such as the physical correspondence between the classical vortex dynamics model and Schwarz's rule and how to extend our 2D vortex dynamics model to changes in 3D topology. Nevertheless, this approach should be examined after checking whether the second approach is valid.

From an engineering perspective, our goal is not to propose a first-principles model, but to enable the simulation of practical bulk liquid helium; we only have to capture the dynamic interactions among vortices accurately enough to reproduce macroscopic quantum effects. In practice, a direct comparison between the internal structure within the healing length of the quantum vortex and the vortex in our method requires at least several hundred times more particles in each direction, even when using a coarse-grained approach. However, this was beyond the scope of this study. At present, this is the first study to numerically reproduce the dynamics of rotated cryogenic 3D liquid helium-4 in bulk form using a classical hydrodynamic approach. Accordingly, it is scientifically significant that we observed multiple spinning nondissipative vortices in the 3D problem of horizontally rotated cryogenic liquid helium-4.

In this paper, we report a 3D simulation of the rotating liquid helium-4 using a two-fluid model with spin-angular momentum conservation using SPH. The conventional two-fluid model was derived by combining the nonlinear Schrödinger equation for bosons with the thermodynamic Gibbs–Duhem equation. In short, the usual two-fluid model is obtained by making macro-corrections to the microscopic relational equations. In contrast, our model is derived from the particle approximation of an inviscid fluid with residual viscosity in part. Despite the fully classical mechanical picture, the resulting system equations were consistent with those of the conventional two-fluid model. Specifically, we consider the bulk liquid helium-4 to be an inviscid fluid, assuming that the viscous fluid component remains at finite temperatures. As the temperature decreased, the amount of the viscous fluid component decreased, ultimately becoming a fully inviscid fluid at absolute zero. Weak compressibility is assumed to express the volume change because some helium atoms do not render fluid owing to BECs or change states owing to local thermal excitation. Using explicit SPH (smoothed-particle hydrodynamics), one can solve the governing equations for an incompressible fluid, simultaneously reproducing density fluctuations and describing the fluid in a many-particle system. We assume the following fluid–particle duality: a hydrodynamic interfacial tension between the inviscid and viscous components or a local interaction force between two types of fluid particles. The former can be induced in the horizontal direction when non-negligible non-uniformity of the particles occurs during forced 2D rotation, and the latter is non-negligible when the former is negligible. We performed a large-scale simulation of 3D liquid helium forced to rotate horizontally using 32 GPUs. Compared with the low-resolution calculation using 2.4 × 10⁶ particles, the high-resolution calculation using 19.6 × 10⁶ particles showed spinning vortices close to those of the theoretical solution. In the future, the introduction of quantum mechanical recombination among vortices will capture the behavior of quantum fluids. We obtained a promising venue to establish a practical simulation method for bulk liquid helium-4.

This study was supported by JSPS KAKENHI (Grant Numbers 22K14177, JST PRESTO, and JPMJPR23O7). The authors thank Editage (www.editage.jp) for the English language editing. The author would also like to express his gratitude to his family for their moral support and encouragement.

The authors have no conflicts to disclose.

Satori Tsuzuki: Conceptualization (lead); Data curation (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 sharing is not applicable to this article as no new data were created or analyzed in this study.

To facilitate understanding, we provide a schematic list of the equations for the derivation of the Navier–Stokes equations with spin angular momentum conservation in Fig. 9. We also provide a schematic summary of SPH in Fig. 10.