Two Dimensional Turbulence in a Massless Fluid with a Relativistic Lattice Boltzmann Model

We investigate a relativistic adaptation of the Lattice Boltzmann Method that reproduces the equations of motion for a turbulent, two-dimensional, massless hydrodynamic system. The classical Lattice Boltzmann Method and its extension to relativistic fluid dynamics is described. The numeric formulation is evaluated using a zero-averaged stirring force introduced into the numerics to induce turbulence, and the flow characteristics produced are compared to properties of a classical turbulent hydrodynamic flow. The model can reasonably be expected to offer quantitative simulations of electron fluid flows in graphene or Kagome lattices.


I. INTRODUCTION
The study of the phenomenon of turbulence in hydrodynamic systems is vital in understanding the dynamics of fundamental systems in the physical world, from atmospheric weather patterns to astronomical systems. Fluid systems with relativistic flows such as the jets that accompany compact astrophysical objects like neutron stars and black holes, or those modeling heavy ion collisions [1] are of interest for their unique hydrodynamic characteristics. The motion of a system of massless particles can be viewed as a relativistic fluid as the dynamics of the individual fermions are described by the Dirac equation. When the particles are considered continuously they form a fermionic fluid governed by relativistic hydrodynamics. It is thought that the flow of the charge carriers in a two-dimensional solid lattice formed by graphene may be considered a massless, relativistic, viscous electronic fluid [2] [3]. The sp 2 hybridized orbit of the electrons in the atoms of the hexagonal lattice form resilient σ bonds with neighboring atoms that create a linear dispersion relation forming so-called Dirac cones in the energy spectrum. Low energy excitations of the particles and the particle holes create massless charge carriers in the lattice whose motion follows Dirac's arXiv:2101.06187v1 [physics.flu-dyn] 15 Jan 2021 equation; Dirac fermions [4]. Near the points of the cones these charge carriers form a massless electronic fluid where the Fermi velocity v F (∼ 10 6 m s ) plays the role of the speed of light. The fluid velocities, bounded from above by v F , are found to be 0.3 percent the speed of light. Despite these non-relativistic velocities, the electron fluid behaves very much like a relativistic fluid since there are no mass-scales in the problem. For this reason, it is expected that the electron velocity probability distribution follows the relativistic Jüttner distribution instead of the classical Maxwell distribution. This niche (non-relativistic velocities yet relativistic velocity distribution) is the focus of the present work. Due to the smaller velocity, the flow in the electronic fluid in graphene is likely pre-turbulent or non-turbulent.
Recently, there has been interest in the lattice structure of a solid known as a Kagome solid that also exhibits a linear dispersion relation in its energy spectrum that forms Dirac cones [5]. This lattice structure appears to support a stronger coupling constant estimated to be 3.2 times that of graphene due to the orbital hybridization. The bulk viscosity caused by phonons and impurities negligible, the dominant coulomb forces drive a larger fluid velocity and a simultaneously smaller shear viscosity. These properties allow the emergence of a turbulent massless fluid that is accessible via experiment. The work by Di Sante et al. [5] uses a holographic approach to determine the shear viscosity of the electronic fluid in a Kagome solid and was able to measure it experimentally.
Due to the non-linearity of the equations of motion, particularly in the turbulent regime, the study of hydrodynamics often requires a numeric approach to find solutions and many different computational methods have been employed. The wide variety of methodologies include mathematical modeling of the Euler, Navier-Stokes, and Boltzmann equations, and modeling with particle-based solvers. Each has its advantages and disadvantages.
The Lattice Boltzmann Method (LBM), born out of a particle-based solver and based on the Boltzmann equation, is one such approach and has had a great deal of success in the modeling of hydrodynamics over the last few decades [6][7] [8]. It is conceptually simple, highly parallelizable, and unconditionally stable in the presence of shocks owing to its reliance on Boltzmann's probability distribution function.
LBM modelers of relativistic hydrodynamic flow are less abundant. A Relativistic Lattice Boltzmann Method (RLBM) described by Mendoza, Boghosian, Hermann and Succi [9] uses parameter matching against conservation laws to fix the form of the equilibrium probability distribution function used in the collision operation. Similar schemes for ultra relativistic flows have been proposed by Mohseni, Mendoza, Succi and Herrmann [10], and Mendoza, Karlin, Succi and Herrmann [11]. A full description of a dimension-independent procedure to design a relativistic lattice Boltzmann scheme is presented by Gabbana, Simeoni, Succi and Tripiccione [12]. A variation of the RLBM proposed by Romatschke, Mendoza and Succi (2011) [13] expands the Jüttner distribution as the equilibrium probability distribution function. An extension of this work, applicable when the chemical potential is not constant, is described by Ambrus and Blaga [14]. The present work extends [13] to model turbulent flows by including a forcing term in the relativistic lattice Boltzmann algorithm. With sufficient viscosity in the collision operator and with large enough forcing magnitudes, this term produces Reynolds numbers in the simulation sufficient to create flows in the turbulent regime. The model is validated through numerical comparison against established turbulent hydrodynamic theory. Using this proposed relativistic turbulent hydrodynamic solver we wish to explore the characteristics of turbulent massless fluid flow.

