Relativistic Lattice Boltzmann Model with Improved Dissipation

0.5 MB·English

by  M. Mendoza

Relativistic Lattice Boltzmann Model with Improved Dissipation
M. Mendoza,1,∗ I. Karlin,2,† S. Succi,3,‡ and H. J. Herrmann1,4,§
1 ETH Zu¨rich, Computational Physics for Engineering Materials,
Institute for Building Materials, Schafmattstrasse 6, HIF, CH-8093 Zu¨rich (Switzerland)
2ETH Zu¨rich, Department of Mechanical and Process Engineering,
Sonneggstrasse 3, ML K 20, CH-8092 Zu¨rich (Switzerland)
3Istituto per le Applicazioni del Calcolo C.N.R., Via dei Taurini, 19 00185, Rome (Italy),
and Freiburg Institute for Advanced Studies, Albertstrasse, 19, D-79104, Freiburg, (Germany)
3 4Departamento de F´ısica, Universidade Federal do Ceara´,
1 Campus do Pici, 60455-760 Fortaleza, Ceara´, (Brazil)
0 (Dated: January 16, 2013)
2
We develop a relativistic lattice Boltzmann (LB) model, providing a more accurate description

n of dissipative phenomena in relativistic hydrodynamics than previously available with existing LB
a schemes. The procedure applies to the ultra-relativistic regime, in which the kinetic energy (tem-
J perature) far exceeds the rest mass energy, although the extension to massive particles and/or low

