TURBULENCE IN TWO-DIMENSIONAL RELATIVISTIC HYDRODYNAMIC SYSTEMS WITH A LATTICE
BOLTZMANN MODEL

by

MARK WATSON

B.S., West Texas A & M University, 1993

B.S., West Texas A & M University, 1996

M.S., University of Colorado Colorado Springs, 2019

A dissertation submitted to the Graduate Faculty of the

University of Colorado Colorado Springs

in partial fulfillment of the

requirements for the degree of

Doctor of Philosophy

Department of Physics and Energy Science

2022

This dissertation for the Doctor of Philosophy degree by

Mark Watson

has been approved for the

Department of Physics and Energy Science

by

Paul Romatschke, Chair

Kathrin Spendier

Robert Camley

Michael Calvisi

Mark Hoefer

Date  May 9, 2022

Watson, Mark (Ph.D., Applied Science — Physics)
Turbulence in Two-Dimensional Relativistic Hydrodynamic Systems with a Lattice
Boltzmann Model
Dissertation directed by Professor Paul Romatschke.

ABSTRACT

Using a Lattice Boltzmann hydrodynamic computational modeler to simulate
relativistic fluid systems we explore turbulence in two-dimensional
relativistic flows. We first a give a pedagogical description of the
phenomenon of turbulence and its characteristics in a two-dimensional system.
The classical Lattice Boltzmann Method and its extension to relativistic fluid
systems is then described. The model is tested against a system incorporating
a random stirring force in k-space and then applied to a realistic sample of
graphene.

Part I: We investigate the relativistic adaptation of the Lattice Boltzmann
Method reproducing a turbulent, two-dimensional, massless hydrodynamic system
with a zero-averaged stirring force randomly generated in momentum space. The
numeric formulation is evaluated and the flow characteristics produced are
compared to properties of classical turbulence. The model can reasonably be
expected to offer quantitative simulations of charged fluid flows in two-
dimensional relativistic fluid systems.

Part II: At low Reynolds numbers, the wind flow in the wake of a single wind
turbine is generally not turbulent. However, turbines in wind farms affect
each other’s wakes so that a turbulent flow can arise. An analogue of this
effect for the massless charge carrier flow around obstacles in graphene is
outlined. We use a relativistic hydrodynamic simulation to analyze the flow in
a sample containing impurities. Depending on the density of impurities in the
sample, we indeed find evidence for a potentially turbulent flow and discuss
experimental consequences.

###### Contents

  1. 1 Introduction
  2. 2 Background
    1. 2.1 Hydrodynamic Equations of Motion
    2. 2.2 Turbulence
    3. 2.3 Charged Quasi-Particle Flow
  3. 3 The Lattice Boltzmann Method
    1. 3.1 The Standard/Classical LBM
    2. 3.2 Relativistic Lattice Boltzmann
    3. 3.3 Boltzmann to Lattice Boltzmann
  4. 4 Part I: Reproduction of a 2D turbulent flow with the RLBM
    1. 4.1 The Driving Term
    2. 4.2 Model and Methods
    3. 4.3 Results
    4. 4.4 Conclusions
  5. 5 Part II: Reproduction of turbulence in graphene
    1. 5.1 Model and Methods
    2. 5.2 Results
      1. 5.2.1 Single Impurity
      2. 5.2.2 Multiple Impurities
    3. 5.3 Conclusions
  6. 6 Summary and Conclusions
  7. A Maxwellian Form of the Equilibrium Distribution Function
  8. B Three Dimensional Orthogonal Polynomials
  9. C Discrete the Forcing Term
  10. D Minimum Simulation Time
  11. E Determination of the Equation of State
  12. F Computational Platform

###### List of Figures

  1. 2.1 Current density spectra
  2. 2.3 Charged Quasi-Particle Flow
  3. 2.4 Fermi surface & Dirac cones
  4. 3.1 D3Q19 cubic lattice
  5. 3.2 Spherical lattice
  6. 4.1 Average fluid velocity
  7. 4.2 Determination of a constant velocity
  8. 4.3 Energy spectrum for $k_{f}>6$
  9. 4.4 Energy spectrum for $k_{f}>14$
  10. 4.5 Energy spectrum for varying velocities
  11. 4.6 Energy spectrum for varying volume
  12. 4.7 Energy spectrum lattice spacing ratio
  13. 5.1 Airflow in a wind farm
  14. 5.2 Charged flow about a single impurity
  15. 5.3 Current density from a single impurity
  16. 5.4 Charged flow about two impurities
  17. 5.5 Current density from two impurities
  18. 5.6 Current density from two close impurities
  19. 5.7 Current density from two distant impurities
  20. 5.8 The effect of separation on current density
  21. D.1 Development of Energy spectrum in time
  22. D.2 Determination of the time step

##  Chapter 1 Introduction

The study of turbulence in fluid systems is vital to the understanding of
dynamics of many fundamental systems in the physical world, from atmospheric
weather patterns to astronomical systems, and even the flow of electricity.
Fluid systems with relativistic flows can be found associated with compact
astrophysical objects such as the jets accompanying neutron stars or black
holes. Relativistic fluid flows arise naturally in the effects of a collision
between two black holes in AdS (Anti-de Sitter) space (see [1]). The fluid
created by heavy ion collisions, quark-gluon plasma, is of interest for its
unique hydrodynamic characteristics [2] [3] [4] [5] [6] [7] [8] [9] [10]. Work
in that fluid system has lead to a number of successful predictions at the
RHIC and the LHC. See, for instance [11], [12], [13] and [14]. The work in
[15] successfully predicted the particle emission patterns in 3He + Au
collisions at the RHIC.

Also of interest is the flow of massless charge carriers in a two-dimensional
solid lattice such as that formed by graphene, which is considered to be a
viscous relativistic charged fluid [16] [17]. This system of particles can be
viewed as a relativistic fluid since the dynamics of the individual fermions
are described by the Dirac equation. When the particles are considered as a
continuous volume they form a fermionic fluid governed by relativistic
hydrodynamics.