II. LATTICE BOLTZMANN METHOD
Some explanation of the standard/classical Lattice Boltzmann Method may be appropriate. The standard LBM grew out of the Lattice Gas model, which is a particle-based solver that models particles moving and colliding within a discrete lattice framework. The LBM approach replaces the particles in the lattice with the conceptualized 6 + 1 dimensional probability distribution function f (t, x, v), quantifying the probability that a particle will be in a given state with regard to position and velocity at a given time. The dynamics of the probability distribution function are described by the well-known Boltzmann equation.
v and m represent the particle's velocity and mass, F is an external force, and C [f ] is the collision operator, originally described by Boltzmann as a complex six-dimensional integral.
Boltzmann's form of the operator is difficult to model computationally so it is commonly replaced by a simpler Bhatnagar-Gross-Krook (BGK) type collision operator [15], In this operator f eq is the local probability distribution function at equilibrium, and τ R is the time to relaxation to the local state of equilibrium. This simplified collision operator ansatz sacrifices the complex microscopic details of the particle interactions, but faithfully preserves the evolution of the fluid's macroscopic moments by maintaining conservation of mass, momentum and energy.
The macroscopic moments, specifically the local mass density ρ, the fluid velocity u, and the total (kinetic plus internal) energy are found through integration of the probability distribution function over velocity space [16].
These definitions also hold true when the flow is at a local equilibrium state. Then the probability distribution function is the equilibrium probability distribution function; f → f eq . The equilibrium probability distribution function can be determined in terms of the macroscopic moments. Consider a fluid at equilibrium, when its microscopic velocity is isotropic in velocity space and is therefore equivalent to the macroscopic velocity u (v = u), and can be considered to be constant (|v| 2 = const). Following the same lines of Maxwell's derivation of the equilibrium distribution function [17] and using the definitions in (3), it is straightforward to find the form of the equilibrium probability distribution function in terms of the macroscopic moments: This form introduces the square of the local sound speed as a function of temperature T , This definition of the equilibrium distribution function is used within the BGK collision operator to complete the equations for the evolution of the probability distribution function.
The sacrifice of collision details in favor of the macroscopic moments permits the dynamics of the system to be modeled accurately and efficiently based on the moments with an appropriately defined computational model. The Lattice Boltzmann Method is a family of modeling methods employed to model the dynamics of the probability distribution function through the Boltzmann equation projected onto a discrete spatial lattice with a superimposed discrete velocity space. This provides the model with the minimum amount of information required to solve the Boltzmann equation accurately. There is an abundance of literature on the implementation and application of LBM models; see for instance [8] and [18].
In an LBM model, space is discretized into a cubic lattice where all the information in the probability distribution function for a discrete volume is contained within the nodes of the lattice, separated from each other by a lattice spacing in one, two, or three dimensions.
The velocity space is also discrete at every spatial lattice node using a quadrature analysis method that identifies only the discrete speeds necessary to preserve the moments. In that way, the LBM lattice can be considered a velocity lattice.
A closer look at how a discrete velocity space is determined using a quadrature technique provides insight into the LBM mechanism. Given a function f (x), expressed as f (x) = ω(x)g(x), where ω(x) is an appropriately chosen weighting function, an integral of the form ∞ −∞ ω(x)g(x)dx can be approximated by a summation of the function g(x) evaluated at a suitable choice of abscissae x i and a suitable choice of weights w i where i = 1...n [19].
is represented using a class of orthogonal polynomials, the weighting function ω(x) is the generating function associated with the polynomial set, and the abscissa x i are the n roots of the polynomials. Furthermore, if the function g(x) is of order N or less, then using n orthogonal polynomials evaluates the integral exactly when N ≤ 2n − 1. Here This implies the straightforward use of an expansion of the equilibrium probability distribution function in terms of Hermite polynomials, and the use of the Gauss-Hermite quadrature integral estimation technique.
In this case the discrete weights are determined by: The Gauss-Hermite quadrature integral relation is generalized to d dimensions as [18] This quadrature technique can be used in a straightforward way to evaluate the velocity space integrals (3) using a Hermite polynomial expansion of the probability distribution function.
ρ =m In these summations f i expands f about a velocity evaluated at the abscissae, and incorporates the weights. A discrete summation to evaluate a continuous integral implies there is a finite number of velocities n that are needed to faithfully reproduce the fluid macroscopic moments (3). That is, the roots of the three-dimensional Hermite polynomials, the abscissae Using this discretization formula, the continuous probability distribution function is then where f i is the probability distribution function at lattice node x j for the discrete speed i, incorporating the weighting factor w i . The principal of dynamic similarity asserts that the fluid equations of motion are the same at all length scales, so it is customary within a LBM formulation to rescale the space and time values to dimensionless quantities with respect to some length R ⊥ and frequency