5 temperatures is conceptually straightforward. In order to improve the description of dissipative ef-
1 fects, theMaxwell-Ju¨ttner distribution is expandedin a basis of orthonormal polynomials, so as to
correctlyrecoverthethirdordermomentofthedistributionfunction. Inaddition,atimedilatation

] is also applied, in order to preserve the compatibility of the scheme with a cartesian cubic lattice.
h
To the purpose of comparing the present LB model with previous ones, the time transformation
p
is also applied to a lattice model which recovers terms up to second order, namely up to energy-
-
p momentumtensor. Theapproach isvalidatedthroughquantitativecomparison betweenthesecond
m andthirdorderschemeswithBAMPS(thesolution ofthefullrelativistic Boltzmannequation),for
moderatelyhighviscosityandvelocities,andalsowithpreviousLBmodelsintheliterature. Excel-
o
lentagreementwithBAMPSandmoreaccurateresultsthanpreviousrelativisticlatticeBoltzmann
c
models are reported.
.
s
c
PACSnumbers: 47.11.-j,12.38.Mh,47.75.+f
i
s Keywords:
y
h
p I. INTRODUCTION agreementwith the solution of the full Boltzmann equa-
[
tion as obtained by Bouras et al. using BAMPS (Boltz-

1 mann Approach Multi-Parton Scattering) [10, 11]. The
v Relativistic hydrodynamics and kinetic theory play a RLB makes use of two distribution functions, the first
3 major role in many forefronts of modern physics, from one modeling the conservation of number of particles,
2 large-scale applications in astrophysics and cosmology, and the second one, the momentum-energy conservation
4
to microscale electron flows in graphene [1–3], all the equation. The model was constructed by matching the
3
way down to quark-gluon plasmas [4–6]. Due to their first and second order moments of the discrete-velocity
.
1 strong non-linearity and, for the case of kinetic theory, distribution function to those of the continuum equilib-
0 high dimensionality as well, the correspondingequations rium distribution of a relativistic gas
3
1 areextremelychallengingevenforthemostpowerfulnu- In a subsequent work, Hupp et al.[12] improved the
: mericalmethods,letaloneanalytics. Recently,apromis- scheme by extending the equilibrium distribution func-
v
ing approach, based on a minimal form of relativistic tion for the number of particles, in such a way as to
i
X Boltzmann equation, whose dynamics takes place in a include second order terms in the velocity of the fluid,
r fully discrete phase-spaceandtime lattice, knownas rel- therebytamingnumericalinstabilitiesforhigherpressure
a ativisticlatticeBoltzmann(RLB),hasbeenproposedby gradients and velocities. However, the model was not
Mendoza et al. [7–9]. To date, the RLB has been ap- able to reproducethe rightvelocity and pressureprofiles
plied to the simulation of weakly and moderately rela- fortheRiemannprobleminquark-gluonplasmas,forthe
tivistic fluid dynamics, with Lorentz factors of γ 1.4, caseoflargevaluesofthe ratiobetweenthe shearviscos-
∼
where γ =1/ 1 v2/c2, c being the speed of light and ity and entropy density, η/s 0.5, at moderate fluid
v the speedofpthe−fluid. Thismodelreproducescorrectly speeds (v/c 0.6). ∼
shock waves in quark-gluon plasmas, showing excellent ∼
Inordertosetupatheoreticalbackgroundforthe lat-
tice version of the relativistic Boltzmann equation, Ro-
matschke et al. [13] developed a scheme for an ultrarel-
ativistic gas based on the expansion on orthogonalpoly-
Electronicaddress: [email protected]
∗ nomials of the Maxwell-Ju¨ttner distribution [14] and, by
Electronicaddress: [email protected]
†
Electronicaddress: [email protected] following a Gauss-type quadrature procedure, the dis-
‡
Electronicaddress: [email protected] crete version of the distribution and the weighting func-
§
2
tions was calculated. This procedure was similar to the It represents a genuine lattice Boltzmann discretization
oneusedforthenon-relativisticlatticeBoltzmannmodel of space and time, with no need of any interpolation
[15, 16]. This relativistic modelshowedverygoodagree- scheme, thereby avoiding the otherwise ubiquitous spu-
ment with theoretical data, althoughit was not compat- rious dissipation. The new lattice Boltzmann model is
ible with a lattice, thereby requiring linear interpolation shown to reproduce with very good accuracy the results
in the free-streaming step. This implies the loss of some of the shock-waves in quark-gluon plasmas, for moder-
keypropertiesofthestandardlatticeBoltzmannmethod, ately high velocities and high ratios η/s.
suchasnegativenumericaldiffusionandexactstreaming. The paper is organized as follows: in Sec. II we de-
Very recently, Li et al. [17] noticed that the equa- scribe in detail the model and the way it is constructed;
tion of conservation for the number of particles, recov- in Sec. III, we implement simulations of the Riemann
eredbytheRLB model[7,8],exhibitsincorrectdiffusive problem in order to validate our model and compare it
effects. Therefore,they proposedan improvedversionof withBAMPSandpreviousrelativisticlatticeBoltzmann
RLB, using a multi-relaxation time collision operator in models; finally, in Sec. IV, we discuss the results and
the Boltzmann equation, showing that this fixes the is- future work.
suewiththe equationfortheconservationofthe number
of particles. The generalized collision operator allows to
tune independently the bulk and shear viscosities, yield- II. MODEL DESCRIPTION
ing results for the Riemann problem closer to BAMPS
[10] when the bulk viscosity is decreased. However, the A. Symmetries of the relativistic Boltzmann
third order moment of the equilibrium distribution still equation
does not match its continuum counterpartand therefore
themodelstillhasproblemstoreproducehighη/s 0.5,
∼ To build our model, we start from the relativistic
for moderately high velocities, β = v/c = 0.6. Thus,
Boltzmannequationfortheprobabilitydistributionfunc-
while surely providing an improvement on the original
tion f:
RLB model, the work [17] did not succeed in reproduc-
ingthevanishingbulkviscosity,whichispertinenttothe p Uµ
ultra-relativisticgas,while allowingthe bulk viscosityto pµ∂µf =− cµ2τ (f −feq) , (1)
vary independently on the shear viscosity.
Note that in the much more studied case of the lat- where the local equilibrium is given by the Maxwell-
tice Boltzmannmodels forthe non-relativisticfluids, the Ju¨ttner equilibrium distribution [14],
question of the choice of the lattice with higher-order
symmetry requirements has only recently been solved, feq =Aexp( pµUµ/kBT) , (2)
−
in the framework of the entropy-compliant constriction
[18, 19]. However, the lattices (space-filling discrete- Intheabove,Aisanormalizationconstant,cthespeedof
velocitysets)foundinthatcasearetailoredtoreproduce light,andkB theBoltzmannconstant. The4-momentum
the moments of the non-relativistic Maxwell-Boltzmann vectors are denoted by pµ = (E/c,~p), and the macro-
distribution, and do not seem to be directly transferable scopic 4-velocity by Uµ = (c,~u)γ(u), with ~u the three-
to the present case of the relativistic (Maxwell-Ju¨ttner) dimensionalvelocityofthefluid. Notethatwehaveused
equilibrium distribution, which has fairly different sym- the Anderson-Witting collision operator[20] (rhs of Eq.
metries as compared to the non-relativistic Maxwell- (1)), making our model compatible with the ultrarel-
Boltzmann distribution. Therefore, the extension of the ativistic regime. Hereafter, we will use natural units,
previous LB models has to be considered anew. c = kB = 1, and work in the ultrarelativistic regime,
ξ mc2/k T 1.
In this paper, we develop a new lattice Boltzmann B
≡ ≪
model capable of reproducing the third order momentof According to a standard procedure [13, 15, 16], we
the continuum equilibrium distribution, and still realiz- first expand the Maxwell-Ju¨ttner distribution in the
ableonacubiclattice. Themodelisbasedonasingledis- rest frame, feq = Aexp( p0/T), in an orthogonal
−
tributionfunctionandsatisfiesconservationofbothnum- basis. Since in the ultrarelativistic regime, p0/T =
ber of particles and momentum-energy equations. The p~2/T2+m2/T2 p/T, being p = p~2, we can write
≃
modelis basedonthe singlerelaxationtime collisionop- pthe equilibrium distribution in sphericpal coordinates,
eratorproposedby Andersonand Witting [14, 20]which
is more appropriatefor the ultra-relativisticregime than d3p ∞ π 2π
ge−p0/T = gpe−p/Tdpsin(θ)dθdφ ,
the Marle model used in the previous works, Thus, the Z p0 Z Z Z
0 0 0
proposed model offers significant improvement on previ- (3)
ousrelativisticlatticeBoltzmannmodelsintworespects: where g is an arbitrary function of momentum. Follow-
(i) It captures the symmetry of the higher-order equilib- ing Romatschke [13], we can expand the distribution in
rium moments sufficiently to reproduce the dissipative each coordinate separately, and subsequently, by using
relativistic hydrodynamics at the level of the Grad ap- a Gauss quadrature, calculate the discrete values of the
proximation to the relativistic Boltzmann equation; (ii) 4-momentumvectors. Thus,thediscreteequilibriumdis-
3
tribution can be written as, Order Polynomial Jk k
0th 1 0
fleq =iX,j,kaijk(Uµ)Pi(θl)Rj(pl)Fk(φl) , (4) 21nsdt (p0−6)pp0√0+−262,,(√ppx02−,4√p)py2x,,√p(zp20−4)py 1,5,2,6,3,74
2√3 2√2 2√2
wthheerdeisttrhiebuctoioenffiocinentthseapiojkly(Unoµm)iaalrsePth(eθ )pro(jepc)tFion(sφo)f, (p02−√42)pz, p−xpp0z2,+4pp√xy22p+z2,pyp2xp,y−p024−√36px2 181,, 91,2,1013
and the discrete 4-momenta are deinolteRdj byl pkµl l= 3rd 112(p0−6)22√p20−22√,2((p02−√1204)√p05+20)px 14, 15
(qpule,nptllcyo,s(tφhle)sdiins(cθrle)t,eplfsoirnm(φlo)fsitnh(eθl)B,oplltczoms(aθnln)).eqCuoantisoen- −214(p0−6)(cid:0)p02−3px2(cid:1), 5px32−4√3p502px 16, 17
takes the form, ((p0−104)√p05+20)py, (p0−46√)3pxpy, px4p√y3pz 18, 19, 20
(p0 6)(p02 px2 2py2) px( p02+px2+2py2)
− − − , − 21, 22
p Uµδt − 8√3 8√3
fl(xµ+pµl/p0lδt,t+δt)−fl(xµ,t)=− lµτp0l (fl−fleq)(5). (ppy0(−−63)pp0x2p2+z4√3,p2x2+pz4(ppy022)−,5p((xp20)−,10(4p)√p005−+62)p0y)ppzz 252,32,62,427
4√3 − 8√30 4√3
However,notethat,inthestreamingprocessontheright- (p0−6)pypz, pz(−p02+px2+4py2) 28, 29
hand-side of Eq.(5), the distribution moves at velocity 4√3 24√2
pµ/p0, which implies that the information travels (in a
l l TABLEI:PolynomialsJk thatareorthonormalontheweight
single time step) from each cell center to a position that function w(p0) in Cartesian coordinates (x,y,z).
belongstothesurfaceofasphereofradiuscδt=1. Fur-
thermore, to represent correctly the third order moment
of the equilibrium distribution, Note that in this Table, the 4-momentum has the no-
tation pµ = (p0,px,py,pz). Since these polynomials are
Pαβλ = feqpαpβpλ , (6) orthonormal,there areno normalizationfactors,andthe
l l l l
Xl Maxwell-Ju¨ttner distribution can be approximated, up
to third order in the momentum space, by an expansion
the number of points needed on the surface of the unit
as simple as
sphere exceeds 6 and 12, which correspond to the first
neighborsforacubicandhexagonalclosedpacked(HCP) 29
feq w(p0)a (T,Uµ)J (pµ) , (7)
lattices, respectively. This implies that, in general, the k k
≃
4-vectors pµ/p0 cannot be embedded into a regular lat- Xk=0
tice, and therefore, an interpolation algorithm has to be where the projections a are calculated by,
k
used. Bydoingthis,wearelosingoneofthemostimpor-
d3p
tant features of lattice Boltzmann models, which is the a = feqJ (pµ) . (8)
k Z k p0
exact streaming. Thus, within this spherical coordinate
representation, the streaming process cannot take place Since the Anderson-Witting model is only compatible
on a regular lattice. withtheLandau-Lifshitzdecomposition[14,20],wemust
calculate the energy density of the fluid by solving the
eigenvalue problem,
B. Moment projection of the equilibrium
TαβU =ǫUα , (9)
β
In this work, we shall use a different approach to the ǫ being the energy density of the fluid, and
quadrature representation. We first expand the equilib-
d3p
rium distribution at rest, w(p0) = feq = Aexp( p0) Tαβ = fpαpβ , (10)
− Z p0
by using Cartesian coordinates, unlike the spherical co-
ordinate system used in Ref. [13], and choose the 4- the momentum-energy tensor. For the particle density,
momentum vectors such that they belong to the lattice we use the relation,
(from now on and without loss of generality, we will use
d3p
the notation p0/T p0). This procedure also avoids n=U fpα , (11)
→ αZ p0
extra terms in the product, P (θ ) (p )F (φ ) for the
i l j l k l
R
sphericalcase,whicharenotnecessaryifweonlyneedto and, by using the equation of state, ǫ = 3nT, we can
recover correctly the first three moments of the equilib- calculate the temperature of the fluid.
rium distribution. This considerably simplifies the dis-
crete equilibrium distribution.
By performing a Gram-Schmidt procedure with the C. Discrete-velocity representation of the
weight w(p0), we construct a set of orthonormalpolyno- quadratures
mials. Theorthonormalpolynomialsincartesiancoordi-
nates up to third order, herefrom denoted by J , where NotethattheabovederivationusingCartesiancoordi-
k
the index k runs from 0 to 29, are shown in Table I. nates still refers to the continuous 4-momenta. In order
4
todiscretizethe abovemomentprojectionofthe equilib-
rium distribution, we must choose a set of 4-momentum
vectors that satisfies the very same orthonormality con-
ditions, namely:
d3p
w(p0)J (pµ)J (pµ) = w J (pµ)J (pµ)=δ ,
Z l k p0 i l i k i lk
Xi
(12)
while, at the same time, pµ/p0 corresponds to lattice
points. Here, we choose to work with a cubic lattice, al-
thoughtheproceduredescribedherealsoappliestoother
ones, e.g. HCP lattice.
Since,duetoitsnature,pµ/p0 leadstovelocityvectors
which belong to a sphere of radius c in the space com-
ponents, using the procedure in Ref. [13] will generally
resultinoff-sitelatticepoints. Forthisreason,weoptfor
another quadrature based on this orthonormality condi-
tion, and impose that the distribution function at rest
frame shouldsatisfy the moments ofthe equilibrium dis-
tribution, up to 6-th order. This is made to ensure that
the 5-th order moment of the equilibrium distribution is FIG. 1: Directions of the velocity vectors ϑ~i to recover up to
recovered(at leastatverylow fluid velocities),which, in
the third order moment of the Maxwell-Ju¨ttner distribution.
thecontextoftheGradtheoryfortheAnderson-Witting The radius of the sphere is R = √41. The points represent
model[14],isarequirementforthecorrectcalculationof lattice sites belonging to thespheresurface.
thetransportcoefficients,namelytheshearandbulkvis-
cosities and thermal conductivity. The condition for the
6-thordermoment,istochoosefromthemultiplelattice However, since this newly acquired dependence remains
solutions,the one that presents the highestsymmetry to local, we shall be able to find a discrete-velocity quadra-
model the Maxwell-Ju¨ttner distribution. In order to use ture which also allows for a lattice Boltzmann-type dis-
general features of classical lattice Boltzmann models, cretization in time and space without any interpolation.
likebounce-backboundaryconditionstoimpose zerove- Indeed, in a cubic cell of length δx = 1 there are only 6
locityonsolidwalls,wewillalsorequirethattheweights neighbors, which are not sufficient to satisfy the orthog-
w correspondingtothediscrete4-momentumvectorspk onality conditions and the third order moment of the
i i
have the same values as the ones corresponding to pk equilibrium distribution. However, by multiplying this
− i
(latin indices run over spatial components). equation by a constant R at both sides, and perform-
In order to generate on-site lattice points, let us first ing a time transformation (dilatation), δt Rδt′ and
analyse the relativistic Boltzmann equation, which can τ Rτ′, we obtain →
→
be written as,
p Uµ
p0∂ f +pa∂ f = pµUµ(f feq) , (13) ∂t′f +ϑa∂af =− τµ′p0 (f −feq) , (16)
t a
− τ −
wherewehavedefinedϑa =Rva. Duetothistransforma-
and in the ultrarelativistic regime, tion,the4-momentumvectorsarereconstructedthrough
the relation
p Uµ
p0(∂ f +va∂ f)= µ (f feq) . (14)
t a − τ − pµ =p0(1,ϑ~/R) , (17)
whereva arethe componentsofthe microscopicvelocity.
Atthisstage,wecanchoosetheradiusofthespheresuch
These microscopic velocities have the same magnitude
that the lattice points that belong to the surface of the
but, in general, different directions. Dividing both sides
sphere and the cubic lattice exhibit enough symmetries
of Eq.(16) by p0, we obtain
to satisfy both conditions. This is equivalent to solving
p Uµ the Diophantine equation,
∂ f +va∂ f = µ (f feq) . (15)
t a − τp0 − n2 +n2+n2 =R2 , (18)
x y z
In other words, in the ultra-relativistic regime, the rela-
tivisticBoltzmannequationcanbecastintoaformwhere where nx, ny, and nz are integer numbers, being ϑ~ =
the time derivative and the propagation term become (nx,ny,nz). Thus, we can determine the components of
the same as in the non-relativistic case, at the price of the discrete version of the velocities ϑ~ which are needed
an additional dependence on p0 in the relaxation term. forthe streamingtermintheBoltzmannequation,lhsof
5
Eq. (16). However, on the rhs of this equation, and for w 0 , (22f)
i
≥
the calculation of the discrete 4-momentum vectors via
Eq. (17), we also need to know the discrete values of p0. and look for any solution for wi that fulfills the above
The4-vectorpµ isneededtocomputetheorthonormality relations. Shouldnonebefound,werepeattheprocedure
with a different value of R. By performing this iteration
conditions given by Eq. (12) and the moments of the
process, we found that R = √41 is sufficient to recover
equilibrium distribution.
Due to the fact that p0 is the magnitude of the 4- up to the third order moment of the Maxwell-Ju¨ttner
momentum, p0 = √pµp , in 3 + 1-dimensional space- distribution, andup to sixth orderofthis distribution in
µ
the Lorentz rest frame.
time, it is natural to assume that its discrete values can
The corresponding discrete velocity vectors ϑ~ are:
be calculated by using the weight function in spherical m
coordinates, w(p) = 4πAp2exp( p), where the angular ( 6, 2, 1), ( 6, 1, 2), ( 2, 6, 1), ( 1, 6, 2),
± ± ± ± ± ± ± ± ± ± ± ±
− ( 1, 2, 6), ( 2, 1, 6), ( 5,0, 4), ( 5, 4,0),
componentshavebeenintegratedout,andusingthezeros
± ± ± ± ± ± ± ± ± ±
(0, 5, 4), ( 4, 5,0), (0, 4, 5), ( 4,0, 5),
of its respective orthonormal polynomial of fourth order
± ± ± ± ± ± ± ±
( 4, 3, 4), ( 3, 4, 4), and ( 4, 4, 3); with
(this is because we are interested in an expansion up to
t±he v±alue±s for p±0 ± 0±.743, 2.572,±5.73±1, ±and 10.95.
third order, so we need one more order to calculate the l ≃
Consequently, this gives a total of 4-momentum vectors
zeros). This fourth order polynomial is given by:
= 384. However, the last condition in Eq. (22)
N
1 allows some weights to become zero. Therefore, in our
(4)(p)= [120+p( 240+p[120+(p 20)p])] .
R 24√5 − − iteration procedure, we have taken the minimal number
(19) of 4-momentum vectors pµ, by imposing the maximum
i
theTiorrseusmpemcatirvizeew,in,woredfierrsttoficxalRcualantdestohleveditshcereetqeupaµitiaonnds onnulmyb1e2r8ovfewctiortsopbµenzeeerdoe.d Ftoorfutlhfiilsl rtehaesocno,ndthiteiorensairne
i i
Eq. (22). In principle, all the velocity vectors ϑ~ are
m
n2x+n2y+n2z =R2 , (20a) needed, but only some of the combinations with p0 are
l
required. The detailed list of the ϑ~ , p0, and pµ, and
m l i
(4)(p)=0 , (20b) their respective discrete weight functions wi are given in
R the Supplementary Material [21].
to obtainthe solutions for n , n , n , and p. With these In Fig. 1 we report the configuration of the veloc-
x y z
values, we build the discrete 4-vectors ity vectors ϑ~ to achieve the third order moment of the
Maxwell-Ju¨ttner distribution function. The points cor-
pµlm =p0l(1,nx,m/R,ny,m/R,nz,m/R) , (21) respond to lattice nodes of a cubic lattice that, at the
sametime, belong to the surfaceofthe respective sphere
wherel=1,...,4denotesthefourzerosofthepolynomial of radius R = √41. The relatively large number of dis-
(4)(p), andm=0,..., the triplets (nx,ny,nz)m that cretevelocitiesshouldnotcomeasasurprise;inthecase
R M
satisfytheDiophantineequation,assumingthat isthe of non-relativistic lattice Boltzmann, the number of dis-
M
number ofsolutions. Here, forsimplicity, we regroupthe crete velocities alsobecomes high(at least41 forachiev-
pair of indexes lm to i, so that we can label the discrete ingcompleteGalileaninvarianceinthenon-thermalcase
4-momentumsaspµi,wherei=1,...,N withN =4×M. and 125 in the thermal case, see [18, 19]). Note that the
Next, we replace these values into the equations, specified values of p0 play the same role in defining the
quadrature as the reference temperature (energy) in the
w(p0)J (pµ)J (pµ)d3p = N w J (pµ)J (pµ)=δ , non-relativistic case [18, 19].
Z l k p0 i l i k i lk Finally, we can write the discrete version of the equi-
Xi
(22a) librium distribution up to third order,
29
w(p0)pµpνpσpλd3p = N w pµpνpσpλ , (22b) fieq =wi an(T,Uµ)Jn(pµi) , (23)
Z p0 Xi i i i i i nX=0
whichis shownindetails inAppendix B,Eq.(B2). Note
thatthisdistributionfunctionrecoversthefirstthreemo-
w(p0)pµpνpσpλpγd3p = N w pµpνpσpλpγ , (22c) ments of the Maxwell-Ju¨ttner distribution in the ultra-
Z p0 i i i i i i relativistic regime,
Xi
d3p 128
feqpµ = feqpµ =Nµ , (24)
d3p N Z p0 i i
w(p0)pµpνpσpλpγpβ = w pµpνpσpλpγpβ , Xi=1
Z p0 i i i i i i i
Xi
(22d) d3p 128
feqpµpν = feqpµpν =Tµν , (25)
wi =wj (if pki =−pkj) , (22e) Z p0 Xi=1 i i i
6
d3p 128 and by using a Maxwellian iteration method [14],
feqpµpνpλ = feqpµpνpλ =Pµνλ , (26)
Z p0 i i i
Xi=1 U Pµνβ U Pµνβ = τ∂ Pµνβ . (35)
µ − µ E − µ E
where
Note that we need at least the third order moment of
Nν =nUν , (27) the equilibrium distribution, PEµνβ, to compute the dis-
sipation coefficients (namely, bulk and shear viscosities
and heat conductivity). This requirement is fulfilled in
Tνµ = nTηνµ+4nTUνUµ , (28) our discrete and continuum expansions of the equilib-
−
rium distribution via Eqs. (24), (25), (26). However, to
being the number of particles 4-flow and the energy- recoverfulldissipation,wewouldalsoneedtorecoverthe
momentum tensor, respectively, and third momentof the non-equilibriumdistribution, which
accordingtothe14momentsGrad’stheory,canbewrit-
Pνµλ = 4nT2(ηνµUλ+ηνλUµ+ηµλUν) ten as,
−
+24nT2UνUµUλ ,
Pµνβ =Pµνβ +b Pµνβα+d Pµνβαλ , (36)
(29) E α E αλ E
with n = 2T3. However, the extension to the case of
where b and d are coefficients that carry the infor-
massive particles is straightforward,by changing the co- α αλ
mation on the transport coefficients [14]. Note that we
efficients, a , in Eqs. (8) and (23).
n need to recover terms up to the fifth order of the equi-
librium distribution. In principle, this could be done by
the procedure described on this paper, but the resulting
D. Discrete relativistic Boltzmann equation
value for R could be unpractically large. Nevertheless,
at low velocities, Uµ (1,0,0,0), the Maxwell-Ju¨ttner
In the model of Anderson-Witting for the collision ∼
distribution can be approximated by the weight func-
operator, the relativistic Boltzmann equation takes the tion w(p0), and in analogy to the discrete case, by w ,
i
form given by Eq. (1),
andthe fourthandfifthorderarerecoveredviaEq.(22).
As a result, at relatively low velocities, we expect the
pµU
pµ∂ f = µ(f feq) . (30) non-equilibrium third order tensor to be also fulfilled.
µ
− τ − Therefore,thetransportcoefficientsforanultrarelativis-
tic gas, i.e. µ=0 for the bulk viscosity, η =(2/3)Pτ for
This collision operator is compatible with the Landau-
the shear viscosity, and λ = (4/5T)Pτ for the thermal
Lifshitz decomposition [14], which implies fulfillment of
conductivity, also apply to our model.
the following relations
To discretize the relativistic Boltzmann equation, we
d3p d3p firstimplement the time transformationdescribedin the
UµNµ =UµZ fpµ p0 =UµNEµ =Z feqpµ p0 , previous section and integrate in time Eq. (30) between
(31a) t′ and t′+δt′. This yields:
pµU
UµTµν =UµZ fpµpνdp30p =UµTEµν =Z feqpµpνdp30p , f(xa+ϑaδt′,t′+δt′)−f(xa,t′)=− τ′p0µ(f−feq)δt′ .
(37)
(31b)
By changing pµ pµ, f f and ϑa ϑa, we obtain
Here, the subscript E denotes the quantities calculated → i → i → i
with the equilibrium distribution. Therefore, upon inte-
pµU
grating Eq. (30) in momentum space, we obtain f (xa+ϑaδt′,t′+δt′) f (xa,t′)= i µ(f feq)δt′ .
i i − i − τ′p0 i− i
i
∂ Nµ =0 , (32) (38)
µ
This relativistic lattice Boltzmann equation presents an
which is the conservation of the number of particles 4- exact streaming at the left hand side, and the collision
flow. By multiplying by pν and integrating, we obtain operatoratthe righthandside looksexactlylikeits con-
the conservation of the momentum energy tensor tinuum version. Therefore, the conservationlaws for the
number of particles density 4-flow, and the momentum-
∂ Tµν =0 . (33) energy tensor, are also fulfilled, as long as they are ob-
µ
tainedbyusingtheLandau-Lifshitzdecomposition. This
In order to calculate the transport coefficients, we need means that, first, we need to calculate the momentum-
the third order moment, so that, upon multiplying energy tensor,
Eq. (30) by pνpβ, we obtain
128
1 Tαβ = f pαpβ , (39)
∂ Pµνβ = (U Pµνβ U Pµνβ) , (34) i i i
µ −τ µ − µ E Xi=1
7
and with this tensor, we solve the eigenvalue problem, 0.7 0.34 BAMPS
RLBD 3rd Order
0.3
TαβU =TαβU =ǫUα , (40) 0.6 RLBD 2nd Order
β E β 0.26 RLB
obtaining the energy density ǫ and the 4-vectors Uα. 0.5 0.22
0.18
Subsequently, the particle density can be calculated by 0.4 --11..22 -1 ----0000....888
c
v/
128 0.3
n=U Nµ =U Nµ = f pµU . (41)
µ E µ i i µ
Xi=1 0.2
0.1
The temperature T is obtained by using the equation of
state for the ultrarelativistic gas, ǫ = 3nT. The trans-
0
port coefficients are the same as in the continuum case, -3 -2 -1 0 1 2 3
with the lattice correction resulting from second order z (fm)
Taylor expansion of the streaming term. All factored in,
1.2
the coefficients take the following expression µ = 0, η = BAMPS
RLBD 3rd Order
(2/3)P(τ′ δt′/2), and λ = (4/5T)P(τ′ δt′/2). Note
1 RLBD 2nd Order
− −
thatrevertingbackthetimetransformation,wecanwrite RLB
thetransportcoefficientsasη =(2/3)P(τ δt/2)/R,and
− 0.8 000...666666
λ=(4/5T)P(τ δt/2)/R.
000...666222
−
Summarizing,the presentmodeldoesnotpresentspu- P0 0.6 0.58
rious dissipation in the number of particle conservation P/ 0.54
equation,incontrastto previousRLB schemes[7, 8,12], 00..55
0.4
and also improves the dissipative terms given by the 000...444666
-1.2 -1 -0.8
multi-relaxationtimescheme[17]. Inaddition,itrealizes
0.2
the expansion of the Maxwell-Ju¨ttner distribution on a
cubic lattice, in contrast to Ref. [13]. We can also con-
0
struct a relativistic lattice Boltzmann model that recov- -3 -2 -1 0 1 2 3
ers only up to second order (momentum-energy tensor), z (fm)
to compare with the third order model and determine
the influence of the third order moment in the expan-
FIG.2: Velocity(top)andpressure(bottom)profilesasfunc-
sion. Details of the second order model can be found in
tionofthez-coordinateforthecaseofashockwaveinquark-
Appendix A. gluon plasma, with η/s=0.1.
III. NUMERICAL VALIDATION
[12]. MRT RLB yielded good agreement with the re-
sults at η/s = 0.1, but presented notable discrepancies
In order to validate our model, we solve the Riemann
for η/s = 0.5. The failure of both RLB and MRT RLB
problem for a quark-gluon plasma and compare the re-
to solvethe Riemannproblemforhighviscousfluids can
sults with BAMPS and two previous relativistic Boltz-
be ascribed to their inability to recover the third order
mannmodels. Thefirstone,proposedby Mendozaetal.
moment of the distribution [12, 17].
[7, 8] and later improved by Hupp et al. [12], which we
Inthissection,wewillstudythecaseofhighη/s 0.1
willdenote simply byRLB,andthe secondone,whichis
≥
in a regime of moderate velocities. We perform the sim-
arecentextensionofthe RLBdevelopedbyLietal. [17]
ulations on a lattice with 1 1 1600 cells, only half of
toincludemulti-relaxationtime,whichwewilldenoteby
× ×
which are represented in our domain owing to symme-
MRT RLB. BAMPS was developed by Xu and Greiner
try condition (the other half is a mirror, in order to use
[10] andapplied to the Riemann problemin quark-gluon
periodic boundary conditions for simplicity). Therefore,
plasma by Bouras et al. [11]. Since BAMPS solves the
our simulation consists of 1 1 800 lattice sites, with
full relativisticBoltzmannequation,we takeits resultas
× ×
δx=0.008fm and δt=√410.008fm/c for RLBD third
areferencetoaccesstheaccuracyofourmodel. However,
order, and δt=0.024fm/c for RLBD second order.
we keep in mind that BAMPS also produces approxi-
mate solutions. The present model is hereafter denoted The initial conditions for the pressure are P0 =
by RLBD (RLB with Dissipation). 5.43GfmeV3 and P1 = 0.339GfmeV3. In numerical units, they
For small ratios η/s, where s is the entropy density, correspond to 1.0 and 0.062, respectively. The initial
RLB and MRT RLB reproduced BAMPS results to a temperature z 0 is T = 200MeV (in numerical units
1
≥
satisfactory degree of accuracy. However, for higher 0.5), and T = 400MeV for z < 0, which corresponds
0
η/s 0.1 and moderately fast fluids, γ 1.3, RLB to 1.0 in numerical units. The entropy density s is cal-
failed≥to reproduce the velocity and pres∼sure profiles culated according to the relation, s = 4n nln(n/neq),
−
8
0.6
0.7 BAMPS BAMPS
RLBD 3rd Order T = 1.0
0.6 RLBD 2nd Order 0.5 T = 2.5
RLB
T = 2.8
MRT RLB
0.5
0.4
0.4
c
v/ v/c 0.3 0.4
0.3
0.3
0.2
0.2 0.2
0.1
0.1 0.1
0
2.4 2.7 3
0
0
-3 -2 -1 0 1 2 3
-3 -2 -1 0 1 2 3
z (fm)
z (fm)
1.2
BAMPS
RLBD 3rd Order FIG.4: Velocityprofileasfunctionofthez-coordinateforthe
1 RLBD 2nd Order
caseofashockwaveinquark-gluonplasma,withη/s=0.5,by
RLB
increasingthereferencenumericaltemperatureinthelattice,
MRT RLB
0.8 leading to a smaller relaxation time τ.
0
P 0.6
P/
that there is again an improvement, including the at-
0.4 tainment of the right value of the maximum velocity (at
z 1.5fm). Inthe pressureprofile,RLBDgets closerto
∼
0.2 BAMPSthan MRT RLB in the regionofthe discontinu-
ity in the initial condition (z 0).
∼
0 Note that there is a staircase shape in the results of
-3 -2 -1 0 1 2 3
RLBD for η/s = 0.5 in Figs. 3. This is due to the large
z (fm)
values taken by the single relaxation time in order to
achieve such shear viscosity-entropy density ratios, τ
FIG.3: Velocity(top)andpressure(bottom)profilesasfunc- 20 40 (in numerical units), which is beyond the hydro∼-
tionofthez-coordinateforthecaseofashockwaveinquark- dyn−amic approximation and therefore higher order mo-
gluon plasma, with η/s=0.5.
ments(fourthandhigherorders)ofthedistributionfunc-
tionwouldberequired,whichisnotfulfilledinourRLBD
model. In order to prove this statement, we have per-
where neq is the density calculated with the equilibrium
formedseparatesimulations,seeFig.4,whereweobserve
distribution, neq = d T3/π2, with d = 16 being the
G G thatbyincreasingthevalueofthe referencetemperature
degeneracy of the gluons.
of the lattice (typically set at T = 1), so as to achieve
The velocity and pressure profiles at t = 3.2fm with
c the same shear viscosity, η = (2/3)nT(τ 1/2)/R, the
viscosity-entropy density ratios of η/s = 0.1, are shown −
valueofτ decreasesandthestaircasedisappears. Inpar-
in Fig. 2. In this figure, we compare the results with
ticular, for T 2.5, the results get closer to the ones
BAMPSandRLB,wherewecanseethatRLBpresentsa ≥
with BAMPS, and become independent of the reference
discontinuityatz =0,whilebothsecondorderandthird
temperature. Unfortunately, due the discretization pro-
orderRLBDgetclosertotheBAMPSsolution. Sincethe
cedure used to develop this model, whenever the refer-
only difference between second and third order RLBD is
ence temperature T > 4 the model becomes unstable,
the third order moment of the distribution, we conclude
mostlylikelybecausethe expandedequilibriumdistribu-
thatatrelativelylowη/s,thethirdorderdoesnotplaya
tion function takes negative values.
crucialroleneither inthe conservativedynamics nordis-
sipative dynamics of the system. However, note that at
z 3fm,the third ordermodel provides anoutstanding
∼ IV. CONCLUSIONS AND DISCUSSIONS
fit of the numerical results by BAMPS.
Onthe other hand, by increasingthe ratio η/s, we see
from Fig. 3 that, while RLB gets worse and the second We have introduced a new relativistic lattice Boltz-
order RLBD fixes the discrepancy only in part, the 3rd mann model with improved dissipation, as compared to
order RLBD improves significantly the accuracy of the RLB and MRT RLB. To this purpose, we have per-
velocity and pressure profiles. formed an expansion of the Maxwell-Ju¨ttner distribu-
In Fig. 3, we also compare the results obtained with tion onto an orthonormal basis of polynomials in the
MRT RLB and BAMPS, for η/s = 0.5. Here, we see 4-momentum space. In addition, in order to make the
9
model compatible with a regular cubic lattice, we have
performedtheexpansionincartesiancoordinatesandap-
plied a time transformation, such that particles travel
justthe distancenecessarytoreachlatticenodes,always
at the speed of light. The time transformationgenerates
a sphere of radius R which intersects the cubic lattice,
the intersection points being lattice nodes by construc-
tion. Inaddition,wehavereproduceduptosecondorder
moment of the equilibrium distribution, and up to third
order moment, finding R = 3 and R = √41 for second
and third order moment compatibility, respectively.
The discrete energy component of the 4-momentum,
p0, has been calculated by using Gaussian quadrature,
the nodes corresponding to the zeros of the next order
polynomial. With this configuration,we need 90 vectors
for recovering second order and 384 for the third order
moment case. However, only 66 and 128, respectively,
are actually needed to calculate the moments correctly.
In order to validate the model, we have compared our
results with BAMPS, as well as previous RLB models.
We have found that for η/s=0.1, our model accurately FIG. 5: Directions of the velocity vectors ϑ~i to recover up
describes the Riemann problem in quark-gluon plasma, tothesecondordermomentoftheMaxwell-Ju¨ttner distribu-
includingtheexpansionuptosecondorder. However,for tion,namelythemomentum-energytensor. Theradiusofthe
sphere is R=3. The points represent lattice sites belonging
the case of η/s = 0.5, the second order model, although
to thesurface of thesphere.
better than RLB, is less accurate than both MRT RLB
and the third order model. The third order model yields
better results than the previous RLB, but it develops a
Appendix A: Second order relativistic lattice
staircaseshape asaconsequenceofthe largevalueofthe
Boltzmann model
single relaxation time, which lies beyond the hydrody-
namic regime. We have shownthat the staircase pathol-
To construct the second order lattice Boltzmann
ogy can be tamed by increasing the reference tempera-
model,weusethe proceduredescribedinthis paper. We
ture in the model. Nevertheless, increasing the reference
have obtained that R = 3 presents enough symmetries
temperaturebeyondT =4hitsagainststabilitylimitsof
to fulfill the conditions in Eqs. (22), and the velocity
the model.
vectors ϑ~ are given by, ( 3,0,0), (0, 3,0), (0,0, 3),
We may envisage that a multi-relaxation time exten- ± ± ±
( 2, 1, 2), ( 1, 2, 2), and ( 2, 2, 1). The val-
sion of the present model would further improve the ac- ± ± ± ± ± ± ± ± ±
uesforthediscretep0comefromthesolutionoftheequa-
curacy of the results. A similar improvement may be
tion,
anticipated by implementing higher order expansions of
the equilibrium distribution. However, since the trans-
1
port coefficients depend on the collision operator, their (3) = p0(p0 6)2 2=0 , (A1)
R 12 − −
calculationwithinamulti-relaxationtimemodelbecomes
increasinglyinvolved. Onthe other hand, by performing instead of (4) for the case of the third order expan-
expansionstoincludehigherordermoments,thevalueof sion. ThisRgives the values p0 0.936, 3.305, and
R might become unpractically large, with several ensu- 7.759. The discrete 4-momentulm≃vectors pµ are con-
i
ingdiscretizationissues. Notwithstandingsuchpotential
structed with Eqs. (17), and (21), and they are in total,
difficulties, these extensions are surely worth being ana-
= 3 30 = 90. However, as in the third order ex-
lyzed in depth for the future. N ×
pansion, we have retained the minimal amount, out of
90, that are necessary to recover the second order mo-
ment, by imposing the maximum number of w to be
i
zero. Thisgivesonly664-momentumvectors. Thevalue
of the weight functions for every momentum vector and
Acknowledgments
the relation with the 30 directions are given in the Sup-
plementary Material[21]. In Fig. 5 we reportthe spatial
We acknowledge financial support from the Euro- configuration of the vectors ϑ~ .
i
pean Research Council (ERC) Advanced Grant 319968- ThediscreteversionoftherelativisticBoltzmannequa-
FlowCCS. Work of I.V.K. was supported by the ERC tion, Eq. (38), still applies and the discrete equilibrium
Advanced Grant 291094-ELBM. distribution function is written in detail in Appendix B,
10
Eq. (B1). However, due to the fact that the third order the description of dissipative effects by performing the
momentisnotsatisfied,ananalyticaltheorytocalculate third order expansion and place it on a cubic lattice, we
the transportcoefficientswouldbe verycomplicatedand are not interested in the bulk viscosity and the thermal
goes beyond the scope of this work. Therefore, we have conductivity for this case, and leave this task for future
calculatednumericallyonlytheshearviscosity,bymatch- work.
ing the results for low velocity with the third order mo-
mentmodel. This,inordertocomparetheresultsofboth
expansionswithothermodelsintheliterature. Thisgives
a shear viscosity η (1/7)P(τ δt/2)/R. We could,
2nd ∼ − Appendix B: Equilibrium Distribution Functions
inprinciple,calculatethe thirdordermomentassociated
withtheequilibriumdistributiongivenbyEq.(B1),and,
by applying the Grad method, compute the other trans- The equilibrium distribution function capable to re-
port coefficients. However, this procedure would need cover the first and second order moments of the equilib-
to be performed entirely numerically, since the weights riumdistributionis calculatedby using upto the second
w and 4-momentum vectors pµ are only known numeri- orderpolynomialsin Eq.(7), namely the 14polynomials
i i
cally. Since the mainpurpose ofthis paperis to improve J with k =0,...,13, obtaining
k
feq = nwi p02 T2 2U02 Ux2 Uy2 1 2TU0+1 +2p0(T(T(U0(pxUx+pyUy+pzUz 4U0)+1)
i 4T (cid:20) i (cid:16) (cid:16) − − − (cid:17)− (cid:17) i i i i −
pxUx pyUy pzUz+7U0) 4)+T2 px2 U02+2Ux2+Uy2+1 +2pxUx(pyUy+pzUz 4U0)
− i − i − i − (cid:18) i (cid:16)− (cid:17) i i i −
(B1)
+py2 U02+Ux2+2Uy2+1 +2pyUy(pzUz 4U0)+8U0(U0 pzUz) 2
i (cid:16)− (cid:17) i i − − i − (cid:19)
+2T(5(pxUx+pyUy+pzUz) 8U0)+12 ,
i i i − (cid:21)
For the case of the third order moment expansion, we (k =0,...,29). This leads to the following expressions:
repeat the same procedure, using all the polynomials