Particle approximation of the two-fluid model for superfluid 4 He using smoothed particle hydrodynamics
aa r X i v : . [ c ond - m a t . o t h e r] A p r Particle approximation of the classical two-fluid model for superfluid He usingsmoothed particle hydrodynamics Satori Tsuzuki ∗ Research Center for Advanced Science and Technology,The University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8904, Japan (Dated: March 2020)This paper presents a finite particle approximation of the classical two-fluid model for liquid He.We propose to use smoothed particle hydrodynamics (SPH), a well-established Lagrangian particlemethod in astrophysics. The model parameters are confirmed to reproduce the phenomenologicalcharacteristics of superfluid He qualitatively. We have thus pioneered a new paradigm of the La-grangian particle mechanics, which will hopefully stimulate further numerical studies of superfluidsin condensed matter physics.
I. INTRODUCTION
The bizarre behavior of liquid helium 4 (hereinafter, He) has attracted the attention of condensed matterphysicists for many years; however, the detailed dynam-ics of this fluid remains a mystery even now, eighty yearsafter its discovery.Figure 1 reviews the basic properties of He through aschematic phase diagram. In 1937, P. L. Kapitza foundthat the viscosity of liquid He decreases upon coolingand reaches almost zero below the λ -point (labelled aspoint c in Fig. 1). The liquid He was later confirmed todisplay viscosity in the high-temperature area over the λ -line, and to lose it in the low-temperature area belowthe λ -line. Once lack of viscosity is achieved in this lattertemperature range, liquid He displays many unexpectedproperties that characterize it as a “superfluid.” In con-trast, liquid He within the high-temperature area overthe λ -line is still classified as a “normal fluid.”The superfluidity condition is characterized by proper-ties such as the following two examples. Superfluid Hecan penetrate a high-dense capillary filter (e.g., gypsum),thanks to the so-called “superleak effect”, whereas this isnot possible for normal fluid He. The “film flow effect”,instead, takes place when superfluid He spontaneouslycomes out of its vessel, due to van der Waals forces be-tween helium atoms and walls being stronger than theforces between He atoms themselves.Many physicists have made efforts to explain these in-teresting properties of superfluids at both the ontologicaland phenomenological level. One year after the discoveryby P. L. Kapitza, F.W. London explained the phase tran-sition of liquid He from a quantum mechanical point ofview, comparing it to a type of Bose-Einstein condensa-tion (BEC) [1]. London’s theory has now been validatedby many experiments. Meanwhile, the “two-fluid model”proposed by L. Tisza[2] and L. Landau [3], which assumescoexistence of superfluidity and normal fluidity in liquid He, has become the most acknowledged phenomenolog-ical description. In this model, the total mass density ρ ∗ The University of Tokyo: [email protected]
FIG. 1. Schematic phase diagram of He. Points a , b , c ,and d correspond to the following (p, T) values, respec-tively: (25bar , , . . , . . , . is expressed as follows: ρ = ρ n + ρ s , (1)where ρ n and ρ s are the mass densities of the normalfluid and of the superfluid, which respectively satisfy thelaws of mass and entropy conservations as follows: ∂ρ∂t + ∇ · ( ρ n v n + ρ s v s ) = 0 , (2) ∂∂t ( ρσ ) + ∇ · ( ρσ v n ) = 0 . (3)Here, v n and v s are the velocities of the normal fluidand of the superfluid, respectively, while σ is the entropydensity.A remarkable achievement of the two-fluid model isthat it connects the macroscopic perspective of thermo-dynamics with quantum mechanics, so that the super-fluid can be described within a framework of basic me-chanics. By solving the Gross-Pitaevskii equation (a non-linear Schr¨odinger equation for boson particles) and theGibbs-Duhem equation simultaneously, we obtain the fol-lowing equations of motions for liquid He under a phe-nomenological perspective [4]: ρ s D v s D t = − ρ s ρ ∇ P + ρ s σ ∇ T , (4) ρ n D v n D t = − ρ n ρ ∇ P − ρ s σ ∇ T + η n ∇ v n , (5)where D {·} /Dt represents a material derivative, v n and v s are the velocities of normal fluid and superfluid. η n isthe viscosity of the normal fluid. T is temperature, P ispressure, and σ is entropy density, respectively.Finite-approximation of the phenomenological equa-tions of motion shown in Eq. (4) and Eq. (5) represents anessential step towards direct numerical simulations. How-ever, until today, only a few studies have followed such anapproach; this is primarily because most previous workstake a non-phenomenological point of view (e.g., numeri-cal simulations of nonlinear Schr¨odinger equations [5, 6],and quantum Monte Carlo calculations [7, 8] aiming toexamine the static responses [9] or to explore quantumvortices [10–12]). Although previous works [13–17] intro-duce direct numerical simulations of the two-fluid model,they focus on mesh-based discretization; when simulatingliquid He, some studies [13–15] developed finite differ-ence approximations, and others [16, 17] devised novelfinite element approximations. Most importantly, noneof the previous studies have reported a finite particleapproximation of these governing equations, despite itsgreat potential to realize flexible and practical simula-tions.The main purpose of this study is to present a schemeof finite particle approximation of the two-fluid model.We propose to discretize the two-fluid model of superfluid He by smoothed particle hydrodynamics (SPH) [18], awell-established Lagrangian particle approximation par-ticularly popular in the field of astrophysics. By com-bining theoretical methodologies typical of two differ-ent branches of physics (low-temperature physics and as-trophysics), we aim to pioneer a new paradigm of La-grangian particle mechanics targeting superfluids.The remainder of this article is structured as follows.Section 2 describes the concept of SPH, provides a sum-mary of wave equations in superfluid He, and discussesthe expression of thermodynamic parameters in the two-fluid model, in preparation for Section 3. In Section 3,we derive a finite particle approximation of liquid He.Section 4 examines the discretized parameters of the de-rived approximation model. Section 5 summarizes ourresults and concludes the paper.
II. PREPARATIONSII.1. A brief overview of SPH
The fundamental concept of SPH is the approximationof the Dirac delta function δ in integral form by a dis-tribution function W , which is called smoothed kernel FIG. 2. Schematic of a kernel function. function, as follows [18]: φ ( r ) = Z Ω φ (´ r ) δ ( r − ´ r ) d ´ r ≃ Z Ω φ (´ r ) W ( r − ´ r , h ) d ´ r , (6)where φ is a physical value at the position r in a domainΩ and the parameter h is called a kernel radius. Thesmoothed kernel function W must have at least the fol-lowing four characteristics; (a) W converges to the deltafunction δ as h gets close to 0 (lim h → W = δ ). (b) W satisfies the normalization condition ( R W d r = 1). (c) W satisfies the symmetricity W ( − r ) = W ( r ). (d) W is suchthat the value of W becomes 0 at a distance kh , which iscalled a compact support condition. Figure 2 representsa schematic view of a kernel function. An example of asmoothed kernel function W is the Gaussian, which issimply expressed as W ( r ) := Ch d e − ( r − ´ r ) /h , (7)where d is the dimension ( d = 1 , , or 3), and C is anormalization factor. The Gaussian kernel is frequentlyused in the field of astrophysics. In the field of incom-pressible SPH, in contrast, polynomial functions are oftenused [19, 20].Equation (6) is further discretized on the basis of sum-mation approximation using the position r j , the infinites-imal volume ∆ V j , density ρ j and mass m j of the j thparticle as follows. φ ( r i ) ≈ N p X j φ ( r j ) W ( | r i − r j | , h )∆ V j = N p X j φ ( r j ) ρ j m j W ( | r i − r j | , h ) , (8)where m j = ρ j ∆ V j and N p is the number of fluid par-ticles to approximate the system. The gradient and theLaplacian of the physical value φ ( r i ) are axiomaticallyobtained from Eq. (8) using a vector analysis, as follows: ∇ φ ( r i ) = ρ i N p X j m j n φ ( r i ) ρ i + φ ( r j ) ρ j o ∇ W ij , (9) ∇ φ ( r i ) = 2 N p X j m j ρ j φ ( r i ) − φ ( r j ) | r i − r j | ( r i − r j ) · ∇ W ij , (10)where W ij = W ( | r i − r j | , h ). Here, we have usedthe mathematical relationship ∇ f = ρ [ ∇ ( f /ρ ) +( f /ρ ) ∇ ρ ] to derive Eq. (9), so that the gradient op-erator has the symmetricity of pressures between the i thand the j th particles. For more details of the operatorsin SPH, refer to [21]. II.2. Wave equations of superfluids
According to the literature [3], the governing equationsof Eq. (2) to Eq. (5) yield wave propagation equations re-garding density ρ and entropy S , which produce multipletypes of sound waves. Under a first approximation thatneglects the viscosity term, the wave equations give twodifferent sound waves expressed as follows [22, 23]: c = (cid:18) ∂P∂ρ (cid:19) / S , (11) c = (cid:18) ρ s ρ n T S C v (cid:19) / , (12)where c and c are called the first sound velocity andthe second sound velocity, respectively. C v representsthe specific heat at constant volume. Equation (12) canbe rewritten by using the number of He atoms N andthe mass of a He atom M as ρ s ρ n = ( C v c ( N M ) T σ ) , (13)where we have introduced a basic relation between theentropy density σ and the entropy S of σ = S/ ( N M ). II.3. Expression of thermodynamic variables in anelementary excitation model
As a consequence of the two-fluid model by L. Landau,it was found that the elementary excitation of liquid Hederives from two components, “phonon” and “roton” [3].The pressure P of the liquid He can thus be expressedas follows [24]: P = P ph + P rot . (14)Here, let us denote the ratio of a circle’s circumferenceas π , and the speed of sound as c . We also introduce theconstant values of µ , p , and ∆, each of which is spec-ified as 1 . × − g, 2 . × − gcms − , and 8 . T ≪
93 K,Eq. (14) can be approximated as [24, 25]: P ≃ π ( k B T ) ~ c + p p ( πµ/ π ~ ( k B T ) / e − ∆ /T , (15)where k B represents the Boltzmann constant, and ~ isthe reduced Planck’s constant. Note that Eq. (15) issimilar to Eq. (5) in the literature [24]. The entropy S is obtained from the first-order partial derivatives of P with respect to temperature T as S := (cid:18) ∂P∂T (cid:19) (16) ≃ π k B T ~ c + (cid:18) k B π (cid:19) / √ µp ∆ ~ e − ∆ /T √ T , (17)where we have ignored the term with the higher order inthe derivation of the second term in Eq. (17).Accordingly, the entropy density σ is obtained usingthe relation σ = S/ ( N M ) as σ = ζT + ξ e − ∆ /T T / , (18) ζ := 2 π k B ~ c N M , (19) ξ := (cid:18) k B π (cid:19) / √ µp ∆ ~ N M . (20)Here, both ζ and ξ are constant parameters, because theircomponents of N, M, π, c, µ, p , ∆ , k B , and ~ are all con-stant and can be given as known parameters. Now wehave prepared all the parameters needed for the deriva-tion of a particle approximation of liquid He.
III. A PARTICLE APPROXIMATION FORSUPERFLUIDS
First, Eq. (4), the motion equation for superfluid com-ponents, is rewritten as:D v s D t = − ρ ∇ P + σ ∇ T . (21)By using Eq. (9), we obtain the particle discretization ofthe right-hand side of Eq. (21), with respect to the i thparticle, as follows: * − ρ ∇ P + i = − X j m j n P i ρ i + P j ρ j o ∇ W ij , (22) * σ ∇ T + i = θ si X j m j n T i ρ i + T j ρ j o ∇ W ij , (23)where θ si is represented using Eq. (18) as θ si = ρ i σ i . (24)Here, the superscript of θ xi represents the first letter of“superfluids” or “normal fluids” ( x = s, n). σ i is a specificcase of Eq. (18) with T = T i .Second, Eq. (5), the motion equation for normal fluidcomponents, is rewritten as:D v n D t = − ρ ∇ P − ρ s ρ n σ ∇ T + η n ρ n ∇ v n . (25)The first term of the right-hand side in Eq. (25) can bediscretized similar to the case of Eq. (22). Meanwhile, thediscretized expression of the second term with respect tothe i th particle is obtained by using Eq. (9) as follows. * − ρ s ρ n σ ∇ T + i = θ ni X j m j n T i ρ i + T j ρ j o ∇ W ij , (26)where the breakdown of θ ni is given as θ ni = χ ρ i T i σ i , (27) χ := C v c ( N M ) . (28)In the end, a substitution of Eq. (10) into the third termdirectly leads to: * η n ρ n ∇ v n + i = 2 ν n X j m j ρ j v ( i ) n − v ( j ) n | r i − r j | ( r i − r j ) · ∇ W ij , (29)where ν n is kinetic viscosity that is equal to η n /ρ n . Nowthe right parts of Eq. (21) and of Eq. (25) are appro-priately discretized. The material derivative of velocitiesin the left parts of Eq. (21) and of Eq. (25) can be dis-cretized using an explicit time-integrating scheme, e.g.,the Verlet algorithm [27], in accordance with the usualmanner in SPH simulations. Consequently, we have suc-cessfully derived a particle approximation of the govern-ing equations of liquid He.
IV. DISCUSSION
In the discretized forms discussed in the previous sec-tion, the variables for the i th particle are given as aset of ( r i , v ( i ) n , v ( i ) s , m i , ρ i , σ i , P i , T i ). The constant pa-rameters characterizing the system are given as a set of( N, M, C v , c, ν n ). The exact constant values that are in-dependent of the systems are ( µ, p , ∆ , k B , ~ ). In incom-pressible SPH simulations, pressure P i is calculated fromdensity ρ i by solving a one-to-one thermodynamic corre-spondence between P i and ρ i given by the Tait equationof state [28]. Furthermore, the entropy density σ is givenas a function of temperature T , as shown in Eq. (18), andthe positions r i are calculated from the velocities. Thus,the independent variable of the i th particle is given as( v ( i ) n , v ( i ) s , m i , ρ i , T i ). On the other hand, the remain-ing parameter, the second sound velocity c should besemiempirically determined according to literature exem-plified by [29, 30]. Note that it is known that c emergeswhen the system reaches temperatures lower than thecritical temperature of approximately 2 .
17 K, increases
FIG. 3. Dependences of (a) the entropy S and (b) ( T S ) − onthe temperature T . as cooling of the system increases, and finally reachesapproximately 135 ms − after experiencing a sluggish-ness [30].Let us discuss the mechanical effect of the discretizedparameters of θ si and θ ni . It can be said from Eq. (23)and Eq. (26) that each of θ si and θ ni functions as a param-eter that determines the degree of the effect of tempera-ture gradient forces. The entropy density σ i is dominantin the behavior of θ si , because the density ρ i fluctuatestypically less than one percent and at most a few per-cent compared to its averaged density in incompressibleSPH [31, 32]. Besides, the behavior of σ corresponds tothat of entropy S because the value of N M is constant;thus, θ si is proportional to the entropy S .Figure 3(a) shows the dependence of entropy S ontemperature T . The following fact is suggested fromFig. 3(a); The effect of the temperature gradient forcecontrolled by θ si on the superfluid component diminishesas the temperature T decreases, and it eventually reacheszero at the absolute zero. We can say from this result thatthe parameter θ si well reflects the characteristic that su-perfluids have almost zero or a very small temperaturegradient corresponding to the extremely high heat con-ductivity.Let us assume in this discussion that c is constant.In this ideal case, the parameter χ becomes constant,and thus θ ni becomes proportional to ( T i σ i ) − . Also,the behavior of σ corresponds to that of entropy S assimilar to the case of θ si . After all, the behavior of θ ni is determined by that of ( T S ) − . Figure 3(b) showsthe dependence of ( T S ) − on the temperature T witha semilogarithmic scale. The effect of θ ni on the fluid Heis confirmed to greatly increase as the temperature T i decreases. Note that the parameter θ ni is the discretizedformula of Eq. (13), expressing the ratio of superfluiddensity to normal density, which generates the secondsound wave. Figure 3(b) suggests that the effect of ρ s /ρ n increases as the temperature T decreases on the premisethat c is constant; the effect is more amplified in actualcases because the velocity c itself increases as mentionedabove [30]. In any case, both θ si and θ ni are confirmed toreflect the phenomenological behavior of superfluid Hequalitatively.
V. CONCLUSIONS
In this paper, we have theoretically derived a finite par-ticle approximation of the two-fluid model of superfluid He using smoothed particle hydrodynamics. Our ap-proximation model builds on previous results, and modelparameters are found to reflect the phenomenologicalfacts of He at the first approximation level.The author understands that there is still room forimprovement before practical use, because it is possiblethat the parameter θ ni is too sensitive for numerical sim-ulations; in that case, it would be necessary to decreasethe sensitivity of this parameter. The effect of the param-eters θ si and θ ni on the stability of numerical simulationsshould also be investigated in future work. Nevertheless,the author believes this theoretical work is meaningful for revealing the overall picture of particle approximation ofthe two-fluid model of superfluid He. The author hopesthat the present study may stimulate further numericalstudies of superfluids in condensed matter physics.
ACKNOWLEDGEMENT [1] Bose, Plancks gesetz und lichtquantenhypothese,Zeitschrift f¨ur Physik , 178 (1924).[2] L. TISZA, Transport phenomena in helium ii, Nature , 913 (1938).[3] L. Landau, Theory of the superfluidity of helium ii, Phys.Rev. , 356 (1941).[4] M. Tsubota, K. Fujimoto, and S. Yui, Numerical stud-ies of quantum turbulence, Journal of Low TemperaturePhysics , 119 (2017).[5] L. Lehtovaara, T. Kiljunen, and J. Eloranta, Efficient nu-merical method for simulating static and dynamic prop-erties of superfluid helium, Journal of ComputationalPhysics , 78 (2004).[6] I. Danaila and F. Hecht, A finite element method withmesh adaptivity for computing vortex states in fast-rotating bose-einstein condensates, Journal of Compu-tational Physics , 6946 (2010).[7] D. CEPERLEY and B. ALDER, Quantum monte carlo,Science , 555 (1986).[8] J. Casulleras and J. Boronat, Unbiased estimators inquantum monte carlo methods: Application to liquid He, Phys. Rev. B , 3654 (1995).[9] S. Moroni, D. M. Ceperley, and G. Senatore, Staticresponse from quantum monte carlo calculations, Phys.Rev. Lett. , 1837 (1992).[10] S. A. Vitiello, L. Reatto, G. V. Chester, and M. H. Kalos,Vortex line in superfluid He: A variational monte carlocalculation, Phys. Rev. B , 1205 (1996).[11] L. Vranjeˇs, J. Boronat, J. Casulleras, and C. Cazorla,Quantum monte carlo simulation of overpressurized liq-uid He, Phys. Rev. Lett. , 145302 (2005).[12] D. E. Galli, L. Reatto, and M. Rossi, Quantum montecarlo study of a vortex in superfluid He and search for avortex state in the solid, Phys. Rev. B , 224516 (2014).[13] Y. S. Ng, J. H. Lee, and W. F. Brooks, Numerical mod-eling of helium-ii in forced flow conditions, Journal ofThermophysics and Heat Transfer , 203 (1989).[14] M. Murakami and K. Iwashita, Numerical computationof a thermal shock wave in he ii, Computers & Fluids ,443 (1991).[15] O. C. Idowu, D. Kivotides, C. F. Barenghi, and D. C.Samuels, Numerical Methods for Coupled Normal-Fluidand Superfluid Flows in Helium II (Springer Berlin Hei- delberg, Berlin, Heidelberg, 2001), pp. 162–176.[16] L. Bottura, C. Darve, N. A. Patankar, and S. W. V.Sciver, A method for the three-dimensional numericalsimulation of superfluid helium, Journal of Physics: Con-ference Series , 012008 (2009).[17] C. M.-T. Darve, Phenomenological and numerical stud-ies of helium ii dynamics in the two-fluid model, a disser-tation at the Northwestern university, Evanston, Illinois(2011).[18] R. A. Gingold and J. J. Monaghan, Smoothed particlehydrodynamics: theory and application to non-sphericalstars, Monthly notices of the royal astronomical society , 375 (1977).[19] M. Desbrun and M.-P. Gascuel, Smoothed particles: Anew paradigm for animating highly deformable bodies,in
Computer Animation and Simulation96 , pp. 61–76,Springer, 1996.[20] M. M¨uller, D. Charypar, and M. Gross, Particle-basedfluid simulation for interactive applications, in
Proceed-ings of the Eurographics symposium on Computer anima-tion , pp. 154–159, Eurographics Association, 2003.[21] J. J. Monaghan, Smoothed particle hydrodynamics, An-nual Review of Astronomy and Astrophysics , 543(1992).[22] R. J. Donnelly, The two-fluid theory and second soundin liquid helium, Phys. Today , 34 (2009).[23] E. LIFSHITZ and F. Leib, 12 - radiation of sound inhelium iireprinted from journal of physics, 8, part 2, 110,1944., in Perspectives in Theoretical Physics , pp. 177 –183, Pergamon, Amsterdam, 1992.[24] I. N. Adamenko, K. E. Nemchenko, I. V. Tanatarov, andA. F. G. Wyatt, Pressure of thermal excitations in su-perfluid helium, Journal of Physics: Condensed Matter , 245103 (2008).[25] A. Schmitt, Introduction to superfluidity, Lect. NotesPhys (2015).[26] K.-H. Bennemann and J. B. Ketterson, Novel superflu-ids volume 1 (OUP Oxford, 2013).[27] L. Verlet, Computer ”experiments” on classical fluids. i.thermodynamical properties of lennard-jones molecules,Phys. Rev. , 98 (1967).[28] J. J. Monaghan, Smoothed particle hydrodynamics, Re-ports on Progress in Physics , 1703 (2005). [29] J. Maynard, Determination of the thermodynamics ofhe ii from sound-velocity data, Phys. Rev. B , 3868(1976).[30] J. Wilks and D. S. Betts, An introduction to liquid he-lium. 2, (1987).[31] J. Monaghan, Simulating free surface flows with sph,Journal of Computational Physics , 399 (1994).[32] Nomeritae, E. Daly, S. Grimaldi, and H. H. Bui, Explicitincompressible sph algorithm for free-surface flow mod-elling: A comparison with weakly compressible schemes,Advances in Water Resources97