The Lattice Boltzmann Equation is then:
The rescaled quantities are then easily converted to the scale of the domain being modeled when they are recorded. With both Euclidean space and velocity space discretized, the Lattice Boltzmann Method then simulates the time evolution of a hydrodynamic system with discrete time steps following this process: 1. Initialize the local equilibrium distribution function f eq i at each node on the grid according to defined initial conditions of the macroscopic variables, and set the new local probability distribution function to be equal to the equilibrium probability distribution 2. "Stream" the new distribution function f * i to the neighboring node connected by the discrete speed, accounting for boundary conditions. This technique is advantageous computationally since, though the dynamics of a hydrodynamic system are non-linear (see (1)), the non-linear computations are carried out locally at each node independently and then streamed linearly to the (non-local) neighboring nodes.
This allows each node to be computed independently and simultaneously with all the other nodes, providing a high degree of parallelism. One has only to ensure the simpler streaming step happens serially after the collision calculation for all nodes has been completed.

A. Relativistic Boltzmann
The classical version of the LBM has seen wide adoption for the simulation of a variety of hydrodynamic systems. Many variants have been adopted to improve accuracy or performance or to address other difficulties that can arise in different domains.  [20].
For this equation, and for the remainder of this treatment, the units are natural; c = k B = = 1. ∇ µ is the covariant derivative, and Γ λ µν is the Christoffel symbols which are given by derivatives of the underlying metric tensor g µν . For this treatment we assume the Minkowski metric for a flat space-time configuration. The macroscopic fluid velocity is is the Lorentz factor and v i is the microscopic particle 3-velocity. In this description, Greek indices refer to 4-vector space-time components, and Latin indices are used for 3-vector spatial components. In the relativistic regime energy, E = m 2 + p 2 , is no longer quadratic in the velocity. This is reflected in the relativistic form of the equilibrium probability distribution function f eq J which is the Jüttner distribution.
Here T is the local temperature and Z parameterizes the number of degrees of freedom, which in this work will be taken to be one.
The macroscopic moments in relativistic hydrodynamics are obtained from the energymomentum (energy-stress) tensor T µν which can be obtained by integration over the 4 components of momentum space [13].
In this definition m is the particle's mass and H is the Heaviside step function. (Note the notation difference for the relativistic probability distribution function, x, p).) Once obtained, the symmetric energy-momentum tensor contains the energy density (T 00 ), the momentum density (T 0i ), and the momentum flux density (T ik ).
Using the Jüttner form of the equilibrium probability distribution function, integration of (11) at equilibrium produces the equilibrium energy-momentum tensor T µν eq with the familiar relation where both energy density and pressure, and P , are functions of temperature. Then the viscosity coefficient, represented in the collision term, is η = τ R +P 5 , with relaxation time where s is the entropy density [21]. We apply the canonical, so-called Landau-Lifshitz condition, where u µ T µν ≡ u ν , to provide a relationship between the energy density and the fluid velocity in the rest frame of reference. Because the particles modeled are considered to be massless, the equation of state used is simply that of an ideal gas, P = /3.

