Gaussian Quadrature and Lattice Discretization of the Fermi-Dirac Distribution for Graphene
GGaussian Quadrature and Lattice Discretization of theFermi-Dirac Distribution for Graphene
D. Oettinger, M. Mendoza, ∗ and H. J. Herrmann
2, 3 ETH Z¨urich, Department of Physics, CH-8093 Z¨urich, Switzerland ETH Z¨urich, Computational Physics for Engineering Materials,Institute for Building Materials, Schafmattstrasse 6, HIF, CH-8093 Z¨urich (Switzerland) Departamento de F´ısica, Universidade Federal do Cear´a,Campus do Pici, 60455-760 Fortaleza, Cear´a, (Brazil) (Dated: August 13, 2018)We construct a lattice kinetic scheme to study electronic flow in graphene. For this purpose, wefirst derive a basis of orthogonal polynomials, using as weight function the ultrarelativistic Fermi-Dirac distribution at rest. Later, we use these polynomials to expand the respective distribution ina moving frame, for both cases, undoped and doped graphene. In order to discretize the Boltzmannequation and make feasible the numerical implementation, we reduce the number of discrete points inmomentum space to 18 by applying a Gaussian quadrature, finding that the family of representativewave (2+1)-vectors, that satisfies the quadrature, reconstructs a honeycomb lattice. The procedureand discrete model are validated by solving the Riemann problem, finding excellent agreementwith other numerical models. In addition, we have extended the Riemann problem to the case ofdifferent dopings, finding that by increasing the chemical potential, the electronic fluid behaves asif it increases its effective viscosity.
I. INTRODUCTION
Since its discovery [1, 2], graphene has shown a seriesof wonderful electrical and mechanical properties, suchas ultra-high electrical conductivity, ultra-low viscosity,as well as exceptional structural strength, combined withmechanical flexibility and optical transparence. Due tothe special symmetries of the honeycomb lattice, elec-trons in graphene are shown to behave like massless chiral ultrarelativistic quasiparticles, propagating at aFermi speed of about v F ∼ m/s [3, 4]. This placesgraphene as an appropriate laboratory for experimentsinvolving relativistic massless particles confined to a two-dimensional space [5].Electronic gas in graphene can be approached froma hydrodynamic perspective [6–9], behaving as a nearlyperfect fluid reaching viscosities significantly smaller thanthose of superfluid Helium at the lambda-point. Thishas suggested the possibility of observing pre-turbulentregimes, as explicitly pointed out in Ref. [10] and laterconfirmed by numerical simulations [11]. All these char-acteristics in graphene open up the possibility of study-ing several phenomena known from classical fluid dy-namics, e.g. transport through disordered media [12],Kelvin-Helmholtz and Rayleigh B´enard instabilities, justto name a few. However, the study of these phenomenaneeds appropriate numerical tools, which take into ac-count both, the relativistic effects and the Fermi-Diracstatistics.Recently, a solver for relativistic fluid dynamics basedon a minimal form of the relativistic Boltzmann equa-tion, whose dynamics takes place in a fully discrete ∗ [email protected] phase-space lattice and time, known as relativistic lat-tice Boltzmann (RLB), has been proposed by Mendoza etal. [13, 14] (and subsequently revised in Ref. [15] enhanc-ing numerical stability). This model reproduces correctlyshock waves in quark-gluon plasmas, showing excellentagreement with the solution of the full Boltzmann equa-tion obtained by Bouras et al. using BAMPS (BoltzmannApproach Multi-Parton Scattering) [16, 17]. In order toset up a theoretical background for the lattice version ofthe relativistic Boltzmann equation for the Boltzmannstatistics, Romatschke et al. [18] developed a schemefor an ultrarelativistic gas based on the expansion in or-thogonal polynomials of the Maxwell-J¨uttner distribu-tion [19] and, by following a Gauss-type quadrature pro-cedure, the discrete version of the distribution and theweight functions was calculated. This procedure was sim-ilar to the one used for the non-relativistic lattice Boltz-mann model [20–23]. This relativistic model showed verygood agreement with theoretical data, although it wasnot compatible with a lattice, thereby requiring linearinterpolation in the free-streaming step. Another modelbased on a quadrature procedure was developed recentlyin order to make the relativistic lattice Boltzmann modelcompatible with a lattice [24]. However, all these modelsare based on the the Maxwell-J¨uttner distribution, whichis based on the Boltzmann statistics, and therefore, theirapplications to quantum systems is limited.In this work, we construct a family of orthogonal poly-nomials by using the Gram-Schmidt procedure using asweight function the ultrarelativistic Fermi-Dirac distri-bution at rest. By applying a Gauss-type quadrature, wefind that the family of discrete (2+1)-momentum vectors,needed to recover the first three moments of the equilib-rium distribution, are fully compatible with a hexagonallattice, avoiding any type of linear interpolation. Thisresult is very convenient, since the crystal of graphene a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a y shares the same geometry, facilitating the implementa-tion of boundary conditions, allowing for instance hav-ing a good approximation for the electronic transportin nanoribbons with armchair or zigzag edges [25, 26]by implementing the typical bounce-back rule for latticeBoltzmann models.The paper is organized as follows: in Sec. II, we de-scribe in details the expansion of the Fermi-Dirac distri-bution in an orthogonal basis of polynomials, and per-form the Gauss-type quadrature. In this section, we alsoexplain the discretization procedure. In Sec. III, we im-plement the validation of our model by simulating theRiemann problem; and in Sec. IV, we perform additionalsimulations for doped graphene. Finally, in Sec. V, wediscuss the results and future work. II. MODEL DESCRIPTION
The electronic gas in graphene can be considered as agas of massless Dirac quasi-particles obeying the Fermi-Dirac statistics in a two-dimensional space. Thus, we de-fine the single-particle distribution function f ( x µ , p µ ) inphase space, being x µ = ( x , x , x ) and p µ = ( p , p , p )the time-position and energy-momentum coordinates, re-spectively. Here x denotes time, (cid:126)x = ( x , x ) spatial co-ordinates, p the energy, and (cid:126)p = ( p , p ) the momentumof the particles. In the ultrarelativistic regime, we get p µ p µ = 0 (in this paper we use the Einstein notation, i.e.repeated indexes denote summing over such indexes). Inour approach, we assume that the distribution function f evolves according to the relativistic Boltzmann-BGKequation [19], p µ ∂ µ f = − p α U α v F τ ( f − f eq ) , (1)where τ is the relaxation time, and f eq the equilibriumdistribution, which in our case, is the relativistic Fermi-Dirac distribution defined by f eq ( x µ , p µ ) = 1e ( p α U α − µ ) /k B T + 1 , (2)with T the temperature, k B the Boltzmann constant, U µ the macroscopic (2+1)-velocity of the fluid [19, 27],and µ the chemical potential. The relation betweenthe Lorentz-invariant U µ and the classical velocity (cid:126)u =( u , u ) is given by U µ = γ ( v F , u , u ), with v F being theFermi speed and γ = 1 / (cid:112) − (cid:126)u /v F . A. Moment expansion
Here, we perform an expansion of the Fermi-Dirac dis-tribution, Eq. (2), in an orthogonal basis of polynomi-als. In our case, since we are interested in the hydro-dynamic regime, we will truncate the expansion preserv-ing only the polynomials up to second order, although achieving higher orders is also possible by using thesame procedure. In particular, we need to reproducethe first three moments of the equilibrium Fermi-Diracdistribution, namely (cid:104) (cid:105) ( eq ) , (cid:104) p α (cid:105) ( eq ) , and (cid:104) p α p β (cid:105) ( eq ) for α, β = 0 , ,
2. The angular brackets denote expectationvalues using the distribution f via (cid:104) Q (cid:105) = (cid:82) d µ Qf , withd µ = d p/ p (2 π ) , and the subscript ( eq ) indicates thatthe equilibrium distribution f eq is taken instead of f .This method was originally introduced by Grad [28]who expanded the Maxwell-Boltzmann distribution inHermite polynomials, based on the fact that they areorthogonal, using as weight function the Maxwellian dis-tribution at rest. In this spirit, we will derive a newbasis of polynomials that are orthogonal with respect tothe Fermi-Dirac distribution at rest, w ( p ) = 1e p /k B T + 1 . (3)For the following derivations it is useful to choose nat-ural units, c = k B = (cid:126) = 1. In addition, we will consideronly the case for µ = 0, although a general approach isstraightforward. By introducing a reference temperature T , we define θ = T /T , ¯ p = p /T , (cid:126)v = (cid:126)p/ | (cid:126)p | , and using p = | (cid:126)p | , we rewrite the equilibrium distribution as f eq,E ( t, (cid:126)x, ¯ p, (cid:126)v ) = 1e ¯ pγ (1 − (cid:126)v · (cid:126)u ) /θ + 1 , (4)where the subscript E stands for “Exact”. The distribu-tion f eq,E is expanded using tensorial polynomials P ( n ) ,for the angular contribution, and F ( k ) , for the radial de-pendence, such that f eq,E ( t, (cid:126)x, p, (cid:126)v ) = 1e ¯ p + 1 ∞ (cid:88) n,k a ( nk ) i ( t, (cid:126)x ) P ( n ) i ( (cid:126)v ) F ( k ) (¯ p ) . (5)Here, the (2+1)-momentum vectors have been expressedin polar coordinates, p µ = (¯ p, ¯ p cos φ, ¯ p sin φ ) with (cid:126)v =(cos φ, sin φ ) being a unit vector that carries the angulardependence φ , and the index i denotes a family of indices i , ..., i n ∈ { , } whose total number equals the order n of the tensor for the angular dependence, i.e. P ( n ) i and a ( nk ) i are tensors of rank n . Such an ansatz has been usedby Romatschke et al. [18] to expand the Maxwell-J¨uttnerdistribution. Employing the Gram-Schmidt procedure,the radial polynomials F ( k ) are constructed satisfying theorthogonality relation (cid:90) ∞ d¯ p π w (¯ p ) F ( k ) (¯ p ) F ( l ) (¯ p ) = Γ ( k ) F δ kl , (6)and the angular ones by satisfying (cid:90) π d φ π P mi ( φ ) P nj ( φ ) = Γ ( m ) P,ij δ mn . (7)The resulting polynomials and Γ-constants up to secondorder are given in Appendix A. With these polynomials f Φ/π f eq,E f eq f Φ/π f eq,E f eq f Φ/π f eq,E f eq f Φ/π f eq,E f eq FIG. 1. Comparison betweem the expanded Fermi diracdistribution f eq and the full version f eq,E as a function of theangular component φ , for ¯ p = 0 (left) and ¯ p = 3 . θ = 1 . θ = 1 . u = 0 . u = 0 . and taking into account Eq. (5), one can show that up tosecond order in n and k , we get a ( nk ) i = g ( n ) Γ ( k ) F T (cid:90) d¯ p π d φ π f eq,E P ( n ) i F ( k ) , (8)with g (0) = 1, g (1) = 2, and g (2) = 4. The explicitform of a ( nk ) is given in Appendix B. Using Eq. (5), thedefinitions of the polynomials, and their orthogonalityrelations it can be easily shown that the moments up tosecond order can be written in terms of the coefficients a ( nk ) i with n, k ≤ f eq up to secondorder becomes f eq = 1e ¯ p + 1 (cid:88) n =0 2 (cid:88) k =0 a ( nk ) i P ( n ) i F ( k ) . (9)This is sufficient to recover the moments (cid:104) p α (cid:105) ( eq ) = nU α , (10a) (cid:104) p α p β (cid:105) ( eq ) = ( (cid:15) + P ) U α U β − P η αβ , (10b)of the full Fermi-Dirac distribution, Eq. (4). In Eq. (10b),we have introduced the Minkowski metric tensor η αβ , theparticle density n = π T , the pressure P = ζ (3) π nT andthe energy density (cid:15) = 2 P , where ζ denotes the Riemannzeta function, ζ (3) ≈ . f eq and the exact f eq,E , for ¯ p ∼
0, is very poor, in contrast with the case, ¯ p ∼ .
5. However, this isnot surprising, since we are dealing with a gas of ultrarel-ativistic particles which are always moving at the Fermispeed, and therefore none of them has energy ¯ p = 0. Onthe other hand, the matching is reasonable for θ = 1,while being off for θ >
1. Thus, we conclude that θ = 1offers the best approximation, and therefore, we will workwith that value. In addition, we have found that θ can-not be chosen far below unity because f eq can presentnegative values. The fact that θ = 1 implies that thereference temperature T should be equal to the temper-ature of the electronic gas T . B. Momentum space discretization
We now need to discretize the momentum space into afinite number N of discrete momentum vectors, p µq (with q = 0 , ..., N ) such that we can replace integrals in thecontinuum momentum space by sums over a small num-ber of discrete momentum (2+1)-vectors. In order to dothat, we use the Gaussian quadrature [29]. As an exam-ple, for the radial dependence of the expansion, in orderto satisfy (cid:90) ∞ d¯ p π w (¯ p ) F ( k ) (¯ p )¯ p l = N (cid:88) q (cid:48) =0 ω (¯ p ) q (cid:48) w (¯ p q (cid:48) ) w (¯ p q (cid:48) ) F ( k ) (¯ p q (cid:48) )¯ p lq (cid:48) , (11)for k, l ≤
2, we should calculate the discrete ¯ p q (cid:48) andrespective radial weights ω (¯ p ) q (cid:48) . By using the Gaussianquadrature theorem, we found the following values:¯ p = 0 . , ω (¯ p )1 = 0 . p = 2 . , ω (¯ p )2 = 0 . p = 6 . , ω (¯ p )3 = 0 . . Note that in fact, ¯ p is always larger than zero, as ex-pected for ultrarelativistic particles, (see Appendix C fornumerical values with higher precision).On the other hand, by following a similar procedure,we can calculate the N (cid:48) discrete angles φ q (cid:48)(cid:48) and angu-lar weights ω ( φ ) q (cid:48)(cid:48) (with q (cid:48)(cid:48) = 1 , ..., N (cid:48) ), such that, for theangular integrals over P ( n ) ( v i ) l ( v j ) m , one gets (cid:90) π d φ π P ( n ) ( v i ) l ( v j ) m = N (cid:48) (cid:88) q (cid:48)(cid:48) =0 ω ( φ ) q (cid:48)(cid:48) P ( n ) ( v i,q (cid:48)(cid:48) ) l ( v j,q (cid:48)(cid:48) ) m , (12)where v i,q (cid:48)(cid:48) denotes v i ( q (cid:48)(cid:48) ). The above expression is re-quired to be an exact quadrature formula for n ≤ l + m ≤
2. The results for the discrete angles andweights functions are φ q (cid:48)(cid:48) = π + ( q (cid:48)(cid:48) − π and ω ( φ ) q (cid:48)(cid:48) = with N (cid:48) = 6.By combining the radial and angular dependenceof the discrete momentum (2+1)-vectors we get a to-tal of 18 discrete lattice vectors p µ q = p µ ( q (cid:48) ,q (cid:48)(cid:48) ) = e (q',1) e (q',2) e (q',3) e (q',4) e (q',5) e (q',6) y = xx = x FIG. 2. The populations f q are moved between the nodes ofa hexagonal lattice which are linked by the vector (cid:126)e q δt . T (¯ p q (cid:48) , ¯ p q (cid:48) cos φ q (cid:48)(cid:48) , ¯ p q (cid:48) sin φ q (cid:48)(cid:48) ), where we have introducedthe index q = ( q (cid:48) , q (cid:48)(cid:48) ). This lattice cell configuration isshown in Fig. 2, where we can observe that for recov-ering hydrodynamics in graphene, we need a hexagonallattice. This is a very convenient result, since due to thefact that it possesses the same honeycomb lattice sym-metries of graphene, we can reproduce with good accu-racy boundary conditions when modeling nanoribbons orother complex structures.The exact quadrature relations, Eqs. (11) and (12),ensure that the moments up to second order are still rep-resented exactly: (cid:104) p α (cid:105) ( eq ) = (cid:88) q ω q w (¯ p q ) f ( eq ) , q p α q , (13a) (cid:104) p α p β (cid:105) ( eq ) = (cid:88) q ω q w (¯ p q ) f ( eq ) , q p α q p β q . (13b)We have expanded and discretized the Fermi-Diracequilibrium distribution for ultrarelativistic particles.Now, we will proceed to discretize the Boltzmannequation and find the evolution equation for the non-equilibrium distribution. C. Lattice Boltzmann algorithm
With the expanded distribution functions and the dis-cretization of momentum space at hand, we may use thefollowing discrete Boltzmann equation [14, 18, 22], f q ( t + δt, (cid:126)x + (cid:126)e q δt ) − f q ( t, (cid:126)x ) = − p α U α p τ ( f q ( t, (cid:126)x ) − f eq, q ( t, (cid:126)x )) , (14)where we have introduced the notations (cid:126)e q = (cid:126)p q /p ,and f q ( t, (cid:126)x ) = f ( t, (cid:126)x, p q ). Note that (cid:126)e q are unit vec-tors, which means that there are effectively 6 different (cid:126)e q . The discrete Boltzmann equation is now embedded into a lattice, and each time step of δt = 1 correspondsto one execution of the following steps:1. Calculate the equilibrium distributions f eq, q ( t, (cid:126)x )from Eq. (9) using the macroscopic variables n = n ( t, (cid:126)x ), (cid:126)u = (cid:126)u ( t, (cid:126)x ), and T ( t, (cid:126)x ). At t = 0, n ( t =0 , (cid:126)x ), T ( t = 0 , (cid:126)x ), and (cid:126)u ( t = 0 , (cid:126)x ) are imposed asinitial conditions.2. Collision: Introducing the post-collisional distribu-tions f (cid:48) q , calculate f (cid:48) q ( t, (cid:126)x ) = f q ( t, (cid:126)x ) − p α U α p τ ( f q ( t, (cid:126)x ) − f eq, q ( t, (cid:126)x )) . At t = 0, take f q = f eq, q .3. Streaming: Move the f (cid:48) q along (cid:126)e q : f q ( t + 1 , (cid:126)x + (cid:126)e q ) = f (cid:48) q ( t, (cid:126)x )4. Calculate the new macroscopic variables. First wecompute the energy density of the system by solv-ing the eigenvalue problem, (cid:104) p α p β (cid:105) U α = (cid:15)U β , ac-cording to the Landau-Lifshitz decomposition [19].From this, we get (cid:15) and U α . Next, we use the rela-tion n = (cid:104) p α (cid:105) U α = n to obtain the particle density.Here, the average values, (cid:104) p α (cid:105) and (cid:104) p α p β (cid:105) , are sim-ply (cid:104) p α (cid:105) = (cid:88) q ω q w (¯ p q ) f q p α q , (cid:104) p α p β (cid:105) = (cid:88) q ω q w (¯ p q ) f q p α q p β q . The streaming step indicates that if we discretize the realspace based on a hexagonal lattice where the sites arelinked by (cid:126)e q δt , as shown in Fig. 2, the values of f q willbe moved between these sites exactly. This is known as“exact streaming” and crucial for the computational ef-ficiency and accuracy of the lattice Boltzmann methods,because it removes any spurious numerical diffusivity.In summary, we have developed a (2+1)-dimensionalrelativistic lattice Boltzmann scheme with the remark-able feature that it takes into account the Fermi-Diracstatistics, while recovering all the moments up to sec-ond order. The discretization is realized on a hexagonallattice such that exact streaming is achieved. The factthat the quadrature corresponds to a hexagonal latticeallows to represent complex boundaries more preciselyin graphene applications. This will be studied in moredetails in future works.Up to now, we are working with undoped graphene, µ = 0. However, by using the same orthogonal polyno-mials, we can easily integrate the Fermi-Dirac statisticsfor the doped case, obtaining the extended formulation.In this work, we will use µ = 0, in order to compare theresults with previous models in the literature that usethe Maxwell-J¨uttner distribution, since transport theoryshows that in the case of undoped Fermi-Dirac statis-tics, the transport coefficients, namely shear viscosityand thermal conductivity, have the same expresions thanfor the Boltzmann statistics [19]. Therefore, the shearviscosity takes the value of η = (3 / P ( τ − δt/
2) [30].Later, we will use the doped case to study the Riemannproblem, which to best of our knowledge has never beenstudied before. However, it is present when, for instance,laser beams are pointed to the graphene sheet in orderto measure transport coefficients [31].
III. VALIDATION: RIEMANN PROBLEM
In order to validate our model, we solve the Riemannproblem for the ultrarelativistic Fermi-Dirac gas. TheRiemann problem is a standard test for both, relativisticand non-relativistic hydrodynamics numerical schemes,because it involves the evolution of two states of the fluidinitially separated by a discontinuity. In our case, we setup an effectively one-dimensional system of L x × L y =3000 × x and y components. Initially, there are two regionswith particle densities, n = 1 (3 L x / > x > L x / n = 0 .
41 ( x ≤ L x / x ≥ L x /
4) creating arectangular plateau of non-zero particle density in thecenter of the simulation zone. Here we consider and ini-tial constant temperature, T = 1. The initial velocityis set to zero and the value of the relaxation time τ iscalculated for two different values of ξ = η/ ( P δt ), with P = ζ (3) π n T . The evolution of the system is displayedin Fig. 3 after 470 time steps, showing the generatedshock wave. We have only plotted the region x > L x / IV. RIEMANN PROBLEM WITH µ (cid:54) = 0 Let us now consider the case when the chemical poten-tial µ is different from zero. For this purpose, we followthe same procedure described before but this time, wekeep µ (cid:54) = 0. The development is straightforward, andtherefore does not deserve a full explanation. The poly-nomials are the same as described in Appendix A, andthe coefficients a ( nk ) i are calculated by using Eq. (8).The hydrodynamic approach of electrons in grapheneworks for low doping, µ/k B T (cid:28) µ/k B T up to third order, neglecting errors of theorder of ( µ/k B T ) . We perform additional simulations ofthe Riemann problem with the same parameters as be-fore, but now, varying the chemical potential. As we canobserve from Fig. 4, increasing the chemical potentialtends to increase also the effective viscosity of the sys-tem, smoothing the profiles of the velocity, pressure and δ x) n / n Exp., ξ = 4.5Exp., ξ = 14.1 ξ = 4.5 ξ = 14.1 δ x) P / P Exp., ξ = 4.5Exp., ξ = 14.1 ξ = 4.5 ξ = 14.1 δ x) U x / c Exp., ξ = 4.5 Exp., ξ = 14.1 ξ = 4.5 ξ = 14.1 FIG. 3. Density, pressure and velocity profile for the solutionof the Riemann problem. Here ξ = η/ ( P δt ) is a dimension-less number. The expected results were calculated using themodel in Ref. [24]. density. This result is very interesting because it sug-gest that, in fact, impurities with soft potentials (small µ/k B T ) in graphene samples can be treated as local mod-ifications in the effective viscosity of the electronic fluid.In other words, this result suggests a promissing way toinclude impurities in the hydrodynamic approach of elec-trons in graphene. Note that in this figure, there is anoise in the profile of the particle density. This numer-ical instability remains with the same amplitude and isalways located at the boundary when n = n , and there-fore, it does not destroy the stability of the simulation. Itcan be due to the relevance of higher order terms whichare not recovered by our expansion. δ x) n / n µ /T = 0 µ /T = 0.1 µ /T = 0.5 δ x) P / P µ /T = 0 µ /T = 0.1 µ /T = 0.5 δ x) U x / c µ /T = 0 µ /T = 0.1 µ /T = 0.5 FIG. 4. Density, pressure and velocity profile of the solutionof the Riemann problem, for different values of the chemicalpotential µ . V. CONCLUSIONS
We have derived a new family of orthogonal polynomi-als using as weight function the Fermi-Dirac distributionfor ultrarelativistic particles in two dimensions. By ap-plying the Gaussian quadrature we have calculated theset of representative momentum (2+1)-vectors, which al-lows us to replace the integrals over the continuum mo-mentum space by sums over such vectors. As a veryinteresting result, we have found that those vectors pos-sess the same symmetries than the honeycomb latticeof carbon atoms in graphene, making possible the accu-rate implementation of complex boundary conditions infuture applications, such as point defects and nanorib-bons. The derivation has been performed by imposingthat the expanded distribution should fulfill at least thefirst three moments of the equilibrium distribution, which are needed to recover the appropriate hydrodynamics.However, higher order moments can also be recovered byusing the same procedure in this paper.In addition, we have developed a new lattice kineticscheme to study the dynamics of the electronic flow ingraphene. The model is validated on the Riemann prob-lem, which is one of the most challenging tests in nu-merical hydrodynamics, presenting excellent agreementwith previous models in the literature. By increasing thechemical potential, we have found that the profiles of ve-locity, particle density, and pressure, change similar tothe case when the viscosity is increased, concluding thatincreasing the Fermi energy results in increasing the effec-tive viscosity of the electronic fluid. This result suggeststhat soft impurities in graphene samples can be treatedas local modifications of the viscosity, however, furtherstudies must be performed in order to confirm this state-ment.The fact that we can propagate the information fromone site to another in an exact way, avoiding interpola-tion, removes any kind of spurious numerical diffusivity.Therefore, we expect this model to be appropriated tostudy many problems in electronic transport in graphenein the framework of the hydrodynamic approach, e.g.turbulence and hydrodynamical instabilities in grapheneflow, just to name a few.Extensions of the present model to take into accounthigher order moments of the Fermi-Dirac equilibrium dis-tribution as well as the inclusion of the distribution anddynamics of holes, will be a subject of future research.
ACKNOWLEDGMENTS
We acknowledge financial support from the Euro-pean Research Council (ERC) Advanced Grant 319968-FlowCCS.
Appendix A: Polynomials and Γ -constants In this section, we write explicitly the family of polyno-mials, which are orthogonal using as weighting functionthe Fermi-Dirac distribution at rest, with their respec-tive normalization factors. For the case of the angulardependence, we have P (0) ( (cid:126)v ) = 1 P (1) i ( (cid:126)v ) = v i P (2) ij ( (cid:126)v ) = v i v j − δ ij with normalization factors,Γ (0) P = 1Γ (1) P,ij = 12 δ ij Γ (2) P,ijkl = 18 ( δ il δ jk + δ ik δ jl − δ ij δ kl ) . For the case of the radial dependence, we have thepolynomials F (0) (¯ p ) = 1 F (1) (¯ p ) = ¯ p − c ,c = π
12 log(2) ,F (2) (¯ p ) = ¯ p − c ¯ p − c ,c = − π log(2) − π ζ (3))5( π −
216 log(2) ζ (3)) ,c = 7 π − ζ (3) π −
216 log(2) ζ (3)) , with ζ denoting the Riemann zeta function. The normal-ization factors for these polynomials are:Γ (0) F = log(2)4 π , Γ (1) F = − π
576 log(2) + 3 ζ (3)8 π , Γ (2) F = 1400 π (cid:18) π log(2) − π ζ (3) + 48600 ζ (3) ) π −
216 log(2) ζ (3)+ 2250 ζ (5) (cid:19) . Appendix B: Coefficients for the expansion of f eq and relation to moments The coefficients of the expansion in Eq. (9) are givenby a (00) = θ , a (10) i = 2 θ γ + 1 u i γ ,a (20) ij = σ ij θ γ + 1) (cid:2) γ ( u i u j − δ ij ) γ + δ ij (cid:3) ,a (01) = α (1) θ [ θγ − ,a (11) i = 2 α (1) θγ + 1 [ θ ( γ + 1) − u i γ ,a (21) ij = 4 α (1) θ ( γ + 1) [(2 − δ ij ) θ ( γ + 2) − σ ij ] [ γ ( u i u j − δ ij / δ ij / ,a (02) = α (2) θ (cid:2) β (2 | ( θγ −
1) + β (2 | ((3 θγ − − θ )+ β (2 | ( θ (3 γ − − (cid:3) ,a (12) i = 2 α (2) θγ + 1 (cid:2) β (2 | ( θ ( γ + 1) − β (2 | θ (3 θγ − γ + 1)+ β (2 | (3 θ γ ( γ + 1) − (cid:3) u i γ ,a (22) ij = 4 α (2) θ ( γ + 1) (cid:2) β (2 | ((2 − δ ij ) θ ( γ + 2) − (2 δ ij + 1) σ ij )+ β (2 | (3 θ ( γ + 1) − − δ ij )( θ ( γ + 2) − δ ij σ ij ))+ β (2 | (3 θ ( γ + 1) − σ ij ) (cid:3) [ γ ( u i u j − δ ij ) γ + δ ij ] , where σ ij = ( − δ ,i δ ,j or( σ ij ) = (cid:18) − (cid:19) , (B1) and, α (1) = 12 π log(2)216 log(2) ζ (3) − π ,α (2) = 5 [2250 ζ (5)(216 log(2) ζ (3) − π ) + 210 π ζ (3) − π log(2) − ζ (3) ] − ,β (2 | = − π log(2) ,β (2 | = − π ζ (3) ,β (2 | = 3240 log(2) ζ (3) , which are approximately, α (1) ≈ . α (2) β (2 | ≈− . α (2) β (2 | ≈ − . α (2) β (2 | ≈ . f eq , weexpressed them in terms of the a ( nk ) i using Eqs. (8),(9), and the expressions in Appendix A, e.g. (cid:104) p (cid:105) = T (cid:16) Γ (1) F a (01) + c Γ (0) F a (00) (cid:17) .Note that for the calculation of the coefficients a ( nk ) i we should use of the integration formula (cid:90) ∞ d x x n − z − e a x + 1 = − z − a − n Γ( n ) Li n ( − z ) , which holds for n > a ∈ R , a >
0. Here, Γ( n ) denotesthe gamma function, which becomes Γ( n ) = ( n − n ∈ N . Li n ( z ) is the polylogarithm which can be definedusing a power series: Li n ( z ) = (cid:80) ∞ k =1 z k k n . If we considerthe chemical potential in the Fermi-Dirac distribution tobe zero, we have z = 1 and the relevant values of thepolylogarithm become Li ( −
1) = − log(2) , Li ( −
1) = − π , Li ( −
1) = − ζ (3). On the other hand, for µ (cid:54) = 0,we take z = e µ/T . Appendix C: Results for radial Gaussian quadrature
When the radial Gaussian quadrature is applied, thefollowing values for the discrete ¯ p i are obtained:¯ p = 0 . , ¯ p = 2 . , ¯ p = 6 . , with its respective weight functions ω (¯ p )1 = 0 . ,ω (¯ p )2 = 0 . ,ω (¯ p )3 = 0 . . [1] K. Novoselov, A. Geim, S. Morozov, D. Jiang, M. Kat-snelson, I. Grigorieva, and S. Dubonos, Nature Letters , 197 (2005).[2] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A.Firsov, Science , 666 (2004).[3] D. P. DiVincenzo and E. J. Mele, Physical Review B ,1685 (1984).[4] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi,Rev. Mod. Phys. , 407 (2011).[5] A. K. Geim and A. H. MacDonald, Phys. Today , 35(2007).[6] M. M¨uller, J. Schmalian, and L. Fritz, Physical ReviewLetters , 025301 (2009).[7] M. M¨uller and S. Sachdev, Phys. Rev. B , 115419(2008).[8] L. Fritz, J. Schmalian, M. M¨uller, and S. Sachdev, Phys.Rev. B , 085416 (2008).[9] M. M¨uller, L. Fritz, and S. Sachdev, Phys. Rev. B ,115406 (2008).[10] M. M¨uller, J. Schmalian, and L. Fritz, Phys. Rev. Lett. , 025301 (2009).[11] M. Mendoza, H. J. Herrmann, and S. Succi, Phys. Rev.Lett. , 156601 (2011).[12] M. Mendoza, H. J. Herrmann, and S. Succi, Sci. Rep. ,1052 (2013).[13] M. Mendoza, B. M. Boghosian, H. J. Herrmann, andS. Succi, Phys. Rev. Lett. , 014502 (2010).[14] M. Mendoza, B. M. Boghosian, H. J. Herrmann, andS. Succi, Phys. Rev. D , 105008 (2010).[15] D. Hupp, M. Mendoza, I. Bouras, S. Succi, and H. J.Herrmann, Phys. Rev. D , 125015 (2011). [16] Z. Xu and C. Greiner, Phys. Rev. C , 064901 (2005).[17] I. Bouras, E. Molnar, H. Niemi, Z. Xu, A. El, O. Fochler,C. Greiner, and D. H. Rischke, Phys. Rev. Lett. ,032301 (2009).[18] P. Romatschke, M. Mendoza, and S. Succi, Phys. Rev.C , 034903 (2011).[19] C. Cercignani and G. M. Kremer, The RelativisticBoltzmann Equation: Theory and Applications (Boston;Basel; Berlin: Birkhauser, 2002).[20] X. He and L.-S. Luo, Phys. Rev. E , 6811 (1997).[21] N. S. Martys, X. Shan, and H. Chen, Phys. Rev. E ,6855 (1998).[22] X. He and L.-S. Luo, Physical Review E , R6333(1997).[23] S. Succi, The lattice Boltzmann equation for fluid dynam-ics and beyond (Clarendon Press and Oxford UniversityPress, Oxford and New York, 2001).[24] M. Mendoza, I. Karlin, S. Succi, and H. J. Herrmann,Phys. Rev. D , 065027 (2013).[25] M. Y. Han, B. ¨Ozyilmaz, Y. Zhang, and P. Kim, Phys.Rev. Lett. , 206805 (2007).[26] V. Barone, O. Hod, and G. E. Scuseria, Nano Letters ,2748 (2006).[27] F. J¨uttner, Z. Physik (Zeitschrift f¨ur Physik) , 542(1928).[28] H. Grad, Communications on Pure and Applied Mathe-matics , 331 (1949).[29] P. J. Davis and P. Rabinowitz, Methods of numerical in-tegration , 2nd ed. (Academic Press, Orlando, 1984).[30] M. Mendoza, I. Karlin, S. Succi, and H. J. Herrmann,JSTAT , P02036 (2013).[31] J.-U. Lee, D. Yoon, H. Kim, S. W. Lee, and H. Cheong,Phys. Rev. B83