Graphene is a single atom thick sheet of graphite with interesting electrical
properties. The two-dimensional solid is made of carbon atoms arranged in a
honeycomb lattice, showing a high conductivity through a simple band structure
[18] [19]. In the charge neutral state each unit in the lattice contains
exactly one electron so that half of the energy levels are occupied. In this
state the $sp^{2}$ hybridized orbit of the electrons in the atoms of the
hexagonal lattice form resilient $\sigma$ bonds with neighboring atoms that
create a nearly linear dispersion relation, forming so-called Dirac cones on
the Fermi surface. The relation’s curves only begin to appear at a very high
temperature; $\approx 10^{5}K$, validated experimentally in [20][21]. The band
gap disappears at points on the surface at the tips of the cones called Dirac
points. Low energy excitations of the particles and particle holes create
highly mobile, massless chiral quasi-particles [22] in the lattice that
transport the electric charge through the conical band structure following
Dirac’s equation of motion [23]. Near the Dirac points these charge carriers
form a relativistic fluid where the Fermi velocity $v_{F}$ $(\sim
10^{6}\frac{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.

The effect of the lack of a mass-scale in the kinetic theory can be seen in
the form of the Maxwell-Jüttner distribution (also called the Jüttner
distribution), which is a probability distribution of the velocities of
particles in a hypothetical gas made up of relativistic particles [24][25].
Its form is:

$f(p^{\alpha})=De^{\left(-\frac{U^{\alpha}p_{\alpha}}{k_{B}T}\right)}.$ (1.1)

In this distribution $U^{\alpha}$ is the fluid 4-velocity, $p_{\alpha}$ is the
4-momentum and $D$ is a normalization constant. The classical analog of this
is the Maxwell-Boltzmann distribution function (also called the Maxwell
distribution),

$f(v)=De^{\left(-\frac{mv^{2}}{2k_{B}T}\right)}.$

It can be seen that the Jüttner distribution becomes the Maxwell distribution
[26] at non-relativistic velocities and non-zero mass.

The transport of the quasi-particles demonstrates a high conductivity that
stays above a minimum value, even when carrier concentrations approach zero
[22]. The viscosity to entropy ratio is determined to be very small ($\sim
0.2\frac{\hbar}{k_{B}}$)[27]; smaller than that of superfluid helium ($\sim
0.8\frac{\hbar}{k_{B}}$)[27][28], and in some cases that of quark-gluon plasma
($0.07\leq\frac{\eta}{s}\leq 0.43\frac{\hbar}{k_{B}}$)[29]. The number density
of the quasi-particles is determined using electrical properties such as the
Hall resistance, which normally assumes an integer value in a two-dimensional
solid, but in graphene it is found in multiples of $1/2$, known as the
fractional Quantum Hall Effect [30]. In the presence of a strong magnetic
field the wave functions of the quasi-particles behave like simple quantum
harmonic oscillators that take on discrete energy levels based on the strength
of the magnetic field. This effect, called the Shubnikov-De Haas effect, is
used to determine the effective mass of the charge carrying quasi-particles.

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 [31]. 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 is negligible, so the Coulomb forces dominate and drive
a larger fluid velocity with 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. [31] used a holographic
approach to determine the shear viscosity and Reynolds number of the charged
fluid in a Kagome solid.

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 fluid dynamics (CFD)
methods have been employed. The Navier-Stokes solvers include the Finite
Volume Method, the Finite Element Method, the Finite Difference Method and the
Spectral Element Method. They solve the Euler and Navier-Stokes equations, and
typically perform fairly well. For open systems, the spectral method ranks
high in performance, though the Lattice Boltzmann Method ranks comparable to
it for a closed system due to its high parallelizability. Particle based
solvers track individual particles in the fluid and the forces against them.
They include the Lattice Gas Method, Dissipation Particle Dynamics, Multi-
Particle Collision Dynamics, Direct Simulation Monte Carlo, and Smooth
Particle hydrodynamics. These solvers provide high accuracy, but large numbers
of particles and velocities are a challenge.

The Lattice Boltzmann Method or Lattice Boltzmann Model (LBM) was born out of
particle based solvers and is based on the Boltzmann equation. It is one
computational approach that has had a great deal of success in hydrodynamic
modeling over the last few decades [32][33][34][35][36][37]. 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 in a work from Mendoza,
Boghosian, Herrmann and Succi [38] 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
[39], and Mendoza, Karlin, Succi and Herrmann [40]. A full description of a
dimension-independent procedure to design a relativistic lattice Boltzmann
scheme is presented by Gabbana, Simeoni, Succi and Tripiccione [41]. A
variation of the RLBM proposed by Romatschke, Mendoza and Succi (2011) [42]
expands a Jüttner form of the equilibrium probability distribution function
(as opposed to the Maxwellian form used in classical LBM). An extension of
this work, applicable when the chemical potential is not constant, is
described by Ambrus and Blaga [43]. The present work employs another extension
of the approach described by [42].

In this thesis submission we present the application of a relativistic
hydrodynamic numeric modeler to the study of the turbulent characteristics of
the unusual flow of massless charge carrying quasi-particles within a two-
dimensional solid, using graphene as the principal platform.

In the first section we begin with a brief description of the classical
hydrodynamic equations of motion and their primary moments. In greater detail
we next define and describe the nature of turbulence, how its presence is
identified and the differences of the characteristics of turbulence within a
two-dimensional vs. a three-dimensional flow.

In the subsequent section we present a detailed description of the
hydrodynamic modeler that is applied in this work, the Lattice Boltzmann
Method, including a mathematical justification of its application, and its
strengths and weaknesses. The next two sections describe the application of
the modeler to the problem of turbulence. With the introduction of a forcing
term in the relativistic lattice Boltzmann algorithm, we explore the energy
propagation resulting from turbulence generated by a two-dimensional stirring
force of random magnitude originating in momentum space. With large enough
forcing magnitudes against small viscosity, this force 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. Next, we model the flow of the massless quasi-
particles in a realistic sample of graphene with one or more obstacles in
order to determine if a turbulent signal can be detected. We conclude with a
discussion and general observations of the work as a whole.

This work describes and validates an effective new formulation of the Lattice
Boltzmann Method that can be used to study the effects of turbulence in two-
dimensional relativistic systems. This new methodology is then applied to the
charge carrying quasi-particles in graphene and demonstrates that turbulence
is possible under specific conditions, even at Reynolds numbers ranging in the
transitional phase.

##  Chapter 2 Background

###  2.1 Hydrodynamic Equations of Motion

Hydrodynamics is an effective theory describing the dynamics of a volume of a
differentially small fluid sometimes referred to as a parcel. The objects that
constitute the parcel are presumed to be so small and numerous that they are
considered collectively as a continuous volume. As it moves, the boundaries of
the parcel can deform and it can collide with other parcels, subject to
conservation laws. The dynamics of the objects that make up the parcel are
sacrificed in favor of the local fluid properties; density, pressure,
temperature and flow velocity, which are defined locally within the fluid and
can vary spatially and temporally throughout [44]. These properties are also
subject to some thermodynamic equation of state that helps to complete the
description of the fluid system. An excellent derivation of fluid dynamics as
a purely continuum theory, without regard to a parcel or its constituent
particles is given in [45].

Taking a Newtonian approach, the fluid’s acceleration expressed in terms of
the fluid velocity $\bm{u}$ (a three-vector which are indicated in bold) is

$\frac{d\bm{u}}{dt}=\frac{\partial\boldsymbol{u}}{\partial
t}+\left(\boldsymbol{u}\cdot\nabla\right)\boldsymbol{u}.$

The first term is the local acceleration of the fluid and the second is the
advection term that tracks the acceleration of the fluid parcel as it moves
from site to site.

If we assume that forces internal to the fluid vary only spatially within it,
and are normal to the parcel boundaries, then we can track an internal driving
term as a pressure gradient force, $\bm{F_{int}}=-\frac{\nabla p}{\rho}$, with
units of force per unit mass. It is the difference of the pressure on opposite
borders of the parcel (e.g. $\delta p=p_{left}-p_{right}$) times the area of
the border (e.g. $\delta y\delta z$), per the mass of the parcel
($m=\rho\delta x\delta y\delta z$). The acceleration of the parcel with this
term for internal forces and some other external forcing term $\bm{g}$ (also
in units of force per unit mass) comprise the momentum equation.

$\frac{\partial\boldsymbol{u}}{\partial
t}+\left(\boldsymbol{u}\cdot\nabla\right)\boldsymbol{u}=-\frac{\nabla
p}{\rho}+\boldsymbol{g}.$ (2.1)

This is the _Euler equation_ for non diffusive and incompressible fluids (i.e
$\nabla\cdot\boldsymbol{u}=0$) [46].

The Navier-Stokes equation of fluid motion adds another forcing term, a
dissipative forcing term, employing a viscous stress tensor $\mathscr{T}$.
This tensor is the sum of two spatial derivative expressions of the fluid
velocity scaled separately with shear viscosity $\mu$ and bulk viscosity
$\lambda$; constants determined experimentally. Using index notation with
Einstein’s summation convention, the viscous stress tensor is defined as

$\tau_{ij}=\mu\left(\frac{\partial u_{i}}{\partial x^{j}}+\frac{\partial
u_{j}}{\partial x^{i}}\right)+\lambda\left(\frac{\partial u_{i}}{\partial
x^{j}}\right)\delta_{ij}.$

The elements of this tensor are presumed to be linear and consistent
throughout the fluid, dissipating the fluid based on the different components
of the fluid velocity. They produce a force that is tracked as a divergence
across the boundaries of the parcel per unit mass,
$\frac{1}{\rho}\nabla\cdot\tau$. Incorporating this term into 2.1 gives the
_Navier-Stokes equation_ :

$\frac{\partial\boldsymbol{u}}{\partial
t}+\left(\boldsymbol{u}\cdot\nabla\right)\boldsymbol{u}=-\frac{\nabla
p}{\rho}+\frac{1}{\rho}\nabla\cdot\mathscr{T}+\boldsymbol{g}.$ (2.2)

This equation is not generally solvable analytically. Where possible,
solutions are often sought using computational methods.

For relativistic systems, the dynamics are expressed in space-time components
using the energy-momentum 4-tensor $T^{\mu\nu}$ [47]. The energy-momentum
tensor is the flux of the momentum 4-vector $p^{\mu}$ over the surface of
$x^{\nu}$. It contains as its elements the density and current of both energy
and momentum. The components are defined by

$T^{\mu\nu}=\left(\epsilon+p\right)u^{\mu}u^{\nu}+g^{\mu\nu}p,$ (2.3)

where $u^{\mu}$ is the fluid velocity 4-vector, $\epsilon$ is the energy
density, and $p$ is the pressure. Please note the subtle difference in
notation between the scalar pressure $p$ (italics), the momentum 3-vector
$\bm{p}$ (boldface), and the momentum 4-vector $p^{\mu}$ (italics with Greek a
upper or lower index). The tensor $g^{\mu\nu}$ is the space-time metric tensor
which, for flat space-time in Minkowski space, has diagonal elements $1,-1,-1$
and $-1$ [47]. (Note, 2.3 is expressed in natural units where $c=1$.) With
this the equations of motion for relativistic fluid dynamics is given simply
as

$\frac{\partial T_{\nu}^{\mu}}{\partial x^{\nu}}=0.$ (2.4)

$T^{i0}$ are the spatial componentized momentum densities, $T^{0i}$ are the
spatial components of energy current, and $T^{ij}$ are the spatial momentum
current components (here $i$, and $j$ are the spatial coordinate indexes
$1-3$). In a co-moving frame of reference (where $u^{0}=0$ and $u^{i}=0$)
$T^{00}$ is interpreted physically as the energy density.

Hydrodynamics describes the dynamics of a fluid system relaxing to a thermal
equilibrium. Considering a relatively large volume, and a long enough time
scale, the dynamics are limited to locally conserved quantities such as
energy, momentum and charge that are determined at the local equilibrium
state. The state principal of equilibrium thermodynamics declares that any
local thermodynamic state of the fluid, pressure, density, temperature,
internal energy and entropy, can be related to any other two with an equation
of state. When defined, this equation can provide an additional relation to
help find solutions to the system. A common example is the equation of state
for an ideal gas [48] ;

$p=\rho RT\frac{1}{M},$

where $R$ is the gas constant and $M$ is the molar mass in kilograms per mole.

Ludwig Boltzmann describes the dynamics of a fluid statistically with the
Boltzmann transport equation. For classical systems, it describes the changes
of a probability distribution function parameterized with position and
momentum across time, $f(t,\boldsymbol{x},\boldsymbol{p})$. (Velocity is
sometimes used in the place of momentum.) A complete time differential of $f$
is shown in the Liouville expansion

$\frac{df}{dt}=\frac{\partial f}{\partial t}+\frac{\bm{p}}{m}\cdot\nabla
f+\bm{F}\cdot\nabla^{p}f,$

where $\bm{F}=\frac{d\bm{p}}{dt}$ represents a summation of some external
driving forces present causing a change in the momentum, and
$(\nabla^{p}\cdot)$ is the divergence across momentum space.

The _Boltzmann transport equation_ accounts for the dynamics of $f$ by setting
it equivalent to a collision term $C[f]$ that represents internal forces
caused by interactions between the particles:

$\frac{\partial f}{\partial t}+\frac{\bm{p}}{m}\cdot\nabla
f+\bm{F}\cdot\nabla^{p}f=C[f].$

The form of the collision operator proposed by Boltzmann is a complicated,
multi-volume integral that is computationally challenging even with a
numerical approach. But if we assume collisions are rare the term can
simplified. A popular approximation of this type is the Bhatnagar-Gross-Krook
(BGK) collision operator [32] [36]. It assumes that all collisions can be
broken down to be between just two particles at a time, and that there are so
many collisions that the resulting incidence angles are roughly equally
distributed. Then the local effect of the collisions is to force the quantity
$f$ toward an equilibrium. It represents the collision term as the separation
of $f$ locally from its statistical equilibrium $f_{eq}$, taken to have a
Maxwellian form classically (see Appendix A), with respect to a relaxation
time $\tau_{R}$.

$\frac{\partial f}{\partial t}+\boldsymbol{u}\cdot\nabla
f\frac{d}{dt}+\boldsymbol{F}\cdot\nabla^{\boldsymbol{p}}f=\frac{1}{\tau_{R}}\left(f_{eq}-f\right)$
(2.5)

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 momentum and
energy.

The integration of the probability function $f$ and another function $g$ over
momentum space will yield the average value of $g$; $<g>=\int gfd^{3}p$ [49].
Applying this to orders of the fluid velocity yields expressions for the
moments of the fluid; the local mass density $\rho$ (order $0$), the
macroscopic fluid velocity $\bm{u}$ (order $1$), and the total (kinetic plus
internal) energy $\epsilon$ (order $2$).

$\displaystyle\rho=\int f(t,\bm{x},\bm{p})d^{3}\boldsymbol{p},$
$\displaystyle\rho\boldsymbol{u}=\int\boldsymbol{u}f(t,\bm{x},\bm{p})d^{3}\boldsymbol{p},$
$\displaystyle\epsilon=\int\frac{1}{2}|u|^{2}f(t,\bm{x},\bm{p})d^{3}\boldsymbol{p}$
(2.6)

These relations hold when the flow is at the local equilibrium state as well,
where the probability distribution function is the _equilibrium_ probability
distribution function $f^{eq}$. Hence, if the macroscopic moments at
equilibrium are known, $f^{eq}$ can be determined in terms of them. The form
of $f^{eq}$ can be determined along the same lines of reasoning followed by
Maxwell (see Appendix A), but it is straightforward to find an expression of
the equilibrium probability distribution function in terms of velocity from
2.1.

$f^{eq}(t,\bm{x},\bm{v})=\frac{\rho}{m}\left(\frac{1}{2c_{s}^{2}(T)\pi}\right)^{\frac{3}{2}}e^{\frac{-\left|\bm{v}-\bm{u}\right|^{2}}{2c_{s}^{2}(T)}}$
(2.7)

This form introduces the square of the local sound speed as a function of
temperature $T$, $c_{s}^{2}(T)=\frac{T}{m}$, and is used within the BGK
collision operator to complete the Boltzmann transport equation.

When close to local equilibrium, the Boltzmann equation can be expanded in
powers of the Knudsen number, $Kn=\frac{\tau_{R}c}{L}$, where $L$ is the
characteristic length scale of the system and $\tau_{R}$ is the relaxation
time, which, combined with $c$, represents the mean free path of the particles
in the fluid system. The first term in the expansion yields the Euler equation
as follows: Using the Knudsen number as the choice of smallness parameter
$\epsilon$, the probability distribution function can be expanded as

$f=f^{(0)}+\epsilon f^{(1)}+\epsilon^{2}f^{(2)}+...\text{.}$

The zero order term in the expansion is the equilibrium distribution function
$f^{eq}$. Substituting the zero and first orders into 2.5 (in one dimension
for simplicity) we have:

$\displaystyle\frac{\partial}{\partial t}\left(f^{(0)}+\epsilon
f^{(1)}\right)+u\frac{\partial}{\partial x}\left(f^{(0)}+\epsilon
f^{(1)}\right)+\frac{F_{x}}{m}\frac{\partial}{\partial
u}\left(f^{(0)}+\epsilon f^{(1)}\right)$
$\displaystyle=\frac{1}{\tau_{R}}\left(\left(f^{(0)}+f^{(1)}\right)-\epsilon
f^{eq}\right).$ (2.8)

Taking only the zero order of $\epsilon$ will give the equilibrium relation

$\frac{\partial f^{eq}}{\partial t}+u\frac{\partial f^{eq}}{\partial
x}+\frac{F_{x}}{m}\frac{\partial f^{eq}}{\partial
u}=\frac{1}{\tau_{R}}\left(f^{eq}-\epsilon f^{eq}\right)=0.$

Multiply the equation by $u$, integrate over $p$ and apply 2.1 to find

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=-\frac{F_{x}}{m}$

which is the one-dimensional Euler equation. The second term in the expansion
will similarly yield the Navier-Stokes equation.

###  2.2 Turbulence

The phenomenon of turbulence has been widely studied, and the literature is
numerous. A few works relevant to this report include (but are not limited
to), [50], [51] and [52]. A turbulent flow is characterized by chaotic changes
in pressure and velocity throughout the fluid as a result of large kinetic
energy in parts of the fluid that overwhelm the viscous damping. It creates
vortices within the flow at different scales; expressed by the curl of the
fluid velocity; $\bm{\omega}=\nabla\times\bm{u}$. The vortex structures cohere
to each other, clockwise rotating structures to other clockwise structures,
and counterclockwise to counterclockwise, eventually detaching at different
scales at periodic intervals, referred to as vortex shedding. When the
periodic intervals cause irregularities in the current, it can potentially be
turbulent.

It is difficult to predict quantitatively if turbulence will be present in a
fluid system. The Reynolds number is a unitless ratio of the characteristic
length $L$ and average (macroscopic) fluid speed of the system $U$ to the
viscous damping term $\nu$; $Re=\frac{LU}{\nu}$. It is commonly used to
identify the boundary between a turbulent and non-turbulent flow. A large
ratio indicates an over-matched damping term, and turbulence is thought to be
present when it is larger than something on the order of $5000$. If the ratio
is smaller than that threshold, but larger than $10$, the flow is thought to
be in the preturbulent regime. Any flow with a Reynolds number around or below
the order of $10$ is not believed to be turbulent. It is the most commonly
accepted metric to predict the presence of turbulence.

A consequence of vortex coherence and shedding is a creation and migration of
modes in the energy spectrum, also referred to as the power spectrum. A vortex
is formed as some scale in k-space, for example a large scale (small $k$), and
new vortices are created at smaller and smaller scales (larger and larger $k$)
within the rotating flow until it dissipates at the molecular scale. This
cascading vortex creation is detectable in the energy spectrum $E(k)$ through
new modes that appear, having drawn energy from the larger scale modes [53].
The energy spectrum is the energy density contained at a given wave number
magnitude (denoted here simply as $k$), or the contribution to the energy
density by a wave number magnitude [54].

$E\left(k\right)\equiv\frac{d\mathcal{E}\left(k\right)}{dk}$

$\mathcal{E}\left(k\right)$ is the energy density in spectral space defined as
a function of $k$:
$\mathcal{E}\left(k\right)\equiv\frac{1}{2}\langle|\tilde{v}(\bm{k})|^{2}\rangle$,
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.

$\mathcal{E}\left(k\right)=\int_{0}^{\infty}E\left(k^{\prime}\right)dk^{\prime}$

The migration of energy to different modes within the energy spectrum, called
the cascade of the power spectrum [55], or the direct cascade, follows a
$-5/3$rd power law. Kolmogorov [55] reasoned that the energy would cascade to
smaller scales at a rate proportional to the $-5/3$rd power of the wave
number. Explicitly, he calculated through dimensional reasoning that the power
spectrum would decay as $\epsilon^{\frac{2}{3}}\ k^{-\frac{5}{3}}$ where
$\epsilon$ is the energy dissipation rate. The form of this power law was
determined by Kraichnan [56], and has since been observed in simulated and
experimental data [57]. Useful literature describing the cascade of the energy
spectrum includes the following: [58], [59], [60], [61], [62], [63], [64] and
[65].

Specific to a two-dimensional system, a cascade of the energy spectrum from
small scales to larger scales, called the inverse cascade of the energy (or
power) spectrum, can also be detected and is also proportional to the $-5/3$rd
power of the wave number [55][56]. Its presence in a two-dimensional system
can be contributed to conservation of enstrophy $Z$, which is the rotational
analog of energy and is expressed as $Z=\int_{S}|\nabla\times\bm{u}|^{2}$. It
is conserved in two-dimensional systems only. In a 2D homogeneous, isotropic
system, the vorticity equation is [66]:

$\frac{\partial\omega}{\partial t}+\bm{u}\nabla\omega=\bm{g}+\nu\Delta\omega$
(2.9)

where $\omega$ is the normal component of the rotation vector
$\omega=\left(\nabla\times\bm{u}\right)_{\hat{\bm{z}}}$, and $\nu$ is the
kinematic viscosity.

In three dimensions the viscous strain can cause the rotation vector
$\bm{\omega}$ to extend, which can cause vorticity to increase. In a 2D
system, with no forcing term $\bm{g}$ and no viscosity, the vorticity is
conserved completely. Where viscous effects are in play, it can only decrease,
and in proportion to the viscosity. Because of this the cascade of the energy
spectrum can only propagate from an initial driving wavenumber to larger
scales, resulting in an indirect energy cascade. However, and enstrophy
cascade from large scales to small scales (the direct enstrophy cascade) is
possible as the dissipation is directly related to viscosity. Kolmogorov
predicted the log scale slope of this cascade would be $-3$. That is, the rate
of enstrophy dissipation would be proportional to $|k|^{-3}$. A decent amount
of literature on turbulence specific to two-dimensional systems is available,
for example the descriptions in the following: [67], [68], [69], [70], [71]
and [72].

As with the power spectrum, turbulence can be identified through its effect on
fluctuations in the spectrum of the current density $\bm{j}$. When the
vortices begin a pattern of coherence and shedding, the flow creates
fluctuations in the current density with a defining spectral signal in both
frequency space and in wave number space. When the flow reaches a state of
turbulence, the modes spread out from initial prominent modes and the spectrum
becomes broadband. If the broadband signal is present in both frequency space
and in k-space the flow is considered turbulent as well as chaotic [73].
However, if the signal is narrow banded in either space, i.e. the modes do not
migrate temporally or spatially, then the flow, though still potentially
chaotic, is not considered turbulent.

Figure 2.1: A sketch depicting how the spectra of the current density can
indicate the presence of a turbulent flow. If the spectrum in both frequency
space and wavenumber space are broadband, the flow can be turbulent.

###  2.3 Charged Quasi-Particle Flow

Of focused interest in this submission is the transport of charge carrying
quasi-particles in solid matter. In particular, we concentrate on the class of
two-dimensional solids (solids that are a single atom thick) with an atomic
lattice structure that produces a linear dispersion relation in its energy
band. Graphene offers one of the best examples of this class of solid because
of its simple atomic structure (a single layer of carbon atoms arranged in a
honeycomb structure) and its accessibility to experiment. It can be wrapped up
into zero-dimensional fullerenes, rolled into single-dimensional carbon
nanotubes, and be stacked together forming the familiar three-dimensional
graphite that can be found in a pencil. The properties of its band structure
were explored by Wallace [74] in 1946, but it wasn’t until relatively recently
that it has been closely examined in experimental settings when it was
discovered that it is reasonably simple to produce a useful sample [23].

The $sp^{2}$ hybridized orbit of the electrons in the atoms create resilient
bonds between in-plane carbon atoms. Two electrons in a $p$ orbit and one in
an $s$ orbit form $\sigma$ bonds with their neighbors, filling up the shell.
The last $p$ orbit is perpendicular to the plane leaving the shell it occupies
half-filled, and the atom charge neutral. The atoms arrange themselves into a
two-dimensional hexagonal pattern with a lattice constant of $a=1.42{\AA}$
[75].

Figure 2.2: Honeycomb lattice of graphene where the carbon atoms are
classified as belonging to either a red or blue sublattice where each carbon
atom’s nearest neighbor is classified in its companion lattice.

Figure 2.3: Brillouin Zone in reciprocal space about a carbon atom within a
honeycomb lattice of graphene. The vector $\boldsymbol{K}$ defines the
boundary of the Brillouin zone and, along with the rest of the atoms in the
lattice, creates a reciprocal lattice.

Wallace determined the band structure of this configuration to be linear,
approximated by

$\epsilon(\boldsymbol{k})=\pm\hbar v_{F}|\boldsymbol{k}|,$ (2.10)

where $v_{F}\approx 10^{6}m/s$ is the Fermi velocity. This dispersion relation
is derived from a ‘second quantization’ approach from the Hamiltonian, given
in [76] as

$H=-t_{1}\sum_{i,j}c_{i\alpha}^{\dagger}c_{j\alpha}-t_{2}\sum_{i}\left(n_{i\uparrow}+n_{i\downarrow}\right).$
(2.11)

The constants $t_{1}$ and $t_{2}$ are short range ‘hopping’ amplitudes from a
site (indexed by $i$) to the nearest neighbor and next-nearest neighbor,
respectively. The operators $c_{i\alpha}$ and $c_{i\alpha}^{\dagger}$ are the
creation and annihilation operators, and the index $\alpha$ indicates the
electrons spin (up $\uparrow$, or down $\downarrow$). The number operators
$n_{i\uparrow}\equiv c_{i\uparrow}^{\dagger}c_{i\uparrow}$ and
$n_{i\downarrow}\equiv c_{i\downarrow}^{\dagger}c_{i\uparrow}$ are interpreted
physically as the electron density. For reference, the commutation relations
are as follows,

$\displaystyle\left\\{c_{i\alpha},c_{j\beta}^{\dagger}\right\\}$
$\displaystyle=c_{i\alpha}c_{j\beta}^{\dagger}+c_{j\beta}^{\dagger}c_{i\alpha}=\delta_{ij}\delta_{\alpha\beta}$
$\displaystyle\left\\{c_{i\alpha},c_{j\beta}\right\\}$
$\displaystyle=c_{i\alpha}c_{j\beta}+c_{j\beta}c_{i\alpha}=0$

With this model the lattice is broken up into two different triangular sub-
lattices so that all nearest neighbors of an individual site are part of the
companion sub-lattice, as shown in figure 2.3 with the two different sub-
lattices shown in red and blue. The basis $e_{n}$ is chosen in the direction
of the nearest neighbors of atoms in the green sub-lattice.

$\begin{gathered}e_{1}=\left(1,0\right),\ \ \
e_{2}=\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right),\ \ \
e_{2}=\left(-\frac{1}{2},-\frac{\sqrt{3}}{2}\right)\\\ e_{i}\cdot
e_{j}=-\frac{1}{2}\text{ for }i\neq j,\ \ \ e_{1}+e_{2}+e_{3}=0\end{gathered}$

A Fourier transform on 2.11 migrates the Hamiltonian to momentum space. The
creation operator transforms as

$\displaystyle
c_{A\alpha}\left(\boldsymbol{k}\right)=\frac{1}{\sqrt{\mathscr{T}}}\sum_{i\in
A}c_{i\alpha}e^{-i\boldsymbol{k}\cdot\boldsymbol{e_{i}}}.$

Applying the transformation to the complete form yields

$\displaystyle
H=-t_{1}\sum_{\boldsymbol{k}}\left(e^{i\boldsymbol{k}\cdot\boldsymbol{e_{1}}}+e^{i\boldsymbol{k}\cdot\boldsymbol{e_{2}}}+e^{i\boldsymbol{k}\cdot\boldsymbol{e_{3}}}\right)c_{A\alpha}^{\dagger}(\boldsymbol{k})c_{B\alpha}^{\dagger}(\boldsymbol{k})+\text{H.c.}$
$\displaystyle-
t_{2}\sum_{\boldsymbol{k}}\left(c_{A\alpha}^{\dagger}(\boldsymbol{k})c_{A\alpha}^{\dagger}(\boldsymbol{k})+c_{B\alpha}^{\dagger}(\boldsymbol{k})c_{B\alpha}^{\dagger}(\boldsymbol{k})\right).$

The ”H.c.” term refers to the Hermitian conjugate of the previous term. This
tensor equation is restated using Pauli spin matrices $\tau^{x}$ and
$\tau^{y}$

$\displaystyle
H=\sum_{k}c^{\dagger}(\boldsymbol{k})[[-t_{2}-t_{1}\left(cos(\boldsymbol{k}\cdot
e_{1})+cos(\boldsymbol{k}\cdot e_{2})+cos(\boldsymbol{k}\cdot
e_{3})\right)\tau^{x}$ (2.12) $\displaystyle+\left(sin(\boldsymbol{k}\cdot
e_{1})+sin(\boldsymbol{k}\cdot e_{2})+sin(\boldsymbol{k}\cdot
e_{3})\right)\tau^{y}]]c(\boldsymbol{k}).$

The energy eigenvalues are found to be

$-t_{2}\pm|e^{i\boldsymbol{k}\cdot e_{1}}+e^{i\boldsymbol{k}\cdot
e_{2}}+e^{i\boldsymbol{k}\cdot e_{3}}|.$ (2.13)

Diagonalizing 2.12 yields the dispersion relation

$\begin{gathered}E_{\pm}(\boldsymbol{k})=\pm
t_{1}\sqrt{3+f(\boldsymbol{k})}-t_{2}f(\boldsymbol{k}),\\\
f(\boldsymbol{k})=2cos(\sqrt{3}ak_{y})+4cos(\frac{\sqrt{3}}{2}ak_{x})cos(\frac{\sqrt{3}}{2}ak_{y}),\end{gathered}$
(2.14)

where recall $a$ is the distance between the carbon atoms.

The Brillouin zone is defined in momentum space around an atom where the
k-vectors begin to repeat themselves, delineating the lattice cells in
reciprocal space. In graphene it forms a hexagon about the atom, separated
from it by a vector $\boldsymbol{K}$ (fig. 2.3) [77]. If we presume the
negligible effects of the next-nearest neighbor on the energy (i.e.
$t_{2}=0$), and we limit the formulation to a single spin (i.e. $\uparrow$),
we can expand 2.14 about a vector $\boldsymbol{q}$ that is small compared to
$\boldsymbol{K}$ such that $\boldsymbol{k}=\boldsymbol{K}+\boldsymbol{q}$ and
$\boldsymbol{q}\ll\boldsymbol{K}$ and the dispersion relation is simplified to

$\displaystyle E_{\pm}(\boldsymbol{k})$ $\displaystyle\approx\pm
t_{1}\sqrt{\left(\frac{3a}{2}q_{x}\right)^{2}+\left(\frac{3a}{2}q_{y}\right)^{2}+O(q_{x}^{3})+O(q_{y}^{3})}$
(2.15) $\displaystyle=\pm v_{F}|\boldsymbol{q}|+O[q^{2}],$

where $v_{F}=\frac{3at_{1}}{2}\approx 10^{6}m/s$. This dispersion indicates
there are two points in the Brillouin zone where the energy is zero. These are
called Dirac points, referred to here as $\boldsymbol{Q}$ and
$-\boldsymbol{Q}$. The surface moving away from these points is linear close
to the points and forms a conical shape called a Dirac cone [78] [79].

Figure 2.4: Fermi surface in the electric band structure of graphene,
featuring Dirac cones near the low energy points called Dirac points.

If we define a continuous Fermi field in the area immediately surrounding the
points, we can define an 8 component continuum Fermi field as follows.

$\displaystyle C_{A1\alpha}(\boldsymbol{k})$
$\displaystyle=\sqrt{G}c_{A\alpha}\left(\boldsymbol{Q}_{1}+\boldsymbol{k}\right)$
$\displaystyle C_{A2\alpha}(\boldsymbol{k})$
$\displaystyle=\sqrt{G}c_{A\alpha}\left(-\boldsymbol{Q}_{1}+\boldsymbol{k}\right)$
$\displaystyle C_{B1\alpha}(\boldsymbol{k})$
$\displaystyle=\sqrt{G}c_{B\alpha}\left(\boldsymbol{Q}_{1}+\boldsymbol{k}\right)$
$\displaystyle C_{B2\alpha}(\boldsymbol{k})$
$\displaystyle=\sqrt{G}c_{B\alpha}\left(-\boldsymbol{Q}_{1}+\boldsymbol{k}\right)$

In this, $G$ is the total area of the honeycomb lattice, and $\boldsymbol{k}$
is assumed to be small. The $\alpha$ index refers to the spin, $A$ and $B$
identify the sub-lattice, and the indexes $1$ and $2$ refer to the ‘valley’,
or conical area in the Dirac cones. Again employing Pauli spin matrices for
spin $(\sigma^{a})$, sub-lattice $(\tau^{a})$ and valley $(\rho^{a})$ the
continuum Hamiltonian is

$H=\int\frac{d^{2}k}{4\pi^{2}}C^{\dagger}(\boldsymbol{k})\left(v_{F}\tau^{y}k_{x}+v_{F}\tau^{x}\rho^{z}k_{y}\right)C(\boldsymbol{k})$
(2.17)

where, again, $v_{F}$ is the Fermi velocity. Diagonalizing 2.17 leads to a
relativistic spectrum

$\pm\sqrt{k_{x}^{2}+k_{y}^{2}}$

corresponding to the values of 2.13 near the Fermi points.

In the region close to the Dirac points in momentum space where the number
density of the electrons approaches zero, one can find long-lived massless
quasi particles that carry an electric charge. The mean free path of the
particles is within the fluid regime, and collectively they may be considered
a massless, relativistic, viscous charged fluid [16] [17]. Note that the
velocity of the quasi-particles is not in the relativistic regime, but they
are relativistic in the sense that they are massless and their momentum
distribution in the fluid is of the Jüttner form. Evidence of this viscous
charged fluid is noted by [80] who identified micrometer sized vortices that
cause an inverted electric field which they attribute to a viscous effect.
This field inversion is most prominent near the Dirac points.

##  Chapter 3 The Lattice Boltzmann Method

###  3.1 The Standard/Classical LBM

The standard/classical 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. It replaces the particles in the lattice with the
probability distribution function described by the Boltzmann equation (2.5),
and takes advantage of the near-equilibrium assumption of the BGK collision
term, sacrificing collision details for the fluid’s macroscopic moments. The
dynamics of the moments are then able to be modeled accurately and efficiently
with an appropriately defined computational model. The Lattice Boltzmann
Method or Lattice Boltzmann Model refers to a family of computational methods
that model the dynamics of the probability distribution function through the
Boltzmann equation projected onto a discrete spatial lattice with a
superimposed discrete velocity lattice. This provides the model with the
minimum amount of information required to solve the Boltzmann equation
accurately and efficiently. There is an abundance of literature on the
implementation and application of LBM models; see for instance [34], [81]
[82], [83], [84] and [85].

In an LBM, the Boltzmann equation is discretized through Chapman-Enskog
analysis [86] [87]. Space is discretized into a spatial lattice where all the
information in the probability distribution function for a discrete volume is
contained within each node, separated from the other nodes by a lattice
spacing in one, two, or three dimensions. The momentum space is also discrete
at every node using a quadrature analysis technique that identifies only the
discrete speeds necessary to preserve the moments. In that way, the LBM
lattice can be considered a momentum lattice.

For each discrete time step the moments of the fluid are evaluated locally at
each spatial node. The moments are obtained through integrals over momentum
(or velocity) space such as those defined in 2.1. The integrals are evaluated
with a quadrature integral evaluation technique. This technique requires
dividing momentum space into discrete momenta, which is how the momentum
lattice is defined. Therefore it is the structure that the LBM is built upon.
Considering that, a closer look at quadrature is appropriate.

For simplicity, consider a one-dimensional function of $x$, $h^{N}(x)$, that
is a non-trivial polynomial of degree $N$, and separable as
$h^{N}(x)=\omega(x)g(x)$. For this expression $\omega(x)$ is some appropriate
weighting function. An integral of the form
$\int_{-\infty}^{\infty}\omega(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...q$ [88] [89] [90].

$\int_{-\infty}^{\infty}\omega(x)g(x)dx\approx\sum_{i=1}^{q}w_{i}g(x_{i})$

The accuracy of the integration depends on the number of abscissae and how
they are chosen. If we choose the abscissae to be the $n$ roots of the $n$th
polynomial $P^{(n)}$ from a class of orthogonal polynomials, and
$n\geq\frac{1}{2}(N+1)$, the resulting summation will evaluate the integral
_exactly_.

$\int_{-\infty}^{\infty}\omega(x)g(x)dx=\sum_{i=1}^{n}w_{i}g(x_{i})$ (3.1)

The weights $w_{i}$ are determined based on the choice of polynomials, which
are, in general, chosen based on the limits of integration. For an integration
over $[-1,1]$, Legendre polynomials are useful. For the range
$\left[0,\infty\right)$, Laguerre polynomials are best suited. And, for the
integral limits described in (3.1), $(-\infty,\infty)$, Hermite polynomials
are the appropriate choice.

When Hermite polynomials are used the discrete weights are determined by

$w_{i}=\frac{n!}{\left(nH^{\left(n-1\right)}\left(x_{i}\right)\right)^{2}},$

where $H^{(n)}$ is the $n$th order Hermite polynomial. The Gauss-Hermite
quadrature integral relation is generalized to $d$ dimensions [81] as

$\int_{-\infty}^{\infty}{d^{d}x\omega\left(\bm{x}\right)g\left(\bm{x}\right)}=\sum_{i=1}^{n^{d}}{w_{i}g\left(\bm{x}_{i}\right)}$
(3.2)

This quadrature technique implies that momentum space can be broken into
discrete momenta at each discrete point in space to evaluate the integrals of
the fluid moments. As a simplified example, consider the evaluation of the
moments in 2.1 with a one-dimensional, classical probability distribution
function.

$\displaystyle\rho=\int fdp=\int w(p)g(p)dp=\sum_{i}w_{i}g_{i}=$
$\displaystyle\sum_{i}f_{i},$ $\displaystyle\rho u=\int fvdp=\int
w(p)g(p)vdp=\sum_{i}w_{i}g_{i}v_{i}=$ $\displaystyle\sum_{i}f_{i}v_{i},$
$\displaystyle\sigma=\int fv^{2}dp=\int
w(p)g(p)v^{2}dp=\sum_{i}w_{i}g_{i}(v_{i})^{2}=$
$\displaystyle\sum_{i}f_{i}v_{i}^{2}$

In these summations $f_{i}$ is probability distribution function for a
discrete momentum $i$, determined using quadrature and incorporating the
weight.

The finite summations faithfully reproduce the moment integrals (2.1) for
sufficiently large $n$, implying the discrete abscissae can be used to
discretize momentum space. That is, the roots of the three-dimensional Hermite
polynomial of at least the highest order of the moments, the abscissae
$x_{i}$, determine the discrete momenta that are required to connect the nodes
of the lattice. The chosen number of abscissae $n$ need only be one less and
half the order of the highest fluid moment; $n=\frac{N-1}{2}$. Different LBM
formulations use different numbers of speeds depending on precision and
computational efficiency requirements. Commonly used values for $n$ with
three-dimensional lattices are $19$ and $27$ in configurations named D3Q19 and
D3Q27, respectively.

Using this discretization formula, the continuous probability distribution
function is then re-expressed in discrete form: $f(t,\bm{x},\bm{p})\rightarrow
f_{i}(t,x_{j})$, 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 Lattice Boltzmann Equation is then:

$f_{i}\left(t+\Delta t,x_{j}+v_{i}\Delta
t\right)-f_{i}\left(t,x_{j}\right)=\Delta
tC\left[f_{i}\left(t,x_{j}\right)\right].$ (3.3)

The collision term on the right-hand side is

$C\left[f_{i}\left(t,x_{j}\right)\right]=\frac{1}{\tau_{R}}\left(f^{eq}_{i}-f_{i}\right),$
(3.4)

where $f^{eq}_{i}$ is the local equilibrium distribution function for lattice
node $i$. Notice that, with this formulation, the full value of the
probability distribution function at each node is never needed. The value at
each discrete momentum is determined independently by the same $f_{i}$ of
neighboring nodes and the equilibrium distribution function $f^{eq}_{i}$,
which is determined by the integration over the $f_{i}$ through quadrature.
However, the quadrature technique is valid only for polynomials, so the
equilibrium probability distribution function must be expanded in terms of the
orthogonal polynomials chosen.

Figure 3.1 depicts a single node of a three dimensional lattice with 19
speeds. Each node contains 19 different probability distribution function
values $f_{i}$, one for each speed at that lattice location.

Figure 3.1: Single node of a D3Q19 configured lattice. The node at the center
of the cube is connected to its neighboring nodes by 19 discrete speeds
(including 0).

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_{\perp}$ and frequency $\omega_{\perp}$ [49].

$\bar{\bm{x}}=\frac{\bm{x}}{R_{\perp}},\ \ \bar{t}=t\omega_{\perp},\ \
\bar{\bm{v}}=\frac{\bm{v}}{R_{\perp}\omega_{\perp}},\ \
\bar{\bm{u}}=\frac{\bm{u}}{R_{\perp}\omega_{\perp}}$ (3.5)

The rescaled quantities are then easily converted to the scale of the domain
being modeled when they are recorded. With both Euclidean space and momentum
space discretized, the Lattice Boltzmann Method then simulates the time
evolution of a hydrodynamic system with discrete time steps following this
process:

  1. 1.

Initialize the local equilibrium distribution function $f_{i}^{eq}$ 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 function, $f_{i}^{*}=f_{i}^{eq}$.

  2. 2.

“Stream” the new distribution function $f_{i}^{*}$ to the neighboring node
connected by the discrete speed, accounting for boundary conditions.

  3. 3.

Calculate the macroscopic moments based on the local state at each node,
optionally recording them for analysis.

  4. 4.

Calculate the equilibrium probability distribution function $f_{i}^{eq}$ based
on the moments.

  5. 5.

Perform the collision calculation locally using the equilibrium probability
distribution function $f_{i}^{eq}$ and the local probability distribution
function $f_{i}$ to find the new probability distribution function:
$f_{i}^{*}=f_{i}+\Delta tC$.

  6. 6.

Repeat from step 2 for each time step in the simulation.

This technique is advantageous computationally since, though the dynamics of a
hydrodynamic system are non-linear (see (2.5)), 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.

###  3.2 Relativistic Lattice 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. Some variants have been proposed to adjust the Lattice
Boltzmann Method framework to model relativistic hydrodynamic flows. These are
less abundant in literature than their classical counterpart, but detailed
descriptions can be found in the following works: [91], [92], [93], [94],
[95], [96], [97], [98], [99] and [100]. One such variation is proposed by
Romatschke, Mendoza and Succi in 2011 [42]. This Relativistic Lattice
Boltzmann Method (RLBM) adaptation is based on the dynamics of the
relativistic version of the probability distribution function,
$f=f(x^{\mu},p^{\nu})$, where $x^{\mu}$ is the position 4-vector and $p^{\nu}$
is the momentum 4-vector. The dynamics of the probability distribution
function are described by the relativistic Boltzmann equation with a
relativistic analog of the BGK collision term [101].

$\left[p^{\mu}\nabla_{\mu}-\Gamma_{\mu\nu}^{\lambda}p^{\mu}p^{\nu}\partial_{\lambda}^{(p)}\right]f=-\frac{p^{\mu}u_{\mu}}{\tau_{R}}\left(f-f_{J}^{eq}\right)$
(3.6)

For this equation, and for the remainder of this treatment, the units are
natural; $c=k_{B}=\hbar=1$. $\nabla_{\mu}$ is the covariant derivative, and
$\Gamma_{\mu\nu}^{\lambda}$ is the Christoffel symbols which are given by
derivatives of the underlying metric tensor $g_{\mu\nu}$. For this report we
assume Minkowski space-time configuration, so that it is space-time is flat
and the second term on the left-hand side is neglected. The macroscopic fluid
velocity is $u^{\nu}=\gamma\left(1,v^{i}\right)$, where
$\gamma=\left(1-v^{2}\right)^{-\frac{1}{2}}$ is the Lorentz factor and $v^{i}$
is the fluid 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=\sqrt{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_{J}^{eq}$ which is the
Jüttner distribution.

$f_{J}^{eq}=Z^{-1}e^{\frac{-p_{\mu}u^{\mu}}{T}}$ (3.7)

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
energy-momentum (energy-stress) tensor $T^{\mu\nu}$ which can be obtained by
integration of $f$ over the 4 components of momentum space [42].

$\displaystyle T^{\alpha\beta}\equiv$ $\displaystyle\int d\chi
p^{\alpha}p^{\beta}f(t,x,p)$ $\displaystyle\equiv$
$\displaystyle\int\frac{d^{4}p}{\left(2\pi\right)^{3}}\delta\left(p^{\mu}p_{\mu}-m^{2}\right)2H\left(p^{0}\right)p^{\alpha}p^{\beta}f(t,x,p)$
(3.8)

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, $f(x^{\mu},p^{\nu})\rightarrow f(t,x,p)$.)

Using the Jüttner form of the equilibrium probability distribution function,
integration of (3.6) at equilibrium produces the equilibrium energy-momentum
tensor $T_{eq}^{\mu\nu}$ with the familiar relation (see (2.3))

$T_{eq}^{\mu\nu}=\left(\epsilon+P\right)u^{\mu}u^{\nu}+Pg^{\mu\nu}.$ (3.9)

Here both energy density $\epsilon$ and pressure $P$ are functions of
temperature. The viscosity coefficient, represented in the collision term, is
$\eta=\tau_{R}\frac{\epsilon+P}{5}$, with relaxation time
$\tau_{R}=5\frac{\eta}{s}T^{-1}$, where $s$ is the entropy density [102]. We
apply the canonical, so-called Landau-Lifshitz condition, where
$u_{\mu}T^{\mu\nu}\equiv\epsilon u^{\nu}$, to provide a relationship between
the energy density and the fluid velocity in the rest frame of reference.

###  3.3 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 are discretized in the same
manner, but the momentum 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 generator function of the Hermite
polynomials; $e^{-\frac{v^{2}}{T}}\ \rightarrow
e^{-\frac{p^{\mu}u_{\mu}}{T}}$. Therefore, the discretization of the velocity
space is less straight forward. To address this, the equilibrium probability
distribution function is re-expressed as an expansion about powers of
$\frac{|\bm{p}|u^{0}}{T_{0}\theta}$.

$e^{-\frac{p\cdot
u}{T}}=e^{-\frac{|p|u^{0}}{T_{0}\theta}}\sum_{n}\left(\frac{\bm{p}}{|\bm{p}|}\right)^{n}\left(\frac{|\bm{p}|u^{0}}{T_{0}\theta}\right)^{n}\frac{\left(\bm{u}\right)^{n}}{\left(u^{0}\right)^{n}n!}$
(3.10)

Here the momentum is expressed as unit vectors,
$\bm{v}=\frac{\bm{p}}{|\bm{p}|}$, and $\theta$ is the scaled temperature with
respect to a reference temperature $T_{0}$, $\theta=\frac{T}{T_{0}}$. In this
form we can express the relativistic probability distribution function in
terms of a (different) set of orthogonal polynomials as:

$f(t,x,p)=e^{\frac{-p^{0}}{T_{0}}}\sum_{n}P_{i_{1}...i_{n}}^{(n)}\left(\bm{v}\right)a_{i_{1}...i_{n}}^{(n)}\left(t,x,\frac{p^{0}}{T_{0}}\right)$
(3.11)

where $p^{0}=|\bm{p}|$. The chosen polynomials
$P_{i_{1}...i_{n}}^{(n)}\left(\bm{v}\right)$ are three dimensional polynomials
about the unit vector velocity determined by orthogonality with respect to the
angular integral $\int\frac{d\Omega}{4\pi}$. Their properties are listed in
Appendix B. The coefficients in (3.11) are determined using orthogonality
conditions:

$\displaystyle a^{(0)}\left(t,x\right)=\frac{1}{\left(\alpha\right)!}\int
d\bar{p}\
\bar{p}^{\alpha}\int\frac{d\Omega}{4\pi}fP^{(0)}\left(\bar{p}\right),$
$\displaystyle a_{i}^{(1)}\left(t,x\right)=\frac{3}{\left(\alpha\right)!}\int
d\bar{p}\
\bar{p}^{\alpha}\int\frac{d\Omega}{4\pi}fP_{i}^{(1)}\left(\bar{p}\right),$
$\displaystyle
a_{ij}^{(2)}\left(t,x\right)=\frac{15}{2\left(\alpha\right)!}\int d\bar{p}\
\bar{p}^{\alpha}\int\frac{d\Omega}{4\pi}fP_{ij}^{(2)}\left(\bar{p}\right)$
(3.12)

where $\bar{p}=p^{0}/T_{0}$ is the scaled momentum magnitude. The parameter
$\alpha$ 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 [42].

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
(3.1). 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 (the connection to neighboring nodes is off-lattice).
This is overcome computationally using a linear interpolation scheme. An
example of a spherical lattice node is shown in figure 3.2.

Figure 3.2: A sketch of a spherical lattice node in an RLBM lattice. The
momentum lattice connecting the nodes is spherical instead of box-like.

It remains to obtain the macroscopic moments of the fluid system from the
probability distribution function by evaluating the energy-momentum tensor.
Using equation (3.3) and (3.2) we find:

$T^{\mu\nu}=\frac{3T_{0}^{4}}{\pi^{2}}\left(\begin{matrix}a^{(0)}&a_{i}^{(1)}\\\
a_{j}^{(1)}&a_{ij}^{(2)}+\frac{1}{3}\delta^{ij}a^{(0)}\end{matrix}.\right)$
(3.13)

The fluid velocity $u^{\mu}$ is determined by the eigenvalue problem set up by
the Landau-Lifshitz condition, $u^{\mu}T_{\mu}^{\nu}=\epsilon u^{\nu}$, with
the energy density as the eigenvalue $\epsilon$. The temperature is determined
from the relation $\epsilon=\frac{3T^{4}}{\pi^{2}}$, and pressure $P$ (denoted
with a capital $P$ in this section) is found from the equation of state, which
is that of an ideal gas in a relativistic system, $P=\frac{\epsilon}{3}$ (see
Appendix E).

##  Chapter 4 Part I: Reproduction of a 2D turbulent flow with the RLBM

We apply the RLBM hydrodynamic modeler to reproduce the flow of an ideal
isotropic fluid system induced to turbulence with a random stirring force,
generated in momentum space. A forcing term is introduced into the
relativistic Lattice Boltzmann equation through the discretization of the
driving term. The new driving term is populated with a stirring mechanism
driving the fluid in random directions and random (bounded) magnitudes in
momentum space at defined frequency ranges to induce an energy cascade,
analogous to forced-turbulence studies using the classical LBM for non-
relativistic fluids. The modeler’s ability to induce a turbulent response in
the flow is evaluated.

###  4.1 The Driving Term

A modification to the RLBM described by [42] is required in order to apply the
stirring force to the model. A turbulent flow is created and sustained by the
introduction of a suitably large externally imposed stirring force, zero
averaged in the spatial component; $F^{\mu}=(0,F^{i})$. The force-included
relativistic lattice Boltzmann equation in flat space-time is

$\left[p^{\mu}\nabla_{\mu}-F^{\mu}\nabla_{\mu}\right]f=-\frac{p^{\mu}u_{\mu}}{\tau_{R}}\left(f-f_{J}^{eq}\right).$
(4.1)

As with the probability distribution function and the equilibrium probability
distribution function, the driving term must be expressed in terms of the
orthogonal polynomials

$F^{\mu}\partial_{\mu}f\approx
e^{-\bar{p}}\sum_{n}{a^{(n)}P^{\left(n\right)}}$ (4.2)

The orthogonality relations again determine the expansion coefficients given
in (3.3).

The projected form of the driving term, evaluated in terms of the probability
distribution function and the fluid velocity is found to be (see Appendix C).

$\displaystyle F_{i}\partial_{p}^{i}f\approx\frac{e^{-\bar{p}}}{T_{0}}$
$\displaystyle\left[\frac{3}{12}F_{i}\int d\chi\
v^{i}f-\frac{1}{2}F^{i}v^{i}\int d\chi\ f\right.$
$\displaystyle+\left.\frac{5}{4}F_{i}v_{j}v_{k}\int d\chi\
fv^{i}v^{j}v^{k}-\frac{5}{2}F_{i}v_{i}v_{k}\int d\chi\ fv^{k}\right],$ (4.3)

where the integration, $\int
d\chi=\int_{0}^{\infty}d\bar{p}\int\frac{d\Omega}{4\pi}$, can be evaluated
using the quadrature technique described. The integrals in (4.1) are then
expressed as expansions with the expansion coefficients in (3.3) 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_{i}^{*}=f_{i}+\Delta tC[f]+\Delta
tF^{\mu}\nabla_{\mu}f$. With this modification in place the RLBM is equipped
to produce turbulence and to model a turbulent relativistic hydrodynamic
system.

###  4.2 Model and Methods

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 with a D3Q27
configuration by limiting the number of nodes in one spatial dimension to 1.
The equation of state used in the numerical scheme, that of a three
dimensional ideal gas, is applicable for this two-dimensional model given the
lattice configuration and the two-dimensional forcing scheme. The boundaries
of the lattice are periodic and the initial state sets the probability
distribution of the momenta at each point as flat; i.e. equivalent at each
node on the lattice. Simulations were conducted with lattice size $N$ set to
$128,256,512,$ and $1024$ and the lattice spacing $\alpha$ set to
$0.32,0.16,0.08$ and $0.04$ lattice units respectively, so that lattice
configurations model systems with constant volume. The viscosity, represented
as the ratio of viscosity to entropy $\frac{\eta}{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}\leq k<k_{max}$, leaving the inertial range as $0<k<k_{f}$. Most tests
were performed with $k_{f}$ set to $6$ in lattice units, 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 applied to the relativistic
lattice Boltzmann equation as a part of the forcing term. Each simulation is
run in $0.05$ time unit increments for $50.0,100.0,200.0$ or $400.0$ lattice
time units. Macroscopic variable calculation is conducted at each time step
and recorded every few time steps.

The fluid velocities recorded are filtered to eliminate the contribution of
the driven spectral range. This provides a better analysis of the effects of
the propagation of energy into the inertial spectral range, but note that
recorded fluid velocities represent only a portion of full average fluid
velocities. The velocities obtained for each test varied based on the size of
the stirring force filter and the maximum magnitude applied. They were
observed to be stable and consistent after a sufficient amount of time (Fig
4.1). The achievement of stable velocities allows the filtered average
velocity to be effectively tuned as a parameter so as to maintain a consistent
value between simulations with different configurations. It was tuned by
selecting an appropriate stirring force maximum magnitude $|F|$. This was done
empirically by comparing the average fluid velocity of simulations with
different forcing configurations for a given lattice and matching it with
those of other lattices. Figure 4.2 shows an example of a line of constant
velocity between four different lattice configurations with different forcing
filters.

Figure 4.1: Filtered average fluid velocity throughout the execution of a
simulation at different maximum driving force magnitudes. The simulations were
executed with $256^{2}$ lattice nodes and a lattice spacing of $0.31$. The
average velocities are stable after about 32 simulated time units. Figure
reproduced from [103]. Figure 4.2: A constant filtered average velocity is
determined across multiple lattice configurations empirically using a plot as
shown. A variation of the lattice spacing in the lattice configuration changes
the number of discrete wave numbers, changing the size of the driving filter
and affecting the average velocity of the simulation. A constant velocity
across multiple lattice configurations is determined by finding a common
constant velocity simulated with varying maximum force magnitudes. Figure
reproduced from [103].

###  4.3 Results

For modest driving force magnitudes, the numerical model produced turbulence
inducing fluid velocities of $\sim 0.01$ (all values in this section are in
lattice units unless otherwise specified) 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. Indeed, in the absence of
viscosity, the vortices can grow unbounded so that, in many computational
fluid dynamic codes modeling an Euler fluid a viscosity is induced numerically
to ensure stability. 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. 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 did not affect the
stability of models with smaller lattice spacings. The maximum viscosity was
found to be $0.18$ at the same lattice spacing, also determined by stability.
Some other instabilities were observed as a result of excessively high fluid
velocities present in the driven range, which are required to induce energy
propagation to the inertial range.

With adequate fluid velocity the RLBM numerical scheme shows evidence of the
spectral energy propagation caused by turbulence, reflecting the slope of the
energy in k-space described by Kolmogorov. Figure [4.3] shows the power
spectrum for systems modeled with different lattice configurations maintaining
a constant volume. Lattice configurations with $N$ set to $128,256,512$, and
$1024$ are shown with lattice spacings of $0.32,0.16,0.08$, and $0.04$
respectively. The driving filter is applied at $k_{f}=6.0$, and at a maximum
forcing magnitude determined for each to produce a constant velocity of
$0.015$. The left panel is the energy spectrum at log scale comparing the
slope of the power spectrum against Kolmogorov’s predicted slope of $-5/3$.
The right panel is a “zoomed in” view of the same spectrum divided by the
Kolmogorov slope in order to highlight the conformity. The energies of the
various configurations are scaled independently so that they align on the plot
for comparison.

Figure 4.3: The energy spectrum for simulations with lattice configurations of
$N\alpha=(128)(0.32),(256)(0.16),(512)(0.08)$, and $(1024)(0.04)$, keeping a
constant volume. The energy for each lattice configuration is scaled to a
similar magnitude so the curves can be easily compared. The left panel shows
the slope of the energy spectrum at log scale as compared to the expected
power law of $-5/3$ in the inertial range. Part of the driven range is visible
at $k>6.0$, marked by a discontinuity in the slope to a much higher energy
value. The right panel is a ”close up” of the same divided by the expected
slope. Close conformity to the $-5/3$ slope is noted, with deviations at
larger scales, and for configurations with larger lattice spacing. Figure
reproduced from [103].

A power spectrum cascade was observed with a reasonable reproduction of
Kolmogorov’s power law in the inertial range for mid and large wave numbers.
Deviations at large scales ($\sim N\alpha$) is noted, and is presumed to be a
product of the limited discrete 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 $\alpha=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 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 4.4 compares the same lattice
configurations as in figure 4.3, but driven with a forcing configuration
leaving the inertial range from $k=0$ to $k=14$. This leaves the $256^{2}$
node configuration simulation driving at only a third of the accessible wave
numbers which disrupted the spectral slope. But the spectra of the larger
lattice configurations show a much better conformity to the power law in a
larger portion of the inertial range. Simulations also showed that the larger
the inertial range, the more kinetic energy injected into the system is
required to achieve the energy propagation. Indeed, a fluid velocity of
$\approx 0.25$ was stable for a $N\alpha=(1024)(0.04)$ configuration with a
forcing configuration of $k_{f}=60$ and $|F|=6000$, where the driven range is
only 30% of the available spectrum.

Figure 4.4: The energy spectrum of the inertial range $0<k<14.0$ for
simulations with lattice configurations of $N\alpha=(256)(0.16),(512)(0.08)$,
and $(1024)(0.04)$, scaled for comparison. The left panel shows the slope of
the energy spectrum at log scale as compared to the expected power law of
$-5/3$, and the right panel is a ”close up” of the same divided by the
expected slope. A larger inertial range shows improved conformity. Figure
reproduced from [103].

The energy spectrum of a larger velocity flow is compared to a smaller
velocity flow for a $N\alpha=(1024)(0.16)$ lattice configuration in figure
[4.5]. 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.

Figure 4.5: The energy spectrum of a simulation with lattice configuration
$N\alpha=(1024)(0.16)$ with large and small filtered average velocity. The
left panel compares the energy spectrum at log scale with the expected slope
of $-5/3$. The right panel is the energy spectrum of the larger velocity
simulation divided by the smaller velocity simulation. The magnitude of the
fluid velocity does not have an effect on the inverse energy cascade, provided
the velocity is sufficient to produce turbulence. Figure reproduced from
[103].

Figure [4.6] shows the energy spectrum for systems of increasing size, modeled
by an increasing number of lattice nodes with constant lattice spacing. The
conformity to the expected slope mimics the simulations performed with a
constant volume implying the model’s propagation of energy is not dependent on
the number of nodes, a non-physical parameter. Figure (4.7) (left) shows the
scaled ratio of $E(k)$ observed for a simulation with a lattice spacing of
$0.32$ to $E(k)$ for a simulation with a lattice spacing of $0.16$. The
absence of a slope for the ratio suggests the power spectrum is not sensitive
to lattice spacing alone. A similar comparison (right) of the energy spectrum
for systems with a constant volume yields a similar ratio. Finally, the energy
propagation was found to be independent of the configuration of the stirring
force (the choice of $k_{f}$, effective $k_{max}$, and $|F|$) provided it is
large enough to induce a sufficient velocity capable of producing turbulence.

Figure 4.6: The energy spectrum of lattice configurations
$N\alpha=(128)(0.16),(256)(0.16),(512)(0.16)$, and $(1024)(0.16)$ (increasing
volume), with a constant filtered average fluid velocity of $0.015$. The left
panel compares the energy spectrum to the expected slope of $-5/3$ at log
scale, and the right is the energy spectrum divided by the expected slope. The
energy spectrum in each panel is scaled for comfortable comparison. The
similarity of the slope of each curve suggests the energy propagation is not
sensitive to volume when the volume is controlled by the number of lattice
nodes. Figure reproduced from [103].

Figure 4.7: Left is the ratio of the power spectrum $E(k)$ for simulations
with lattice spacing $0.32$ and $0.16$ and a constant $N$. The ratio at each
$k$ trends to constant, suggesting the simulated results are not affected by
lattice spacing alone. The $\sim 3/4$ ratio is attributed to an additional
mode available with a tighter lattice. The lattice with the smaller volume is
scaled with a factor of $2^{2}$ to compensate for the larger energy capacity
of the larger volume system. Right is the same power spectrum ratio for
simulations with lattice configurations $N\alpha=(128)(0.032)$ and
$(128)(0.016)$. The ratio also tends to be constant, suggesting the simulated
results are not affected by the non-physical lattice parameters $N$ and
$\alpha$. Figure reproduced from [103].

###  4.4 Conclusions

The energy spectrum of the modeled flow demonstrates agreement with
Kolmogorov’s predicted slope of $-5/3$ in the inertial range indicating an
inverse cascade of energy expected for a turbulent flow. The error is large at
large scales and shows evidence of a pileup, particularly for lattice
configurations with a smaller number of scales in the infrared range. But the
slope is reasonably reproduced in the inertial range above the large scales up
to the driven scales, implying the model is a reasonable reproduction of
turbulent two-dimensional flows of a relativistic fluid. The model is shown to
be insensitive to non-physical parameters including lattice size except with
respect to the effects of a large lattice spacing.

##  Chapter 5 Part II: Reproduction of turbulence in graphene

Graphene is a single atom thick sheet of graphite made up of carbon atoms in
an honeycomb lattice with interesting conductive properties. The linear
dispersion relation of its band structure transports the electric signal
through massless quasi-particles, and the flow of the particles can be
analyzed through a relativistic hydrodynamic mechanism. A sample of graphene
will contain impurities embedded within the lattice structure or within the
supporting substrate that affect the current based on their number density,
size, placement and electric properties. In a hydrodynamic flow, if the
average fluid speed around these obstructions is large enough compared to the
viscous damping effect, a turbulent flow is possible.

In this section we note a dependency of turbulence in the Dirac fluid in
graphene on the position of an obstruction relative to other impurities
present in the sample. The effect is similar to the susceptibility of a wind
farm to turbulence based on the placement of the turbines. An engineer must
take into account the effect the wake a turbine has on downwind turbines when
determining their placement. At high wind velocities a single turbine can
create a vortex street, or regularly sized and spaced vortices, within the
flow of its wake (fig. 5.1 left). If another turbine is placed relatively
close to its upstream neighbor near its wake, the combined effect may produce
a turbulent flow (fig. 5.1 right) which will affect the efficiency of the farm
[104] [105].

Figure 5.1: The left shows a vortex street effect a single wind turbine can
have on the air flow in its wake when the wind velocity is large enough. The
right sketches a turbulent wind flow caused by additional wind turbines
positioned in the wake of another. The vortex street can turn into a turbulent
flow.

Though charge neutral graphene is an exceptionally conductive material and the
charged quasi-particles are able to obtain a large velocity compared to the
viscosity, it is believed that it is unlikely to be capable of producing a
turbulent flow under normal circumstances [106]. However, a recent study by
Mendoza, Herrmann, and Succi [107] modeled a 5 $\mu m$ impurity obstructing a
quasi-particle flow of velocity $10^{5}$ $m/s$ using a two-dimensional
relativistic hydrodynamic numerical model [108]. They calculate the Reynolds
number to be on the order of $100$ and proposed that, under circumstances such
as these, graphene could potentially produce a flow that is in the
preturbulent range. The present report employs a similar hydrodynamic modeler
based on the lattice Boltzmann method to reproduce the flow of a charged Dirac
fluid in an experimentally realistic sample of graphene, incorporating
impurities in the substrate or embedded within the honeycomb lattice. Though
the Reynolds number calculation in this work differs from what was obtained in
[107], we find a signal of a possible preturbulent flow under specific
conditions related to the placement of the impurities.

We employ a hydrodynamic modeler similar to [107] and described in 3.3 to
reproduce the flow of a charged Dirac fluid in an experimentally realistic
sample of graphene, incorporating impurities embedded within the honeycomb
lattice or in the substrate. Using experimentally realistic parameters, though
the estimate for the Reynolds number in this work is smaller, we find a signal
of possible preturbulent flow, corroborating the results of [107]. What is new
in this work is that we also consider _multiple_ impurities in a single sample
and investigate if a relativistic flow around these can give rise to turbulent
(as opposed to preturbulent) flow signals.

###  5.1 Model and Methods

The Reynolds number provides the most commonly accepted metric to predict the
presence of turbulence. Classically the Reynolds number is determined, as
stated, using the kinematic viscosity $\nu$ defined as the ratio of the
dynamic viscosity $\eta$ to the fluid’s mass density $\rho$;
$\nu=\frac{\eta}{\rho}$. For fluid systems consisting of massless particles,
formulations of the Reynolds number commonly replace the mass density with the
entropy density $s$. Therefore, the expression for the dimensionless viscosity
is $\frac{\eta}{s}$, and the Reynolds number is expressed as
$Re=\frac{UL}{(c^{2}/T)(\eta/s)}$, where $\frac{c^{2}}{T}$ balances the units.
The Reynolds number found in [107] is determined with this formulation.
Alternatively, one is able to retain the classical expression of kinematic
viscosity for a fluid of massless particles by defining the mass density in
terms of the number density of the quasi-particles and their “effective” mass;
$\rho=nm_{e}$. The number density is determined experimentally using
techniques such as the measurement of the Hall resistance, while the effective
electron mass is determined through methods such as the Shubnikov-De Haas
effect. The resulting Reynolds number is a ratio incorporating the velocity of
a volume of massless particles with respect to inertial mass to the kinematic
viscosity defined in terms of the effective mass with respect to the
particles’ electric properties. This more classical form of the Reynolds
number is more readily compared to that of traditional fluids.

Turbulence can also be identified through its effect on fluctuations in the
current density. A turbulent flow creates vortices within the fluid creating
multiple traveling modes in Fourier space producing a broadband signal.
Therefore, if the spectrum in frequency space and wavenumber space are both
broadband, the flow can be considered turbulent.

To explore the possibility of a turbulent flow of the charged quasi-particles
in a sample of graphene under realistic conditions we simulate the
hydrodynamic equations of motion using the described adaptation of the
relativistic lattice Boltzmann method defined by Romatschke, Mendoza and Succi
[42] (see 3.3). The RLBM model reproduces a two-dimensional quasi-particle
flow in a $40\mu m$ by $10\mu m$ sample of graphene with one or two rigid
impurities obstructing the flow. The sample is simulated with a $1024\times
256$ node lattice with a spacing of $0.038\mu m$. Initial tests incorporate a
single circular impurity with a diameter $D=0.5\mu m$ placed within the sample
in a region near the inflow boundary, referred to as the “obstruction region”.
Subsequent tests are conducted with a second impurity placed at controlled
distances from the first within the same region. At normal temperatures the
impurities within a graphene sample are believed to be charged, creating
charge puddles in the surrounding region that affect the electric flow. The
impurities are largely sourced from the substrate [109] [110], but can also be
embedded within the sample itself. The size of the impurities can vary greatly
depending on the foreign material, but they typically stay below approximately
$0.5\mu m$ [106]. The size and placement of the impurities are difficult to
control in an experimental setting, but the model seeks to simulate the
effects of one or two quasi-isolated impurities on the current density in
ideal but realistic conditions in order to determine if the detection of a
turbulent signal is possible. Therefore the diameter of the obstacle is chosen
to maximize its turbulence producing potential while maintaining a realistic
size. The velocity of the charged flow in the Dirac liquid can be relatively
high owing to a large effective electric coupling constant. Flow speeds on the
order of $10^{2}$ $m/s$ are common [111] [112], but the flow can approach
velocities as large as 10% of the Fermi velocity, $v_{F}\approx 1.1\times
10^{6}m/s$ [106]. In order to maximize the possibility of a turbulent signal,
the model introduces the largest realizable fluid velocity of $10^{5}$ $m/s$
($0.1$ $v_{F}$) into the sample at the same magnitude along the inflow border.
The borders that are perpendicular to the inflow border use periodic boundary
conditions, effectively simulating an infinitely wide sample with multiple,
regularly placed obstacles, but at a distance where the wakes created by the
obstacles cannot affect each other. Each lattice node in a region occupied by
an impurity implements bounce-back boundary conditions.

The Reynolds number is determined for each test to predict the presence of
turbulence using the kinematic viscosity based on the number density as
described. The kinematic viscosity in graphene is found to be $0.000132$
$\frac{\hbar c^{2}}{eV}$ in natural units ($c=k_{B}=\hbar=1$) [27], readily
obtained from the value of the diffusion constant for momentum
$\eta=2.633\times 10^{3}\frac{eV^{2}}{c^{2}\hbar}$ found in [27], the number
density $n=38.93\frac{eV^{2}}{\hbar^{2}c^{2}}$ determined in [113], and the
effective mass of a charged Dirac quasi-particle participating in the electric
flow, usually given in terms of the mass of an electron $m_{e}$, where
$m_{e}=5.11\times 10^{5}\frac{eV}{c^{2}}$. The effective electron mass is
taken here to be $1.0$ $m_{e}$. The diameter of the embedded impurity is the
most appropriate choice for the system’s characteristic length for the
Reynolds number formulation; $U=D$.

The state data for each node in the lattice is collected throughout the
simulation and the local macroscopic moments are determined and recorded. The
current density $\boldsymbol{j}$ is calculated along the lattice nodes at the
outflow border and the frequency of the fluctuations is determined in Fourier
space against time. The spatial fluctuation of the current density is
determined along the lattice nodes in the region of the sample down stream
from the current inflow referred to as the “current density sampling region”.
The spectrum of the current density fluctuations in frequency space and in
k-space are recorded and plotted for qualitative inspection to look for
evidence of mode generation and migration, indicative of a turbulent flow.

###  5.2 Results

####  5.2.1 Single Impurity

A simulation of a Dirac fluid flow in a sample of graphene containing a single
impurity of diameter $0.5\mu m$ shows a breaking of longitudinal symmetry that
develops into vortices in the wake, and forms a flow pattern known as a von
Kármán vortex street. The vortices form on both sides of the obstacle’s wake
at a roughly consistent size and placement, alternating on either side of the
wake, and shedding at regular intervals (fig. 5.2). The vortex creation,
coherence and shedding produces temporal fluctuations in the current density
that are detected at the outflow border and spatial fluctuations found
throughout the sampling region (fig. 5.3). The fluctuations produce one or
perhaps two prominent modes in the frequency spectrum, but they appear to
broaden to create a slope of about $\omega^{-5/3}$ in small portions of the
spectrum. There is, however, a single prominent mode in the wave number
spectrum indicating the flow around the obstacle does not produce a turbulent
signal in the current density. The Reynolds number for this system is readily
calculated in terms of modified natural units based on the Fermi velocity (see
Appendix B) using the diameter of the impurity as the characteristic length
($L=D=0.5\mu m=3562.5\frac{\hbar v_{F}}{eV}$), the average flow speed
($U=0.1v_{F}$), and the number-density-dependent kinematic viscosity
($\nu=11.8\frac{\hbar v_{F}^{2}}{eV}$). It is found to be

$Re=\frac{LU}{\frac{\eta}{\rho}}=\frac{\left(3562.5\frac{\hbar
v_{F}}{eV}\right)\left(0.1v_{F}\right)}{\left(11.8\frac{\hbar
v_{F}^{2}}{eV}\right)}=6.37.$ (5.1)

It is near the range for preturbulence, but a von Kármán vortex street in the
flow is a preturbulence phenomenon. Therefore we conclude the single impurity
in this system creates a preturbulent flow, but not turbulence.

Figure 5.2: A velocity magnitude heatmap depicting the flow of charge carrying
quasi-particles on a simulated sample of graphene at $0.1v_{F}$ around an
impurity of diameter $0.5\mu m$ in the obstruction region. A large velocity
around the obstacle creates a regular pattern of vortices in the wake
affecting current density fluctuations detected in the sampling region. The
von Kármán vortex street flow pattern is too regular to be considered
turbulent, producing only a single mode in wave number space. Figure 5.3: The
effects of a single impurity of diameter $0.5\mu m$ on the current density
with respect to time $j(t)$ (top left) and space $j(x)$ (bottom left). The
frequency spectrum of $j(\omega)$ (top right in log scale) shows one or two
prominent modes that are beginning to migrate. One prominent mode is seen in
the wave number spectrum of $j(k_{x})$ (bottom right in log scale) indicating
the chaotic flow does not show signs of turbulence. The current density $j$ is
in units of the Fermi velocity $v_{F}$, and $t,x,\omega,$ and $k_{x}$ are in
modified natural units (see Appendix B).

####  5.2.2 Multiple Impurities

A second obstacle is introduced in a subsequent model with the same impurity
size and lattice spacing as in the initial test, $0.5\mu m$ and $0.038\mu m$
respectively. The new obstacle is placed next to, but slightly offset behind
the first with respect to the flow in the obstruction region at a separation
distance of about $3.38\mu m$ (fig. 5.4). The Reynolds number for this system,
found to be $12.74$, is similar to the ratio determined for the initial,
single-obstacle test and falls to the edge of the lower boundary of what can
be considered preturbulent. However, we see a less regular vortex pattern and
the von Kármán vortex street flow pattern is no longer present. The
fluctuations in the current density show emerging modes in both frequency
space and in wave number space which appear to broaden, forming a spectral
slope conforming to the power of $-5/3$ as they migrate (fig. 5.5). There is
one or two slightly prominent outlier modes in frequency space, but the higher
wave number modes are broadband. The spectrum in k-space conforms well to a
$|k|^{-5/3}$ slope that one might expect for mode creation in a nonlinear
system. Comparing the current density fluctuations in fig. 5.2 with fig. 5.4,
we conclude that the resulting flow in this arrangement of obstacles, despite
a smaller Reynolds number, may be considered turbulent.

Figure 5.4: A velocity magnitude heatmap of a quasi-particle flow around two
impurities of size $0.5\mu m$, showing an irregular pattern of vortices in the
wake of the two obstacles suggesting turbulence may be present in the flow.
Vortex creation and shedding is evident as in the same system modeling a
single impurity, but the irregularity of the coherence of the vortex
structures cause a much less regular pattern in the fluctuations in the
current density. Figure 5.5: The current density of a quasi-particle flow
around two impurities of size $0.5\mu m$ with respect to time $j(t)$ (upper
left), with respect to space along the axis parallel to the flow $j(x)$ (lower
left), with respect to frequency $j(\omega)$ (upper right in log scale), and
with respect to wave number $j(k_{x})$ (lower right in log scale). The current
density in frequency space, $j(\omega)$, shows the formation of multiple modes
propagating from a single mode creating a rough conformity to a $-5/3$ slope
in log scales. The current density in k-space, $j(k_{x})$, is broadband with a
close adherence to a $|k|^{5/3}$ slope expected for migrating modes in a non-
linear system. The broadband spectra in frequency space and in k-space hint at
a potential turbulent signal. The current density $j$ is in units of the Fermi
velocity $v_{F}$, and $t,x,\omega,$ and $k_{x}$ are in modified natural units
(see Appendix B).

We find the presence of turbulence in a two-dimensional solid such as graphene
to be sensitive to obstacle placement. A series of subsequent tests place the
second obstacle at various other distances from the first in directions both
parallel and perpendicular with respect to the inflow. The second obstacle is
positioned such that the flow on the front side of the obstacle is impacted to
some degree by the wake created by the first. When the trailing obstacle is
positioned within approximately $2\mu m$ of the leading obstacle, no turbulent
or preturbulent signal is detected in $j$. At this range the current density’s
wave number spectrum shows a single dominant mode, and the frequency spectrum
shows a few distinct prominent modes with initial signs of broadening, but is
not broadband (fig. 5.6). The spectra are very similar to those created by a
single obstacle (see fig. 5.3), implying the close proximity of the obstacles
has the same effect on the current density as a single, larger obstacle.
Additionally, there is no evidence of a turbulent signal when the obstacles
are situated at a large distance with respect to the characteristic length of
the system. For models examined in this work, a turbulent signal is not
detected when the obstacles are separated at distances greater than $5\mu m$.
The current density fluctuations recorded in a simulation of a model with
obstacles $6.02\mu m$ apart show a single prominent mode in wave number space
but a somewhat broadband spectrum in the high frequency range in frequency
space, with the exception of a single outlier mode; also similar to the
current density fluctuations caused by a single obstacle (fig. 5.7). The
contribution of the non-linear effects caused by the wake of first obstacle
evidently dissipates at larger distances.

Figure 5.6: Current density fluctuations recorded from a modeled flow
obstructed by two impurities positioned relatively close. The proximity
produces a single dominant mode in wave number space (shown in log scale)
similar to what is seen in $j(k_{x})$ from a single obstacle (fig. 5.3),
implying the effect of closely placed obstacles on the current density is
similar to that caused by a single obstacle. The multi-modal spectrum in
frequency space (shown in log scale) is also similar the $j(\omega)$ spectrum
for a single obstacle. The units of the plot are indicated in Appendix B.
Figure 5.7: Current density fluctuations from a model simulating a flow around
two impurities situated at a large distance relative to their diameter. In
this arrangement two discrete modes appear prominent in the wave number
spectrum of $j$ shown in log scale, and multiple modes are visible in the
frequency spectrum of $j$, also in log scale. A turbulent signal is not
detected for obstacles placed too far apart. The units of the plot are
indicated in Appendix B.

Fig. 5.8 depicts the current density in frequency space and in k-space of five
runs with increasing obstacle placement. At the smallest distance, $1.68\mu
m$, the frequency space spectrum is multi-modal while the k-space spectrum has
a single prominent mode, resembling the spectral effects of a single obstacle.
At a slightly larger distance, $2.34\mu m$, the k-space spectrum becomes
broadband and conforms to the slope of $|k|^{-5/3}$. More modes are present in
the frequency spectrum, but it is not definitively broadband. A broadband
signal is present in both frequency space and wave number space for objects
placed at a distance of $3.38\mu m$ (also shown in fig. 5.4), indicating
turbulence. As the distance increases from $3.38\mu m$ to $6.02\mu m$, the
frequency spectrum reverts back to multiple modes, and the broadband spectrum
in k-space dissipates into a single mode. The turbulent signal is gone when
the objects are separated at this distance. At the largest separation tested,
$8.39\mu m$, a small number of prominent modes are present in both spectra of
$\bm{j}$ so that it also resembles the current density spectra created by a
single obstacle. Further tests investigating different placement
configurations (not shown) indicate that detection of turbulence is sensitive
to other positional features such as alignment. Because one obstacle must be
within the influence of the other’s wake, when the objects are adjacent and
too far apart a turbulent signal is not detected.

Figure 5.8: The effect of object separation of two impurities on the current
density spectra. A turbulent signal is detected when the spectrum in both
frequency space ($j(\omega)$ in log scale in the second column) and wave
number space ($j(k_{x})$ in the third column; log scale) are broadband. There
is no turbulent signal when the impurities are less than approximately $2\mu
m$ apart as evidenced by the lack of a broadband signal in frequency space. A
broadband signal is similarly absent when the obstacles are separated by
greater than approximately $5\mu m$.

Note that the identification of a broadband spectrum contains an element of
judgment, and a determination for the presence of turbulence using the
spectrum is not rigorous. Because each model configuration shows multiple
modes in frequency space, and because the most prominent mode for each in
k-space could be considered an outlier, it is possible that each depicted
obstacle configuration can be considered capable of producing a turbulent
signal. Increasing the number of lattice nodes in the model will produce a
higher resolution spectrum and subsequently more confidence in the
interpretation, but it would also be more computationally expensive.

###  5.3 Conclusions

To summarize, we were able to detect a potential turbulent flow in the
relativistic hydrodynamics of massless charged quasi-particles in an ultra
pure, idealized sample of graphene. Though the Reynolds number of the modeled
systems were on the lower edge of the preturbulent range, evidence of
turbulence emerges as a result of interactions in the wake of multiple
obstacles placed in the sample. For two similarly sized impurities a signal of
turbulence is dependent upon their separation and position, while apparently
independent of size or shape. This dependence creates a more complex
formulation of turbulence that cannot be represented by the Reynolds number
metric alone, and has consequences for the conductivity of graphene. The
sensitivity to turbulence based on impurity placement is similar to the
sensitivity of the placement of wind turbines in a wind farm. A notable
difference, though, is the importance of higher altitudes above the turbines
on the air flow through the farm. The electric flow within two-dimensional
graphene is not affected by out of plane effects.

A preturbulent signal was detected in a test modeled after [107], in agreement
with the findings of that work. Further investigation is needed to study the
effect of semi-rigid obstacles and momentum relaxation on the production of a
turbulent signal. The resonance of a semi-rigid obstacle caused by the flow
would amplify the current density fluctuations, and the effect would likely be
multiplied by the presence of additional semi-rigid obstacles. However, the
effect would be balanced or mitigated by momentum relaxation effects caused by
the obstacle and the lattice itself. The lattice structure of another two-
dimensional solid such as a Kagome metal should also be investigated with the
employed RLBM numerical modeler. The two-dimensional molecular lattice of a
Kagome metal creates a linear dispersion relation at low energies, similar to
that of graphene, but it has a much stronger electric coupling constant that
enables a faster Dirac fluid flow [114]. With a similar viscosity, this
implies a higher likelihood of detecting a preturbulent or turbulent flow.

##  Chapter 6 Summary and Conclusions

In Chapter 4 we explored the properties of a massless fluid system
incorporating a stirring force into a relativistic modification of the Lattice
Boltzmann Method, driving the simulation of a two-dimensional fluid system at
small spectral scales. The energy spectrum of the modeled flow demonstrated
the slope expected for an inverse cascade of energy propagation. The detection
of an energy cascade lends confidence that the RLBM modeler, along with the
addition of a forcing term, is capable of a reasonable reproduction of
turbulent flows for a fluid of massless particles.

In Chapter 5 we found that a turbulent flow in a sample of graphene in a
realistic but ideal setting is possible, although it requires the interaction
between multiple impurities and is dependent upon their separation and
position. Though, it is apparently only weakly dependent on size or shape at
reasonable scales. We were able to reproduce the preturbulent flow of massless
charge carriers within graphene demonstrated by [107], but were able to
produce a fully turbulent flow with an additional impurity incorporated within
the model at specific separation distances.

The sensitivity of turbulence to impurity placement is similar to the
placement sensitivity of wind turbines in a wind farm. An engineer must take
into account the effect the wakes created by each turbine have on their
neighbors when designing the arrangement of the turbines in order to minimize
wind turbulence that will disrupt the farm’s efficiency.

The RLBM modeler is readily applied to other two-dimensional solid systems
exhibiting a low energy linear dispersion relation such as a Kagome solid. In
a Kagome solid the massless fluid flow is believed to have a similar viscosity
to graphene, but can reach larger velocities owing to a stronger electric
coupling constant [114], and is also experimentally accessible. Further
applications of the RLBM in the exploration of turbulent flows should continue
to be investigated. An investigation of the energy cascade of a three
dimensional system will help validate the model’s ability to reproduce
turbulence. Different turbulence-generating forcing schemes should be
considered in order to achieve higher velocities.

The effect of semi-rigid obstacles and momentum relaxation on the production
of a turbulent signal can be investigated with a modest alteration of the
boundary handling in the RLBM code. The resonance of a semi-rigid obstacle
caused by the flow would amplify the current density fluctuations, and the
effect would likely be multiplied by the presence of additional semi-rigid
obstacles. However, the effect would be balanced or mitigated by momentum
relaxation effects caused by the obstacle and the lattice itself.

Further examination of the turbulent flow created by two or more impurities at
transitional Reynolds number ranges should be carried out to determine, in
better detail, the relationship between obstacle placement and the resulting
turbulence. This might include an objective quantification of the presence of
a broadband signal in the spectra of the current density and how it relates to
the position of a second particle in two dimensions. The approach should then
be extended to three dimensions. The vortex dynamics involved in a turbulent
flow differs significantly between systems in two and three dimensions, thus
the relationship can be expected to be very different, if it indeed exists.

Having successfully reproduced the turbulent flow of graphene’s charged quasi-
particles, this formulation of the Relativistic Lattice Boltzmann Method can
be assumed to reasonably reproduce turbulence in relativistic flows that are
less accessible to experiment. It therefore may be used, for example, to
describe the group dynamics of the particles on the outside of the event
horizon of a black hole, which are thought to make up a relativistic fluid. It
can also be applied to other relativistic flows like the quark-gluon plasma
that is created just after a high energy particle collision, or the
relativistic jets created by compact astrophysical objects such as a neutron
star or black hole.

## Bibliography

  * [1] Hans Bantilan, Takaaki Ishii, and Paul Romatschke.  Holographic heavy-ion collisions: Analytic solutions with longitudinal flow, elliptic flow and vorticity.  Physics Letters B, 785:201–206, oct 2018.
  * [2] Luciano Rezzolla and Olindo Zanotti.  Relativistic Hydrodynamics.  Oxford University Press, 2013.
  * [3] David Radice and Luciano Rezzolla.  Universality and intermittency in relativistic turbulent flows of a hot plasma.  The Astrophysical Journal Letters, 766(1):L10, 2013.
  * [4] Paul Romatschke and Ulrike Romatschke.  Viscosity information from relativistic nuclear collisions: How perfect is the fluid observed at RHIC?  Physical Review Letters, 99(17), oct 2007.
  * [5] K. Dusling and D. Teaney.  Simulating elliptic flow with viscous hydrodynamics.  Physical Review C, 77(3), mar 2008.
  * [6] Björn Schenke, Sangyong Jeon, and Charles Gale.  Elliptic and triangular flow in event-by-event D = 3 + 1 viscous hydrodynamics.  Physical Review Letters, 106(4), jan 2011.
  * [7] Huichao Song, Steffen A. Bass, and Ulrich Heinz.  Spectra and elliptic flow for identified hadrons in $2.76$A TeV Pb + Pb collisions.  Physical Review C, 89(3), mar 2014.
  * [8] Jonah E. Bernhard, J. Scott Moreland, and Steffen A. Bass.  Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma.  Nature Physics, 15(11), nov 2019.
  * [9] F. Mohseni, M. Mendoza, S. Succi, and H. J. Herrmann.  Relativistic effects on the richtmyer-meshkov instability.  Physical Review D, 90(12), dec 2014.
  * [10] Ryan D Weller and Paul Romatschke.  One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+ p, p+ pb and pb+ pb collisions at s= 5.02 tev.  Physics Letters B, 774:351–356, 2017.
  * [11] Pasi Huovinen, Peter F Kolb, Ulrich Heinz, PV Ruuskanen, and Sergei A Voloshin.  Radial and elliptic flow at rhic: Further predictions.  Physics Letters B, 503(1-2):58–64, 2001.
  * [12] Matthew Luzum and Paul Romatschke.  Viscous hydrodynamic predictions for nuclear collisions at the lhc.  Physical review letters, 103(26):262302, 2009.
  * [13] M Habich, JL Nagle, and P Romatschke.  Particle spectra and HBT radii for simulated central nuclear collisions of C + C, Al + Al, Cu + Cu, Au + Au, and Pb + Pb from $\sqrt{s}=62.4$ $s=62.4$ – $27602760$ GeV.  The European Physical Journal C, 75(1):1–6, 2015.
  * [14] Paul Romatschke.  Light-heavy-ion collisions: a window into pre-equilibrium qcd dynamics?  The European Physical Journal C, 75(7):1–15, 2015.
  * [15] JL Nagle, A Adare, S Beckman, T Koblesky, J Orjuela Koop, D McGlinchey, P Romatschke, J Carlson, JE Lynn, and M McCumber.  Exploiting intrinsic triangular geometry in relativistic he 3+ au collisions to disentangle medium properties.  Physical Review Letters, 113(11):112301, 2014.
  * [16] O. Furtmaier, M. Mendoza, I. Karlin, S. Succi, and H. J. Herrmann.  Rayleigh-Bénard instability in graphene.  Physical Review B: Solid State, 91(8):085401, Feb 2015.
  * [17] A. Gabbana, M. Polini, S. Succi, R. Tripiccione, and F. M. D. Pellegrino.  Prospects for the detection of electronic preturbulence in graphene.  Physical Review Letters, 121(23):236602, Dec 2018.
  * [18] Markus Müller, Jörg Schmalian, and Lars Fritz.  Graphene: A nearly perfect fluid.  Physical review letters, 103(2):025301, 2009.
  * [19] Lars Fritz, Jörg Schmalian, Markus Müller, and Subir Sachdev.  Quantum critical transport in clean graphene.  Physical Review B, 78(8):085416, 2008.
  * [20] Yuanbo Zhang, Yan-Wen Tan, Horst L. Stormer, and Philip Kim.  Experimental observation of the quantum hall effect and berry’s phase in graphene.  Nature, 438:201–204, January 2005.
  * [21] R. Krishna Kumar, D. A. Bandurin, F. M. D. Pellegrino, Y. Cao, A. Principi, H. Guo, G. H. Auton, M. Ben Shalom, L. A. Ponomarenko, G. Falkovich, K. Watanabe, T. Taniguchi, I. V. Grigorieva, L. S. Levitov, M. Polini, and A. K. Geim.  Superballistic flow of viscous electron fluid through graphene constrictions.  Nature Physics, 13:1182–1185, January 2017.
  * [22] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov.  Two-dimensional gas of massless dirac fermions in graphene.  Nature, 438:197–200, Nov 2005.
  * [23] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov.  Two-dimensional gas of massless Dirac fermions in graphene.  Nature, 438(7065):197–200, Nov 2005.
  * [24] Sybren Ruurds De Groot.  On relativistic kinetic theory.  In Fundamental problems in statistical mechanics. Volume 5. 1980\.
  * [25] Carlo Cercignani and Gilberto Medeiros Kremer.  Relativistic boltzmann equation.  In The Relativistic Boltzmann Equation: Theory and Applications, pages 31–63. Springer, 2002.
  * [26] A Aliano, L Rondoni, and GP Morriss.  Maxwell-jüttner distributions in relativistic molecular dynamics.  The European Physical Journal B-Condensed Matter and Complex Systems, 50(1):361–365, 2006.
  * [27] Markus Müller, Jörg Schmalian, and Lars Fritz.  Graphene: A nearly perfect fluid.  Phys. Rev. Lett., 103:025301, Jul 2009.
  * [28] VG Bar’yakhtar, ED Belokolos, ON Golubeva, and AD Sukhanov.  A ratio of the shear viscosity to the density of entropy for helium.  Ukrainian Journal of Physics, 59(2):155–155, 2014.
  * [29] Matthew Luzum and Jean-Yves Ollitrault.  Extracting the shear viscosity of the quark-gluon plasma from flow in ultra-central heavy-ion collisions.  Nuclear Physics A, 904:377–380, May 2013.
  * [30] R. Willett, J. P. Eisenstein, H. L. Störmer, D. C. Tsui, A. C. Gossard, and J. H. English.  Observation of an even-denominator quantum number in the fractional quantum hall effect.  Phys. Rev. Lett., 59:1776–1779, Oct 1987.
  * [31] Domenico Di Sante, Johanna Erdmenger, Martin Greiter, Ioannis Matthaiakakis, Rene Meyer, David Rodriguez Fernand ez, Ronny Thomale, Erik van Loon, and Tim Wehling.  Turbulent hydrodynamics in strongly correlated Kagome metals.  Nature Communications, 11:3997, 2020.
  * [32] P. L. Bhatnagar, E. P. Gross, and M. Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  Phys. Rev., 94:511–525, May 1954.
  * [33] Shiyi Chen and Gary D. Doolen.  Lattice boltzmann method for fluid flows.  Annual Review of Fluid Mechanics, 30(1):329–364, 1998.
  * [34] Sauro Succi.  The lattice Boltzmann equation: for fluid dynamics and beyond.  Oxford university press, 2001.
  * [35] Raúl Machado.  On the generalized hermite-based lattice boltzmann construction, lattice sets, weights, moments, distribution functions and high-order models.  Frontiers of Physics, 9(4):490–510, 2014.
  * [36] Xiaowen Shan and Xiaoyi He.  Discretization of the velocity space in the solution of the boltzmann equation.  Physical Review Letters, 80(1):65, 1998.
  * [37] FJ Higuera, S Succi, and R Benzi.  Lattice gas dynamics with enhanced collisions.  EPL (Europhysics Letters), 9(4):345, 1989.
  * [38] M. Mendoza, B. M. Boghosian, H. J. Herrmann, and S. Succi.  Derivation of the lattice Boltzmann model for relativistic hydrodynamics.  Physical Review D, 82(10):105008, Nov 2010.
  * [39] F. Mohseni, M. Mendoza, S. Succi, and H. J. Herrmann.  Lattice Boltzmann model for ultrarelativistic flows.  Physical Review D, 87(8):083003, April 2013.
  * [40] M. Mendoza, I. Karlin, S. Succi, and H. J. Herrmann.  Relativistic lattice Boltzmann model with improved dissipation.  Physical Review D, 87(6):065027, March 2013.
  * [41] A. Gabbana, D. Simeoni, S. Succi, and R. Tripiccione.  Relativistic lattice Boltzmann methods: Theory and applications.  Physics Reports, 863:1–63, June 2020.
  * [42] P. Romatschke, M. Mendoza, and S. Succi.  Fully relativistic lattice Boltzmann algorithm.  Physical Review C, 84(3):034903, Sep 2011.
  * [43] Victor E. Ambrus and Robert Blaga.  High-order quadrature-based lattice Boltzmann models for the flow of ultrarelativistic rarefied gases.  Physical Review C, 98(3):035201, September 2018.
  * [44] Doug McLean.  Continuum fluid mechanics and the navier-stokes equations.  Understanding Aerodynamics: Arguing from the Real Physics, pages 13–78, 2012.
  * [45] Paul Romatschke and Ulrike Romatschke.  Relativistic fluid dynamics in and out of equilibrium: and applications to relativistic nuclear collisions.  2019\.
  * [46] Holger Babinsky.  How do wings work?  Physics Education, 38(6):497, 2003.
  * [47] L.D. Landau, E.M. Lifšic, J.B. Sykes, and W.H. Reid.  Course of Theoretical Physics.  Teoretičeskaja fizika. Butterworth-Heinenann, 1987.
  * [48] Michael J Moran, Howard N Shapiro, Daisie D Boettner, and Margaret B Bailey.  Fundamentals of engineering thermodynamics.  John Wiley & Sons, 2010.
  * [49] Jasmine Brewer, Miller Mendoza, Ryan E. Young, and Paul Romatschke.  Lattice Boltzmann simulations of a strongly interacting two-dimensional Fermi gas.  Physical Review A, 93(1):013618, Jan 2016.
  * [50] Alexander Migdal.  Universal area law in turbulence.  arXiv preprint arXiv:1903.08613, 2019.
  * [51] Geoffrey Ingram Taylor.  The spectrum of turbulence.  Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences, 164(919):476–490, 1938.
  * [52] HA Rose and PL Sulem.  Fully developed turbulence and statistical mechanics.  Journal de Physique, 39(5):441–484, 1978.
  * [53] Arjun Berera and Richard DJG Ho.  Chaotic properties of a turbulent isotropic fluid.  Physical review letters, 120(2):024101, 2018.
  * [54] Uriel Frisch.  Turbulence: The Legacy of A. N. Kolmogorov.  Cambridge University Press, 1995.
  * [55] Andrei Nikolaevich Kolmogorov, V. Levin, Julian Charles Roland Hunt, Owen Martin Phillips, and David Williams.  The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers.  Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, 434(1890):9–13, 1991.
  * [56] Robert H. Kraichnan.  Inertial Ranges in Two-Dimensional Turbulence.  Physics of Fluids, 10(7):1417–1423, Jul 1967.
  * [57] T. Dubos, A. Babiano, J. Paret, and P. Tabeling.  Intermittency and coherent structures in the two-dimensional inverse energy cascade: Comparing numerical and laboratory experiments.  Physical Review E, 64(3):036302, Sep 2001.
  * [58] Pablo Dmitruk and David C. Montgomery.  Numerical study of the decay of enstrophy in a two-dimensional navier–stokes fluid in the limit of very small viscosities.  Physics of Fluids, 17(3):035114, 2005.
  * [59] Philip E. Merilees and Helen Warn.  On energy and enstrophy exchanges in two-dimensional non-divergent flow.  Journal of Fluid Mechanics, 69(4):625–630, 1975.
  * [60] RK Scott.  Nonrobustness of the two-dimensional turbulent inverse cascade.  Physical Review E, 75(4):046301, 2007.
  * [61] Chuong V Tran and John C Bowman.  On the dual cascade in two-dimensional turbulence.  Physica D: Nonlinear Phenomena, 176(3-4):242–255, 2003.
  * [62] JR Chasnov.  On the decay of two-dimensional homogeneous turbulence.  Physics of Fluids, 9(1):171–180, 1997.
  * [63] Chuong V Tran and Theodore G Shepherd.  Constraints on the spectral distribution of energy and enstrophy dissipation in forced two-dimensional turbulence.  Physica D: Nonlinear Phenomena, 165(3-4):199–212, 2002.
  * [64] George K Batchelor.  Computation of the energy spectrum in homogeneous two-dimensional turbulence.  The Physics of Fluids, 12(12):II–233, 1969.
  * [65] Ragnar Fjørtoft.  On the changes in the spectral distribution of kinetic energy for twodimensional, nondivergent flow.  Tellus, 5(3):225–230, 1953.
  * [66] Patrick Tabeling.  Two-dimensional turbulence: a physicist approach.  Physics reports, 362(1):1–62, 2002.
  * [67] Annalisa Bracco and James C McWilliams.  Reynolds-number dependency in homogeneous, stationary two-dimensional turbulence.  Journal of fluid mechanics, 646:517–526, 2010.
  * [68] PATRICK F. CUMMINS and GREG HOLLOWAY.  Reynolds stress and eddy viscosity in direct numerical simulations of sheared two-dimensional turbulence.  Journal of Fluid Mechanics, 657:394–412, 2010.
  * [69] D Biskamp, E Schwarz, and A Celani.  Nonlocal bottleneck effect in two-dimensional turbulence.  Physical review letters, 81(22):4855, 1998.
  * [70] Ch H Bruneau and H Kellay.  Experiments and direct numerical simulations of two-dimensional turbulence.  Physical Review E, 71(4):046305, 2005.
  * [71] Robert H Kraichnan and David Montgomery.  Two-dimensional turbulence.  Reports on Progress in Physics, 43(5):547, 1980.
  * [72] Guido Boffetta and Robert E Ecke.  Two-dimensional turbulence.  Annual review of fluid mechanics, 44:427–451, 2012.
  * [73] Mark Hoefer.  personal communication.
  * [74] P. R. Wallace.  The band theory of graphite.  Phys. Rev., 71:622–634, May 1947.
  * [75] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim.  The electronic properties of graphene.  Rev. Mod. Phys., 81:109–162, Jan 2009.
  * [76] Subir Sachdev.  The Landscape of the Hubbard Model, pages 559–620.
  * [77] Qinlong Luo.  The elementary electronic properties of graphene.  Department of Physics and Astronomy, University of Tennessee, 2010\.
  * [78] Bo Zhen, Chia Wei Hsu, Yuichi Igarashi, Ling Lu, Ido Kaminer, Adi Pick, Song-Liang Chua, John D Joannopoulos, and Marin Soljačić.  Spawning rings of exceptional points out of dirac cones.  Nature, 525(7569):354–358, 2015.
  * [79] Khuong K Huynh, Yoichi Tanabe, and Katsumi Tanigaki.  Both electron and hole dirac cone states in ba (feas) 2 confirmed by magnetoresistance.  Physical review letters, 106(21):217004, 2011.
  * [80] Sayanti Samaddar, Jeff Strasdas, Kevin Janßen, Sven Just, Tjorven Johnsen, Zhenxing Wang, Burkay Uzlu, Sha Li, Daniel Neumaier, Marcus Liebmann, et al.  Evidence for local spots of viscous electron flow in graphene at moderate mobility.  Nano Letters, 21(22):9365–9373, 2021.
  * [81] Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo Silva, and Erlend Magnus Viggen.  The lattice Boltzmann method: Principles and practice.  Graduate Texts in Physics. Springer International Publishing, 2017.
  * [82] Takaji Inamuro, Masato Yoshino, and Fumimaru Ogino.  A non-slip boundary condition for lattice boltzmann simulations.  Physics of fluids, 7(12):2928–2930, 1995.
  * [83] S Succi, G Bella, and F Papetti.  Lattice kinetic theory for numerical combustion.  Journal of Scientific Computing, 12(4):395–408, 1997.
  * [84] Guillaume Tauzin, Luca Biferale, Mauro Sbragaglia, Abhineet Gupta, Federico Toschi, Andreas Bartel, and Matthias Ehrhardt.  A numerical tool for the study of the hydrodynamic recovery of the lattice boltzmann method.  Computers and Fluids, 172:241–250, 2018.
  * [85] Gauthier Wissocq, Pierre Sagaut, and Jean-François Boussuge.  An extended spectral analysis of the lattice boltzmann method: modal interactions and stability issues.  Journal of Computational Physics, 380:311–333, 2019.
  * [86] Sydney Chapman and Thomas George Cowling.  The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases.  Cambridge university press, 1990.
  * [87] Kallash Kumar.  The chapman? enskog solution of the boltzmann equation: A reformulation in terms of irreducible tensors and matrices.  Australian Journal of Physics, 20(3):205–252, 1967.
  * [88] Stegun Abramowitz.  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover, New York, ninth dover printing, tenth gpo printing edition, 1964\.
  * [89] Milton Abramowitz and Irene A Stegun.  Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55.  US Government printing office, 1948.
  * [90] Shaozhong Deng.  Quadrature formulas in two dimensions.  University of North Carolina-Charlotte Math, 5172, 2010.
  * [91] Paul Romatschke.  Relativistic (lattice) boltzmann equation with nonideal equation of state.  Physical Review D, 85(6):065012, 2012.
  * [92] LR Weih, A Gabbana, D Simeoni, L Rezzolla, S Succi, and R Tripiccione.  Beyond moments: relativistic lattice boltzmann methods for radiative transport in computational astrophysics.  Monthly Notices of the Royal Astronomical Society, 498(3):3374–3394, 2020.
  * [93] Miller Mendoza, BM Boghosian, Hans Jürgen Herrmann, and S Succi.  Derivation of the lattice boltzmann model for relativistic hydrodynamics.  Physical Review D, 82(10):105008, 2010.
  * [94] M Mendoza, BM Boghosian, Hans Jürgen Herrmann, and S Succi.  Fast lattice boltzmann solver for relativistic hydrodynamics.  Physical review letters, 105(1):014502, 2010.
  * [95] Robert Blaga and Victor E Ambruş.  Quadrature-based lattice boltzmann model for relativistic flows.  In AIP Conference Proceedings, volume 1796, page 020010. AIP Publishing LLC, 2017.
  * [96] Q Li, KH Luo, and XJ Li.  Lattice boltzmann method for relativistic hydrodynamics: Issues on conservation law of particle number and discontinuities.  Physical Review D, 86(8):085044, 2012.
  * [97] Alessandro Gabbana, M Mendoza, S Succi, and Raffaele Tripiccione.  Towards a unified lattice kinetic scheme for relativistic hydrodynamics.  Physical Review E, 95(5):053304, 2017.
  * [98] D Hupp, M Mendoza, I Bouras, S Succi, and HJ Herrmann.  Relativistic lattice boltzmann method for quark-gluon plasma simulations.  Physical Review D, 84(12):125015, 2011.
  * [99] David Oettinger, Miller Mendoza, and Hans J Herrmann.  Gaussian quadrature and lattice discretization of the fermi-dirac distribution for graphene.  Physical Review E, 88(1):013302, 2013.
  * [100] Lorenzo Bazzanini, Alessandro Gabbana, Daniele Simeoni, Sauro Succi, and Raffaele Tripiccione.  A lattice boltzmann method for relativistic rarefied flows in (2+ 1) dimensions.  Journal of Computational Science, 51:101320, 2021.
  * [101] Carlo Cercignani and Gilberto Medeiros Kremer.  The relativistic Boltzmann equation: theory and applications.  Editorial Universitaria de la Republica de Cuba, 2002.
  * [102] Paul Romatschke and Ulrike Romatschke.  Relativistic Fluid Dynamics In and Out of Equilibrium – Ten Years of Progress in Theory and Numerical Simulations of Nuclear Collisions.  arXiv e-prints, page arXiv:1712.05815, Dec 2017.
  * [103] Mark Watson.  Two-dimensional turbulence in a massless fluid with a relativistic lattice boltzmann model.  AIP Advances, 11(3):035333, 2021.
  * [104] Richard JAM Stevens and Charles Meneveau.  Flow structure and turbulence in wind farms.  Annual review of fluid mechanics, 49:311–339, 2017.
  * [105] Tobias Ahsbahs, Nicolai Gayle Nygaard, Alexander Newcombe, and Merete Badger.  Wind farm wakes from sar and doppler radar.  Remote Sensing, 12(3):462, 2020.
  * [106] Andrew Lucas and Kin Chung Fong.  Hydrodynamics of electrons in graphene.  Journal of Physics: Condensed Matter, 30(5):053001, Jan 2018.
  * [107] M. Mendoza, H. J. Herrmann, and S. Succi.  Preturbulent regimes in graphene flow.  Phys. Rev. Lett., 106:156601, Apr 2011.
  * [108] Alessandro Gabbana, Miller Mendoza, Sauro Succi, and Raffaele Tripiccione.  Numerical evidence of electron hydrodynamic whirlpools in graphene samples.  Computers & Fluids, 172:644–650, 2018.
  * [109] Kentaro Nomura and A. H. MacDonald.  Quantum transport of massless dirac fermions.  Phys. Rev. Lett., 98:076602, Feb 2007.
  * [110] Shaffique Adam, E. H. Hwang, V. M. Galitski, and S. Das Sarma.  A self-consistent theory for graphene transport.  Proceedings of the National Academy of Sciences, 104(47):18392–18397, 2007.
  * [111] Inanc Meric, Melinda Y. Han, Andrea F. Young, Barbaros Ozyilmaz, Philip Kim, and Kenneth L. Shepard.  Current saturation in zero-bandgap, top-gated graphene field-effect transistors.  Nature Nanotechnology, 3(11):654–659, 2008.
  * [112] Vincent E. Dorgan, Myung-Ho Bae, and Eric Pop.  Mobility and saturation velocity in graphene on sio2.  Applied Physics Letters, 97(8):082112, 2010.
  * [113] Daniel E. Sheehy and Jörg Schmalian.  Quantum critical scaling in graphene.  Phys. Rev. Lett., 99:226803, Nov 2007.
  * [114] Domenico Di Sante, Johanna Erdmenger, Martin Greiter, Ioannis Matthaiakakis, Rene Meyer, David Rodriguez Fernand ez, Ronny Thomale, Erik van Loon, and Tim Wehling.  Turbulent hydrodynamics in strongly correlated Kagome metals.  arXiv e-prints, page arXiv:1911.06810, Nov 2019.
  * [115] J. C. Maxwell M.A.  V. illustrations of the dynamical theory of gases.—part i. on the motions and collisions of perfectly elastic spheres.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 19(124):19–32, 1860.
  * [116] Paul Romatschke.  New developments in relativistic viscous hydrodynamics.  International Journal of Modern Physics E, 19(01):1–53, 2010.

##  Appendix A Maxwellian Form of the Equilibrium Distribution Function

The form of the equilibrium probability distribution function was determined
by Maxwell using logic described in [115] and is summarized as follows. When
two particles collide (not point particles) their resulting trajectories are
highly sensitive to the relative angles of collision, so that when considering
a large number of particles, the number of collisions results in angular
distribution of the particle directions. After a time, the particles within
the fluid will reach an equilibrium so that, on average, the velocities of the
particles in the fluid match the mean fluid velocity $(\bm{v}=\bm{u})$. That
is, the fluid is isotropic in velocity space in a reference frame moving with
the speed of $\bm{u}$. Then, the equilibrium distribution function is

$f_{eq}(t,|v|,x)$

For this line of logic, we focus specifically on the distribution of the
velocity parameter. We assume that the distribution function is separable in
velocity.

$f_{eq}(|v|^{2})=f_{eq}(v_{x}^{2}+v_{y}^{2}+v_{z}^{2})=f_{eq,x}(v_{x}^{2})f_{eq,y}(v_{y}^{2})f_{eq,z}(v_{z}^{2})$

At equilibrium, the velocity is isotropic and constant.

$f_{eq}(|v|^{2})=C$

That being the case, it is also true that $ln(f_{eq}(|v|^{2}))=ln(C)$ and
$ln(C)=C$. Then

$ln(f_{eq,x}(v_{x}^{2})f_{eq,y}(v_{y}^{2})f_{eq,z}(v_{z}^{2}))=C$

and

$ln(f_{eq,x}(v_{x}^{2}))+ln(f_{eq,y}(v_{y}^{2}))+ln(f_{eq,z}(v_{z}^{2}))=C$
(A.1)

This equation is fulfilled if each separated one-dimensional equilibrium
distribution is of the form $ln(f_{eq,x}(v_{x}^{2}))=a+bv_{x}^{2}$, where $a$
and $b$ are both constants. Plugging that in to A.1 we have

$ln(f_{eq,x}(v_{x}^{2}))+ln(f_{eq,y}(v_{y}^{2}))+ln(f_{eq,z}(v_{z}^{2}))=3a+b\left(v_{x}^{2}+v_{y}^{2}+v_{z}^{2}\right).$

This works because the right-hand side is constrained to be constant, and $a$,
$b$ and $|v|=v_{x}^{2}+v_{y}^{2}+v_{z}^{2}$ are all identified as constant.
Then, the form of the equilibrium distribution function must be