B. Boltzmann to Lattice Boltzmann
This continuous formulation of a relativistic hydrodynamic system is projected onto a discrete lattice for computational modeling following a prescription similar to the classical LBM. Space and time x µ are discretized in the same manner, but the momentum p µ must be addressed differently. The relativistic form of the equilibrium distribution function is not quadratic in the exponent, and therefore has a different form than the Hermite polynomial Hence, the discretization of the velocity space is less straight forward. To address this, the equilibrium probability distribution function is reexpressed as an expansion about powers of |p|u 0 T 0 θ .
Here the momentum is expressed as unit vectors, v = p |p| , and θ is the scaled temperature with respect to a reference temperature T 0 , θ = T T 0 . In this form we can express the relativistic probability distribution function in terms of a (different) set of orthogonal polynomials as: where p 0 = |p|. The chosen polynomials P wherep = p 0 /T 0 is the scaled momentum magnitude. The parameter α is chosen to be 3, reflective of the number of dimensions and conveniently mapping the coefficients to elements of the energy momentum tensor. The angular portion is made discrete using a similar process. A more complete description of the derivation of the expansion coefficients and the subsequent discretization of the momentum space is described in the work [13].
Having determined the set of orthogonal polynomials to be used to expand the probability distribution function, a discrete set of momenta can now be determined from their roots using the same quadrature technique described in (II). As the polynomials are orthogonal with respect to the solid angle, the resulting discrete three dimensional momentum lattice is spherical, and is not space filling in general (off-lattice). This is overcome computationally using a linear interpolation scheme. An example of a spherical lattice node is shown in figure   2. It remains to obtain the macroscopic moments of the fluid system from the probability distribution function by evaluating the energy-momentum tensor. Using equation (17) and (13) we find: The fluid velocity u µ is determined by the eigenvalue problem set up by the Laudau-Lifshitz condition, u µ T ν µ = u ν , with the energy density as the eigenvalue . The temperature is determined from the relation = 3T 4 π 2 , and pressure P is found from the ideal equation of state, P = 3 .

C. Turbulence Inducing Driving Term
The relativistic modifications to the standard LBM having been determined, the RLBM numerical algorithm follows the standard process of (a The force is incorporated in the RLBM algorithm in the forcing term of the force-included Relativistic Boltzmann Equation.
As with the probability distribution function and the equilibrium probability distribution function, the driving term must be expressed in terms of the orthogonal polynomials The orthogonality relations again determine the expansion coefficients, and are given in (17).
It is straightforward to show the projected form of the driving term evaluated in terms of the probability distribution function and the fluid velocity as (see appendix B): where the integration, dχ = ∞ 0 dp dΩ 4π , can be evaluated using the quadrature technique described. The integrals in (21) are then expressed as expansions with the expansion coefficients in (17) up to the third order. The driving term is calculated and applied in the same step of the LBM process as the collision calculation, and is added to the new probability distribution function, f * = f + ∆tC + F µ ∇ µ f . With this modification in place the RLBM is equipped to produce turbulence and to model a turbulent relativistic hydrodynamic system.

IV. CHARACTERISTICS OF 2D TURBULENCE
When the velocity of a fluid flow reaches a critical speed relative to the viscosity of the fluid, i.e. it attains a critical Reynolds number (Re = uL ν where u is the fluid velocity, L is the characteristic linear dimension and ν is the kinematic viscosity), in the regime of low Mach number the flow will tend to be turbulent. An interesting and widely studied feature of turbulent flow is its tendency to create vortices at different scales within the fluid, expressed by ω = ∇ × v. Due to the non-linearity of the equations of motion obeyed by the flow, the energy forming the vortices at an initial Fourier scale tends to migrate to larger and smaller scales, forming more vortices at the different scales. This is a well established phenomenon known as the cascade of the power spectrum. The power spectrum E(k) is defined as the energy density contained at a given wave number magnitude (denoted simply as k), or the contribution to the energy density by a wave number magnitude [22].
is the energy density in spectral space defined as a function of k: E (k) ≡ 1 2 |ṽ(k)| 2 , where the brackets represent an average over the angular direction. The contribution of the total energy density at all scales determines the energy of the system.
When energy is injected into the system at smaller scales (larger wave numbers), Kolmogorov [23] reasoned that the energy would cascade to larger scales (smaller wave numbers) at a rate proportional to the wave number at the −5/3rd power. This is the so-called −5/3rd power law of the inverse cascade of the power spectrum [23]. Specifically, Kolmogorov calculated through dimensional reasoning that the power spectrum would decay as where is the energy dissipation rate. The form of this power law was derived explicitly by Kraichnan [24], and has since been observed in simulated and experimental data [25].
When energy is injected into the system externally at and/or above a driving scale k f , an inverse cascade of the power spectrum is observed in the inertial range between k f and the largest scales available to the system, 0 < k < k f , and the slope of the power spectrum in that range will be −5/3.
It is noteworthy that a 2006 work by R. K. Scott adjusts the power law to be proportional to the wave number at power −2 ( [26]) in a flow with a non-zero rate of enstrophy transfer.
Then the power spectrum decays as This spectral slope was observed in numerical studies by Westernacher-Schneider [27] [28].
However, because the forcing scheme described in the present work prevents an enstrophy transfer, Kolmogrov's slope of −5/3 is used as the theoretical comparison.

V. MODEL, METHOD AND RESULTS
We consider a two-dimensional isotropic hydrodynamic flow of massless particles under a random, zero-averaged stirring force with infinite boundaries. The RLBM numerical scheme models the system using a spherical D3Q27 momenta configuration connecting N 2 lattice nodes, each containing 27 discrete momenta. The two-dimensional system is represented as the ratio of viscosity to entropy η s , is set to 0.01 in the collision term. Simulations begin initially at rest, and turbulence is introduced into the system as a response to a stirring force applied in Fourier space in zero-averaged random (spectral) directions and at a defined range of wave numbers. A filter in Fourier space restricts the driving wave number range from a chosen wave number magnitude k f to the largest wave number magnitude permitted by the discrete Fourier lattice. Therefore, the driving range is defined as k f ≤ k < k max , leaving the inertial range as 0 < k < k f . Most tests were performed with k f set to 6, and the maximum wave number k max varying with the lattice spacing. At each node in the spectral lattice within the driving range, a driving term is applied at a random magnitude up to a maximum |F |, and in a random spectral direction within the plane so that it averages to zero throughout the lattice. The maximum forcing magnitude is chosen to produce an adequate fluid velocity in the inertial range to induce a turbulent flow. The forcing spectral lattice is then transformed to configuration space and For modest driving force magnitudes, the numerical model produced turbulence inducing fluid velocities of ∼ 0.01 in the inertial range, but yielded no turbulent signal when the driving magnitude was too small. In an inviscid system even small velocities can induce turbulence, but the model simulations suggest the existence of a minimum velocity to show energy propagation across the spectrum in a system with non-zero viscosity. In the RLBM scheme, as viscosity approaches zero, the collision operation becomes infinitely strong (see eq. (2) noting τ R ∝ η), so a non-zero viscosity is required. For a lattice spacing of 0.01, a minimum viscosity of 0.005 was needed to keep the evaluation of the probability distribution function numerically stable, but a smaller viscosity will not affect the stability of models with smaller lattice spacings. The maximum viscosity was found to be 0.18 at the same spectrum causing a "pileup" at the largest discrete scales. The effect is eased by changing the lattice spacing to create more small k modes below the driving range, giving the energy more modes to which it can propagate. Because of the larger errors, the results at that scale of the system are not considered in the conclusions. The simulation configured with N = 128 and α = 0.32 shows a significant pileup at the large scales and lower conformity to the −5/3 power law at small scales. The smaller number of nodes requires a larger lattice spacing to maintain a the same volume, which maximizes the distortion of a discrete lattice and impedes its ability to accurately model a continuous system.
As the inertial range of the spectrum expands the power spectrum forms a better resemblance to the −5/3 power law and the distortion at the large scales is eased, implying a closer conformity at the continuum limit. That is, as the number and range of discrete wave numbers in the inertial range grows the slope of the power spectrum approaches a closer resemblance to what is expected in a continuous system. Figure 6 compares the same lattice configurations as in figure 5, but driven with a forcing configuration leaving the inertial range N α = (1024)(0.16) lattice configuration in figure [7]. The left panel is the energy spectrum at log scale, and the right panel is the energy spectrum of the larger velocity system divided by the smaller one which is rendered as a zero-slope line. No discernible differences in the energy spectrum's slope due to velocity is detected, providing evidence that the numerical model's energy propagation is not dependent on the average velocity, provided it is significant enough induce energy propagation. Further applications of the RLBM in the exploration of turbulent flows should continue to be investigated. In particular, a three dimensional model should be tested, and different turbulence generating schemes should be considered in order to achieve higher velocities.
We thank P. Romatschke for his instruction, helpful discussions and for access to an RLBM source code. We'd also like to thank John Ryan Westernacher-Schneider and Luciano Rezzolla for their helpful comments and insight. This work was supported by the Department The first few orthogonality relations are dΩ 4π

Appendix B: Discrete the Forcing Term
The forcing term is projected onto orthogonal polynomials P (n) up to the second order.
The projection coefficients a (n) are dependent on position and time, and the external force F i is independent of momentum. The divergence of the Boltzmann probability distribution function f is with respect to momentum. The time element of momentum,p = p 0 T 0 , is also the magnitude of the momentum. The first three multi-dimensional orthogonal polynomials are: The coefficients are obtained with the orthogonality relationship, For the zero order coefficient, the orthogonality relationship is, The integrals are regarded as spherical volume integrals in velocity space where Ω is the solid angle, andp is the unit-less magnitude of the momentum representing the radius extending to infinity. We obtain, where the volume integral is represented as V dV = ∞ 0 dpp 2 dΩ 4π . The divergence theorem at an infinite boundary was used, and the divergence of the magnitude of the momentum, given by ∂ i pp = v i T 0 , was also used.
The orthogonal expression for the first order coefficient is, where we evaluated the angular portion of the integral with the relationship dΩ 4π v i v j = 1 3 δ ij , and we used the divergence of the velocity relation ∂ i p v j = 1 p 0 (δ ij − v i v j ). Finally, we used the angular integral dΩ 4π v i v j v l v m = 1 15 (δ il δ jm + δ im δ jl + δ ij δ lm ) to determine the second order projection coefficient.
ij P (2)ij P (2)lm , a (2) ij For this calculation we note the trace of the velocity is unity, containing only angular information. Also, the second order orthogonal polynomial is trace-less, as is its coefficient a (2) ij , so that a (2) ij δ ij = 0. Then the second order projection coefficient is, The forcing term projected onto orthogonal polynomials up to the second order is then,

Appendix C: Minimum Simulation Time
Convergence to a stable energy spectrum for a simulation was achieved before 32 simulated time units. Figure (10) shows the power spectrum for a simulation at t = 2 to 200 time units and indicates a stable spectral slope was realized after t = 32 and remained stable to large time windows. The companion plot depicts the same values divided by the predicted slope providing a higher resolution picture. The plot depicts large errors at the large scales; a result of numerical limitations for small wave numbers in a discrete spectrum. However, the error for the spectrum for the shorter time span engulfs the spectrum for longer time spans, displaying an evolving stability.