# The stochastic relativistic advection diffusion equation from the Metropolis
algorithm

Gökçe Başar<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
North Carolina, Chapel Hill, North Carolina 27599, USA  Jay Bhambure
<EMAIL_ADDRESS>Center for Nuclear Theory, Department of Physics
and Astronomy, Stony Brook University, Stony Brook, New York, 11794-3800, USA
Rajeev Singh<EMAIL_ADDRESS>Center for Nuclear Theory, Department
of Physics and Astronomy, Stony Brook University, Stony Brook, New York,
11794-3800, USA  Department of Modern Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  Derek Teaney
<EMAIL_ADDRESS>Center for Nuclear Theory, Department of Physics
and Astronomy, Stony Brook University, Stony Brook, New York, 11794-3800, USA

###### Abstract

We study an approach to simulating the stochastic relativistic advection-
diffusion equation based on the Metropolis algorithm. We show that the
dissipative dynamics of the boosted fluctuating fluid can be simulated by
making random transfers of charge between fluid cells, interspersed with ideal
hydrodynamic time steps. The random charge transfers are accepted or rejected
in a Metropolis step using the entropy as a statistical weight. This procedure
reproduces the expected strains of dissipative relativistic hydrodynamics in a
specific (and non-covariant) hydrodynamic frame known as the _density frame_.
Numerical results, both with and without noise, are presented and compared to
relativistic kinetics and analytical expectations. An all order resummation of
the density frame gradient expansion reproduces the covariant dynamics in a
specific model. In contrast to all other numerical approaches to relativistic
dissipative fluids, the dissipative fluid formalism presented here is strictly
first order in gradients and has no non-hydrodynamic modes. The physical
naturalness and simplicity of the Metropolis algorithm, together with its
convergence properties, make it a promising tool for simulating stochastic
relativistic fluids in heavy ion collisions and for critical phenomena in the
relativistic domain.

relativistic advection-diffusion equation, stochastic dynamics, density frame,
metropolis algorithm

††preprint:

##  I Introduction

###  I.1 Physical motivation

Nuclear collisions at high energy exhibit remarkable collective flows, which
have been analyzed with considerable success using relativistic hydrodynamics
[1]. For large nuclei, ideal hydrodynamics provides a reasonable description
of the observed flows. Viscous corrections are then incorporated by simulating
(a version of) the relativistic Navier-Stokes equations, improving the
description of the data and clarifying the theoretical consistency of the
simulations. These simulations fit the shear viscosity to entropy ratio of QCD
around the crossover temperature. Current Bayesian fits give $\eta/s\simeq
2\,\hbar/4\pi k_{B}$ [2, 3, 4, 5], which indicates that the medium is
remarkably strongly coupled, with relaxation rates of the order $\sim
k_{B}T/\hbar$.

The hydrodynamic description of heavy ion collisions can be tested by
examining the collisions of light nuclei, such as ${\rm d+Au}$ and ${\rm
He}^{3}+{\rm Au}$ at the Relativistic Heavy Ion Collider (RHIC) and proton-
nucleus collision at the Large Hadron Collider (LHC). Remarkably, these events
also exhibit correlations which are indicative of the collective flow
[*[SeeforexampleSec.3.2.2andSect3.2.3in:][.]Arslandok:2023utm]. However, it
must be emphasized, that the hydrodynamic description of these events is
breaking down. This is in part because (a suitably defined) mean free path
$\ell_{\rm mfp}$ has become comparable to the system size, and in part because
the total number of particles produced in these events $N_{\rm ch}$ is
becoming small, which leads to large fluctuations. One of the motivations for
the current paper is to analyze thermal fluctuations in relativistic
dissipative systems, with the ultimate goal of describing small colliding
systems.

The current manuscript is also motivated by two classical $2^{\rm nd}$ order
phase transitions in QCD, which may be an observable in heavy ion collisions.
In both of these transitions incorporating thermal fluctuations into the
hydrodynamic description is essential to describing the underlying physics in
the critical region. The first phase transition is the $O(4)$ chiral
transition of QCD, which is a $2^{\rm nd}$ order phase transition in the limit
of two massless quark flavors [7, 8]. There is a strong motivation from
lattice QCD to look for signatures of the $O(4)$ transition in the heavy ion
data at the LHC [9, 10]. At lower temperatures and higher baryon density,
strong theoretical arguments suggest that there should be a Ising critical
point in the $(T,\mu_{B})$ plane [11], with $T$ and $\mu_{B}$ being the
temperature and the baryon chemical potential, respectively. Currently at
RHIC, there is an ongoing search for the Ising critical point where the beam
energy is scanned in an effort to tune the baryon chemical potential to the
critical region of the phase diagram [12].

###  I.2 The Metropolis algorithm for relativistic hydrodynamic fluctuations

We are also motivated to study thermal fluctuations in relativistic fluids by
algorithmic developments and the mathematical structure of stationary
stochastic processes. In statistical mechanics, the Metropolis algorithm is
used to generate field configurations of a field $\phi$ from a known
probability distribution, $P[\phi]\propto e^{S[\phi]}$, where $S[\phi]$ is the
entropy [13]. In the algorithm a proposal is made for a change in the fields,
$\phi\rightarrow\phi+\Delta\phi$. This proposal is accepted or rejected
according to the magnitude and sign of the change in the entropy, $\Delta
S\equiv S[\phi+\delta\phi]-S[\phi]$. The Metropolis steps are guaranteed to
converge to the required equilibrium distribution. Recently, in the context of
simulating the $O(4)$ critical point, we simulated the stochastic diffusion of
a conserved charge coupled to the order parameter using a variant of the
Metropolis algorithm [14]. A proposal is made for the transfer of charge
between the fluid cells, and this proposal is then accepted or rejected based
on the change in the entropy of the system. For small enough time steps the
Metropolis updates naturally reproduce the Langevin dynamics of the diffusion
equation. The equivalence of the Metropolis and Langevin dynamics for small
time steps $\Delta t$ has been repeatedly observed over the years [15, 16].

The advantages of a Metropolis based approach is that detailed balance and the
Fluctuation-Dissipation-Theorem (FDT) are automatically preserved,
independently of $\Delta t$, which guarantees that Markov-chain will
equilibrate to a specific action. In non-linear theories this simplifies the
renormalization of the theory and clarifies discretization ambiguities that
arise in non-linear Langevin equations [17, 18]. For fluids at rest, the
Metropolis algorithm has been used to implement the non-linear Langevin
dynamics in a number of challenging applications, leading to the study of
sphaleron transitions of hot QCD [16], the real time dynamics of $O(4)$
critical point in QCD [14, 19], and the dynamics of Model B in the Halperin-
Hohenberg classification of dynamical critical phenomena [20].

Our principal task in this and a companion paper is to generalize the
Metropolis approach to relativistic fluids in general coordinates, where the
(somewhat complicated) form of the dissipative strains should arise naturally
from the accept/reject steps of the Metropolis algorithm. As a first step, in
this paper we will consider the diffusion of charge in a relativistic fluid.

###  I.3 Causality, second order hydrodynamics, and the Metropolis algorithm

To understand the issues that arise with relativity, consider the relativistic
advection-diffusion equation in flat spacetime in the Landau-Lifshitz frame,
momentarily neglecting the stochastic noise for simplicity. We are considering
a fluid moving with three velocity $v^{i}$ in the lab frame and following the
diffusion of a dilute conserved charge within the fluid,
$\partial_{\mu}J^{\mu}=0$. The four velocity is $u^{\mu}=(\gamma,\gamma{\bf
v})$ and the local charge density in the rest frame of the fluid is
$n_{\rm\scriptscriptstyle LF}=-u_{\mu}J^{\mu}$. Landau and Lifshitz define a
hydrodynamic frame where the chemical potential is given by the value of $n$
and ideal equation of state to all orders in the gradient expansions, i.e.
$n_{\rm\scriptscriptstyle LF}=\chi\,\mu_{\rm\scriptscriptstyle LF}$ where
$\chi$ is the charge susceptibility [21].

The problem with the covariant Landau-Lifshitz approach is that the diffusive
current, which was spatial in the rest frame of the fluid
$j_{D}^{i}=-D\partial^{i}n$, involves time derivatives in an arbitrary frame.
This leads to equations which are second order in time, which in turn leads to
runaway solutions and other pathological behavior [22, 23]. One way to correct
this pathology is to promote the diffusive current to an additional dynamical
field which relaxes on collisional timescale to the expected form. This
procedure results in Maxwell-Catteneo or Israel Stewart type equations.

There is merit to the Israel-Stewart approach – it provides an effective way
to realize the dynamics of the relativistic diffusion equation and the
relativistic Navier Stokes equations more generally. Indeed, almost all large
scale simulations of flow in heavy ion collisions are based on this approach.
However, the Israel-Stewart formulation involves fast variables whose physical
significance should be questioned. Indeed, there have been many reformulations
of viscous hydrodynamics in the relativistic domain [24, 25, 26, 27, 28, 29],
and each of these reformulations involve some additional variables (or non-
hydrodynamic modes) which relax quickly to the form constrained by “first-
order” hydrodynamics [30]. A kind of theorem has emerged, which states that it
is impossible to construct a causal and stable relativistic theory of
hydrodynamics without incorporating non-hydrodynamic modes [31, 32, 33].

In this paper we will investigate an alternative to the Israel-Stewart
approach developed to describe fluids without an underlying boost symmetry
[34, 35, 36]. In particular, we found the “density frame” discussed by Armas
and Jain a clarifying formalism when implementing Metropolis updates [36]. The
density frame has no non-hydrodynamic modes and no additional parameters
compared to the Landau theory of first order hydrodynamics, at the price of
not being fully boost invariant. In hydrodynamics without boosts the
constitutive relations are written down for setups where the underlying
interactions are not Lorentz invariant and thus there is a preferred “lab”
frame111For a physical example, consider a fluid flowing over a table and
study the diffusion of charge in this background fluid flow. The table sets a
preferred lab frame.. If the microscopic interactions happen to be Lorentz
invariant, the additional boost symmetry imposes relations between the
coefficients of the gradient expansion, which is an expansion in lab-frame
spatial derivatives. Indeed, the equations of motion in the density frame
follow from the Landau ones if the ideal equations are used to rewrite lab-
frame time derivatives appearing in the dissipative strains as spatial
derivatives. With this rewrite the equations are strictly first order in time
and are stable.

The procedure amounts to a non-covariant choice of hydrodynamic frame where
the relation between the chemical potential and the lab frame charge per
volume $J^{0}$ is given by ideal hydrodynamics at all orders in the gradient
expansion, i.e. $J^{0}=\chi\mu u^{0}$. The frame choice and the resulting
equations of motion in the density frame are not invariant under Lorentz
transformations; but, they are invariant under Lorentz transformations
followed by a change of hydrodynamic frame, which reparametrizes the
hydrodynamic fields. The results obtained in different Lorentz frames will
vary, but the variation is beyond the accuracy of the diffusion equation.
Different choices of Lorentz invariant hydrodynamic frames, such as the
Landau, Eckart or BDNK222The acronym is short for Bemfica, Disconzi, and
Noronha [28] and Kovtun [29]. These authors studied a class of Lorentz
invariant frames where the temperature is related to the Landau choice up to
derivatives. choices, will also give different (if Lorentz invariant) results
to this accuracy. The density frame approach was known to “work” in simple
cases, but the work on hydrodynamics without boosts formalized the procedure.
Finally, Armas and Jain made an important connection of this approach to
modern treatments of hydrodynamics [37, 38] where the conserved charge and the
corresponding canonical conjugates (for instance a $U(1)$ charge $Q$ and the
associated phase $\varphi$) play a dual role [36]. The numerical utility of
the symplectic structure inherent in this duality remains to be fully
exploited.

The structure of the current paper is as follows. Section II discusses the
relativistic advection-diffusion equation in the Landau and density frames,
and derives the density frame constitutive relation from covariant kinetic
theory. Section III compares the density frame relativistic advection-
diffusion equation with relativistic kinetics. As discussed above, the density
frame is not Lorentz invariant, but it is invariant under Lorentz
transformations followed by a reparametrization of the hydrodynamic variables.
In a specific test problem discussed in Section III, we show that, in the
regime of validity of hydrodynamics, the deviations of the density frame from
a underlying covariant microscopic theory are controllably small, even for
highly boosted fluids. We also study the convergence of the gradient expansion
of the density frame, making connections with Lorentz covariant approaches.
Stochastic dynamics in the density frame is studied in Section IV. Since the
hydrodynamics equations are defined using a given foliation of space-time, the
Metropolis updates discussed above are very natural and easy to implement. One
simply makes random transfers of charge in between ideal advective steps.
These charge transfers are then accepted or rejected using the entropy defined
on the spatial slice as the statistical weight. We present a first study of
numerical correlation functions from the Metropolis algorithm for an
equilibrated boosted fluid in Section IV. Although we have used the Metropolis
algorithm for the relativistic hydrodynamics in the density frame, it should
be useful for other approaches to stochastic relativistic hydrodynamics, e.g.
approaches based on BDNK [39, 40] or Israel-Stewart [41]. Finally, in Section
V we conclude with a short discussion of the next steps.

##  II The advection diffusion equation

###  II.1 Setup and first order hydrodynamics in the Landau frame

We are considering the advection and diffusion of a charge in fluid moving at
relativistic speeds in flat spacetime, $\eta_{\mu\nu}=(-,+,+,+)$. The charge
density is low and the temperature and flow velocity
$u^{\mu}=(\gamma,\gamma{\bf v})$ may be considered fixed. In first order
hydrodynamics in the Landau frame the conserved current obeys [21]

$\partial_{\mu}J^{\mu}=0\,,\qquad J^{\mu}\equiv n_{{\rm\scriptscriptstyle
LF}}u^{\mu}+j_{{\rm\scriptscriptstyle D,LF}}^{\mu}\,,$ (1)

where the first term in $J^{\mu}$ is ideal advection and the second term is
the diffusive correction, expressed as

$j_{{\rm\scriptscriptstyle
D,LF}}^{\mu}\equiv-T\sigma\Delta^{\mu\nu}\partial_{\nu}\hat{\mu}_{{\rm\scriptscriptstyle
LF}}\,.$ (2)

Here $T$ is the temperature, $\sigma$ is the conductivity,
$\hat{\mu}_{{\rm\scriptscriptstyle LF}}\equiv\mu/T$ is the scaled chemical
potential, thermodynamically conjugate to charge density
$n_{\rm\scriptscriptstyle LF}$. $\Delta^{\mu\nu}$ is the spatial projector

$\Delta^{\mu\nu}=\eta^{\mu\nu}+u^{\mu}u^{\nu}\,,$ (3)

and satisfies
$\Delta^{\mu\rho}\Delta_{\rho\nu}=\Delta^{\mu}_{\phantom{\nu}\nu}$. Since the
density is low, $n_{{\rm\scriptscriptstyle
LF}}=\chi\mu_{{\rm\scriptscriptstyle LF}}$ where $\chi$ is the (temperature
dependent) susceptibility.

In the next subsections we will briefly review the density frame pointing out
the differences with the Landau frame. To keep the presentation self contained
and pedagogical, sections II.2 and II.3 review a small portion of [36] with a
focus on the diffusion equation. Section II.4 describes how the density frame
constitutive relation arises in relativistic kinetic theory.

###  II.2 Thermodynamics of a boosted fluid

Consider a portion of a fluid in perfect global equilibrium within a lab frame
measurement volume, $V_{0}\equiv\int{\rm d}\Sigma_{0}$. The entropy, energy-
momentum, and charge on this slice are

$\mathcal{S}=V_{0}S\,,\qquad\mathcal{P}_{\mu}\equiv
V_{0}T^{0}_{\phantom{\nu}\mu}\,,\qquad\mathcal{N}=V_{0}J^{0}\,.$ (4)

We will notate the charge density with $N\equiv J^{0}$, and the temporal
components of the energy momentum tensor with
$(E,M^{i})\equiv(T^{00},T^{0i})$, reserving the calligraphic symbols
$\mathcal{S}$, $\mathcal{P}_{\mu}$ and $\mathcal{N}$ for the charges, as
opposed to the charge densities, $S,E,M,N$. When a boost symmetry is
ultimately adopted, the entropy will take the form $S=su^{0}$ where $s$ is a
Lorentz scalar.

Given the conserved charges $\mathcal{P}_{\mu}$ and $\mathcal{N}$, the
temperature, velocity and chemical potential can be determined from the micro-
canonical equation of state $S(\mathcal{P},\mathcal{N},V_{0})$

$\displaystyle{\rm d}\mathcal{S}=$ $\displaystyle-\beta^{\mu}{\rm
d}\mathcal{P}_{\mu}-\hat{\mu}\,{\rm d}\mathcal{N}+p\beta^{0}{\rm d}V_{0}\,,$
(5) $\displaystyle=$ $\displaystyle-\beta^{0}{\rm
d}\mathcal{P}_{0}-\beta^{0}v^{i}{\rm d}\mathcal{P}_{i}-\hat{\mu}\,{\rm
d}\mathcal{N}+p\,\beta^{0}{\rm d}V_{0}\,,$ (6)

where $\hat{\mu}\equiv\mu/T$ and $v^{i}\equiv\beta^{i}/\beta^{0}$ is the fluid
velocity. The Gibbs-Duhem relation follows from the extensivity of the system

$S=-\beta^{\mu}T^{0}_{\phantom{\nu}\mu}-\hat{\mu}J^{0}+p\,\beta^{0}\,.$ (7)

Then for small charge densities the entropy density as a function of $N$ takes
the form

$S(N)=\underbrace{S_{1}(E,M)}_{\mbox{const}}-\frac{\beta^{0}}{2\chi^{00}}N^{2}\,.$
(8)

So far we have not used boost symmetry, and the parameters, such as
$\beta^{0}$, $\mu_{0}$ and $\chi^{00}$ are functions of $E$ and $M$. After
imposing boost symmetry in the next paragraph, $\beta^{0}$ will be a Lorentz
scalar $\beta$ times $u^{0}$, and $\chi^{00}$ will be a Lorentz scalar
$\chi(\beta)$ times $u^{0}u^{0}$

$\beta^{0}=\beta u^{0}\,,\qquad\chi^{00}=\chi(\beta)u^{0}u^{0}\,,$ (9)

justifying the notation a posteriori.

When the fluid has an underlying Lorentz symmetry the dependence on the
velocity is determined by the symmetry. Before imposing the symmetry, the
conservation laws take the form

$\displaystyle\partial_{t}N+\partial_{i}J^{i}=$ $\displaystyle 0\,,$
$\displaystyle\partial_{t}E+\partial_{i}T^{i0}=$ $\displaystyle 0\,,$
$\displaystyle\partial_{t}M^{j}+\partial_{i}T^{ij}=$ $\displaystyle 0\,.$ (10)

At zeroth order in the gradient expansion (ideal hydrodynamics) the three
current $J^{i}$, the energy flux $T^{i0}$, and the spatial stress tensor
$T^{ij}$ are algebraically related to the conserved charges $E$, $M^{i}$, and
$N$. Using the conservation laws, the thermodynamic relations (eqs. (5) and
(7)), the symmetry of the stress tensor imposed by Lorentz invariance
$T^{i0}=M^{i}$, and requiring that $\partial_{t}S+\partial_{i}S^{i}$ be non-
negative, fixes the form of fluxes $J^{i}$, $M^{i}$, $T^{ij}$ and entropy
current $S^{i}$ to the form of ideal hydrodynamics [36]. In particular, the
charges and entropy are parametrized by $\beta^{\mu}\equiv\beta u^{\mu}$ and
the chemical potential $\hat{\mu}$

$E=(e+p)(u^{0})^{2}-p,\qquad M^{i}=(e+p)u^{0}u^{i},\qquad N=nu^{0},\qquad
S=su^{0}.$ (11)

Here $e$, $p$, $n$ and $s$ are scalar functions of
$\beta\equiv\sqrt{-\beta^{\mu}\beta_{\mu}}$ and $\hat{\mu}$ as given by the
equilibrium equation of state. At small chemical potentials, the local charge
density and entropy takes the form

$n=\chi\,\mu\,,\qquad
s(\beta,\mu)=s_{1}(\beta)-\tfrac{1}{2}\beta\chi\mu^{2}\,.$ (12)

where $\chi(\beta)$ is a function of temperature. Comparing the definitions in
eq. (8) and eq. (12), we find $\chi^{00}/\beta^{0}=T\chi u^{0}$ as claimed
above.

In the density frame the familiar relations of ideal hydrodynamics, eqs. (11),
serve to define the temperature, chemical potential, and flow velocity in
terms of the lab frame charges, $E$, $M^{i}$ and $N$, at every order in the
gradient expansion, i.e. the relation between the charges $E,M^{i},N$ and
their conjugates do not receive viscous corrections. In particular, the
chemical potential in the density frame is defined at all orders in the
gradient expansion as

$\mu=\frac{J^{0}}{\chi u^{0}}\,.$ (13)

This definition should be contrasted with the chemical potential in the
Lorentz invariant Landau frame where

$\mu_{{\rm\scriptscriptstyle LF}}=-\frac{u_{\mu}J^{\mu}}{\chi}\,.$ (14)

With the density frame definition a fiducial observer needs to count the
charge in a given measurement volume in order to determine the chemical
potential. With the Landau frame definition, the three current $J^{i}$ also
needs to be measured. Thus, the Landau frame involves counting the charges at
different times in order to define the chemical potential.

###  II.3 The advection diffusion equation in the density frame

Here we will derive the density frame equations of motion by first considering
fluids without a boost symmetry, and then specializing the equations to
Lorentz covariant fluids. An example of a two dimensional non-Lorentz
invariant fluid is a fluid flowing over a flat surface at relativistic speeds.
The diffusion of a charge in this fluid depends on the speed of the fluid
relative to the surface.

The advection-diffusion equation in the density frame consists of the
conservation law

$\partial_{t}N+\partial_{i}J^{i}=0\,,$ (15)

together with a constitutive relation for the diffusive current $J_{D}^{i}$

$J^{i}\equiv Nv^{i}+J_{D}^{i}\,.$ (16)

The diffusive current is expanded in spatial gradients of the conserved
charge, or its thermodynamic conjugate $\hat{\mu}$. The most general form of
$J_{D}^{i}$ at first order in gradients of $\hat{\mu}$ is

$J_{D}^{i}=-\frac{\sigma_{\parallel}(\beta^{0},v)}{\beta^{0}}\,\hat{v}^{i}\hat{v}^{j}\partial_{j}\hat{\mu}-\frac{\sigma_{\perp}(\beta^{0},v)}{\beta^{0}}\left(\delta^{ij}-\hat{v}^{i}\hat{v}^{j}\right)\partial_{j}\hat{\mu}\,,$
(17)

where $\hat{\mu}\equiv\partial S/\partial N=\beta^{0}N/\chi^{00}$ and
$\hat{v}^{i}=v^{i}/|v|$ is a flow unit vector. The first and second terms on
the right hand side of eq. (17) capture the diffusion parallel and
perpendicular to the fluid motion respectively. Demanding that entropy
production be positive leads to the requirement that
$\sigma_{\parallel}(\beta^{0},v)>0$ and $\sigma_{\perp}(\beta^{0},v)>0$, but
no further constraints can be derived on general grounds.

Lorentz invariant fluids can be treated as a special case of eq. (17). Indeed,
the boost symmetry determines the dependence of $\sigma_{\parallel}$ and
$\sigma_{\perp}$ on the velocity. The easiest way to derive this relation is
to return momentarily to the Landau frame. Comparison to the density frame
form gives

$N=n_{{\rm\scriptscriptstyle LF}}u^{0}+j_{\rm\scriptscriptstyle D,LF}^{0}\,,$
(18)

and

$J_{D}^{i}=J^{i}-Nv^{i}=(\Delta^{i}_{\phantom{\nu}\alpha}-v^{i}\Delta^{0}_{\phantom{\nu}\alpha})j_{\rm\scriptscriptstyle
D,LF}^{\alpha}\,.$ (19)

We now use the lowest order equation of motion,

$\partial_{t}\hat{\mu}\simeq-v^{j}\partial_{j}\hat{\mu}\,,$ (20)

to approximate the Landau frame expression for the diffusive current

$j_{{\rm\scriptscriptstyle D,LF}}^{\alpha}\simeq-T\sigma\left(\Delta^{\alpha
j}-\Delta^{\alpha 0}v^{j}\right)\partial_{j}\hat{\mu}\,.$ (21)

Substituting eq. (21) into eq. (19) gives

$J_{D}^{i}=-T\,{\sigma}^{ij}\,\partial_{j}\hat{\mu}\,,$ (22)

where we have defined a frequently occurring matrix

$\displaystyle T\sigma^{ij}=$ $\displaystyle
T\sigma\left(\Delta^{i}_{\phantom{\nu}\alpha}-v^{i}\Delta^{0}_{\phantom{\nu}\alpha}\right)\left(\Delta^{j}_{\phantom{\nu}\beta}-v^{j}\Delta^{0}_{\phantom{\nu}\beta}\right)\,\Delta^{\alpha\beta}\,,$
(23) $\displaystyle=$ $\displaystyle
T\sigma\left(\delta^{ij}-v^{i}v^{j}\right)\,.$ (24)

Comparison with the general form in eq. (17) shows that

$\frac{\sigma_{\parallel}(\beta^{0},v)}{\beta^{0}}=\frac{T\sigma(\beta)}{\gamma^{2}}\,,\qquad\qquad\frac{\sigma_{\perp}(\beta^{0},v)}{\beta^{0}}=T\sigma(\beta)\,.$
(25)

In summary, in the density frame the equation of motion is

$\partial_{t}N+\partial_{i}(Nv^{i})=\partial_{i}\left(T\sigma^{ij}\partial_{j}\hat{\mu}\right)\,,$
(26)

and when $\hat{\mu}$ is written in terms of the charge $N=\chi\mu u^{0}$, we
arrive at an advection-diffusion equation

$\partial_{t}N+\partial_{i}(Nv^{i})=\partial_{i}\left(D^{ij}\partial_{j}N\right)\,,$
(27)

with a diffusion matrix

$D^{ij}=\frac{D}{\gamma}\left(\delta^{ij}-v^{i}v^{j}\right)\,.$ (28)

Here $D=T\sigma/\chi$ is the scalar diffusion coefficient of the Landau frame.
Apart from the tensor structure the equation is numerically similar to the
non-relativistic advection-diffusion equation and can be solved by familiar
numerical methods333In our numerical tests without noise we used an IMEX
scheme with a standard advection step and a Crank-Nicholson like implicit step
using the PETSc library [48]. .

The $\gamma$ factors in the diffusion matrix can be easily understood
physically. The diffusion coefficient has units of distance squared per time.
The rate of transverse diffusion is suppressed relative to a fluid at rest by
one factor of $\gamma$ due to time dilation. The rate of longitudinal
diffusion is suppressed by three factors of $\gamma$ due to time dilation and
length contraction, i.e. each spatial step in the random walk is length
contracted by $\gamma$ and the steps add in square.

###  II.4 The density frame from relativistic kinetics

In this subsection we will briefly describe how the density frame constitutive
relation arises naturally in relativistic kinetic theory. Specifically, we
will show how the conductivity matrix $T\sigma^{ij}$ follows from covariant
kinetics and find how this matrix is determined by the particles through the
first viscous correction $\delta f$ to the phase-space distribution function.
For simplicity, we will assume a relaxation time approximation and consider
single species of classical relativistic particles, which carry the charge of
the system

$p^{\mu}\partial_{\mu}f=-{\mathcal{C}}_{p}\,\delta f\,.$ (29)

Here ${\mathcal{C}}_{p}$ is a momentum dependent parameter controlling the
collision rate in the rest frame of the medium.

In global equilibrium the phase space distribution function is characterized
by a constant chemical potential, temperature, and flow velocity. If the
density of the charged particles depends slowly on space and time then the
parameter $\hat{\mu}(t,{\bf x})$ is no longer a constant but reflects this
dependence

$f_{0}(t,{\bf x},{\bf p})=e^{\hat{\mu}(t,{\bf x})}e^{\beta^{\mu}p_{\mu}}\,.$
(30)

In the density frame $\hat{\mu}(t,{\bf x})$ is adjusted to reproduce the
charge density in the lab frame $J^{0}$, while in the Landau frame
$\hat{\mu}(t,{\bf x})$ is adjusted to reproduce the charge in the rest frame,
$n(t,{\bf x})=-u_{\mu}J^{\mu}$. These two definitions of the chemical
potential agree when gradients are neglected, and in this case $f_{\rm
eq}(t,{\bf x},{\bf p})$ is a solution to the Boltzmann equation.
$\hat{\mu}(t,{\bf x})$ obeys the equations of ideal advection equation at
lowest order

$u^{\mu}\partial_{\mu}\hat{\mu}\simeq 0\,.$ (31)

Following a standard procedure to find the first viscous correction [42], we
parameterize $f=f_{0}+\delta f(t,{\bf x},{\bf p})$ and solve for $\delta f$
order by order in the gradients

$f_{0}\,p^{\mu}\partial_{\mu}\hat{\mu}=-{\mathcal{C}}_{p}\,\delta f\,.$ (32)

In the Landau frame one decomposes the gradient into its temporal and spatial
components as

$\partial_{\mu}\hat{\mu}=-u_{\mu}u^{\alpha}\partial_{\alpha}\hat{\mu}+\Delta_{\mu}^{\phantom{\nu}\alpha}\,\partial_{\alpha}\hat{\mu}\,.$
(33)

We neglect the temporal term in Eq. (33) exploiting the lowest order equations
of motion, Eq. (31). Then we substitute into Eq. (32), which leads to the
familiar form of the first viscous correction in the Landau frame

$\delta f_{\rm\scriptscriptstyle
LF}=-{\mathcal{C}}_{p}^{-1}\,f_{0}\,p^{\alpha}\nabla_{\alpha}\hat{\mu}_{\rm\scriptscriptstyle
LF}\,,$ (34)

where $\nabla_{\alpha}=\Delta_{\alpha}^{\phantom{\nu}\mu}\partial_{\mu}$.
Evaluating the diffusive current

$j^{\mu}_{{\rm\scriptscriptstyle
D,LF}}=\int_{p}\frac{d^{3}p}{(2\pi)^{3}}\,\frac{p^{\mu}}{p^{0}}\delta
f_{\rm\scriptscriptstyle LF}\,,$ (35)

yields the expected form of the current in the Landau frame

$j_{\rm\scriptscriptstyle
D,LF}^{\mu}=-T\sigma\Delta^{\mu\nu}\partial_{\nu}\hat{\mu}_{\rm\scriptscriptstyle
LF}\,.$ (36)

The conductivity in this expression is defined from the transport integrals

$T\sigma\Delta^{\mu\nu}\equiv\Delta^{\mu}_{\phantom{\nu}\alpha}\Delta^{\nu}_{\phantom{\nu}\beta}I^{\alpha\beta}\,,\qquad
I^{\alpha\beta}\equiv\int\frac{d^{3}p}{(2\pi)^{3}p^{0}}\,{\mathcal{C}}_{p}^{-1}f_{0}\,p^{\alpha}p^{\beta}\,.$
(37)

In the density frame one proceeds similarly, but uses the lowest order
equations in the lab frame

$\partial_{t}\hat{\mu}=-v^{i}\partial_{i}\hat{\mu}\,,$ (38)

which yields the form of the first viscous correction

$\delta
f=-{\mathcal{C}}_{p}^{-1}f_{0}(p^{i}-p^{0}v^{i})\,\partial_{i}\hat{\mu}\,.$
(39)

The diffusive current $J_{D}^{i}$ in the density frame is the difference
between the current and the ideal advection $J^{0}v^{i}$

$J_{D}^{i}=\int\frac{d^{3}p}{(2\pi)^{3}p^{0}}\left(p^{i}-p^{0}v^{i}\right)f\,.$
(40)

Substituting the approximate distribution function $f_{0}+\delta f$ leads to
an appealing positive definite symmetric matrix which determines the mean
current444As discussed below, the matrix $\mathcal{K}^{ij}$ also determines
functional form of the noise in the density frame. In Ref. [36] the matrix
$\mathcal{K}$ is written as $D_{jj}^{\rho\sigma}$. This section shows how this
form arises in a microscopic theory.

$J_{D}^{i}=\mathcal{K}^{ij}(-\partial_{j}\hat{\mu})\,,\qquad\mathcal{K}^{ij}\equiv\int\frac{d^{3}p}{(2\pi)^{3}p^{0}}\,{\mathcal{C}}_{p}^{-1}f_{0}\,(p^{i}-p^{0}v^{i})(p^{j}-p^{0}v^{j})\,.$
(41)

Noting that

$p^{i}-p^{0}v^{i}=\left(\Delta^{i}_{\phantom{\nu}\alpha}-v^{i}\Delta^{0}_{\phantom{\nu}\alpha}\right)\Delta^{\alpha}_{\phantom{\nu}\beta}\,p^{\beta}\,,$
(42)

we find that $\mathcal{K}^{ij}$ has the expected density frame form

$\displaystyle\mathcal{K}^{ij}=$
$\displaystyle\left(\Delta^{i}_{\phantom{\nu}\alpha}-v^{i}\Delta^{0}_{\phantom{\nu}\alpha}\right)\left(\Delta^{j}_{\phantom{\nu}\beta}-v^{j}\Delta^{0}_{\phantom{\nu}\beta}\right)\,T\sigma\Delta^{\alpha\beta}\,,$
(43a) $\displaystyle=$ $\displaystyle
T\sigma\left(\delta^{ij}-v^{i}v^{j}\right)\,,$ (43b)

where we used the integrals defined in eq. (37). To summarize this section, we
have shown how the form of the dissipative current in the density frame arises
in covariant kinetic theory (eqs. (41) and (43)). This form is not covariant
although the underlying kinetic theory is covariant. This arises because we
are trying to approximate the full results of kinetic theory in a specific
frame.

##  III Comparison with a kinetic model

###  III.1 A random walk of massless particles: static case

In this and the next subsections we will study an analytically tractable
covariant kinetic model in $1+1$ dimensions and investigate how the current in
this model approaches the density frame constitutive relation. The model
consists of massless “particles” moving with the speed of light along a one-
dimensional line. The particles experience Poissonian random kicks which
changes the direction of their velocities. The dynamical evolution of this
system can be mapped to the famous telegraph equation from which an all-orders
gradient expansion can be derived. We will then compare the results of this
all-order gradient expansion with the predictions from the density frame
formalism.

In this subsection we will analyze this model in a static background case and
then in the next subsection consider a fluid background moving with constant
velocity. The results in this subsection are not new and can be found in many
places, but they will serve to define the terms for the boosted case.

Let us denote the respective densities of the left and right moving particles
as $n_{-}$ and $n_{+}$ (see Fig. 1).

Figure 1: Particles moving left and right with density $n_{-}$ and $n_{+}$,
respectively.

Due to the random kicks, the particles change direction with the rate
$\Gamma\equiv 1/2\tau_{R}$ where $\tau_{R}$ is the relaxation time. The
density of the left/right movers obey the following kinetic equation

$\displaystyle\partial_{t}n_{+}+c\,\partial_{x}n_{+}$ $\displaystyle=$
$\displaystyle-\Gamma\left(n_{+}-n_{-}\right)\,,$ (44)
$\displaystyle\partial_{t}n_{-}-c\,\partial_{x}n_{-}$ $\displaystyle=$
$\displaystyle-\Gamma\left(n_{-}-n_{+}\right)\,.$ (45)

Adding eqs. (44) and (45) and rearranging leads to the conservation equation

$\displaystyle\partial_{t}n+\partial_{x}j=0\,,$ (46)

where $n=n_{+}+n_{-}$ is the density and $j=c(n_{+}-n_{-})$ is the current.
Similarly, subtracting eqs. (44) and (45) and rearranging gives the relaxation
equation for the current

$\partial_{t}j+c^{2}\partial_{x}n=-\frac{j}{\tau_{R}}\,.$ (47)

As usual, this set of first order equations can be expressed as a second order
equation for $n$

$\displaystyle\partial_{t}^{2}n+\frac{1}{\tau_{R}}\partial_{t}n-c^{2}\partial^{2}_{x}n=0\,,$
(48)

which is known as the Telegraph equation. The exact Green functions associated
with this system are known analytically and can be expressed in terms of
modified Bessel functions which are given in the Appendix A for convenience.
For simplicity, we will set $c=1$ for the remainder of this section.

Consider some initial state at $t=0$ (not necessarily in equilibrium) that is
specified by two independent functions $n(t=0,x)=n_{0}(x)$ and
$j(t=0,x)=j_{0}(x)$. Let us further assume that the initial current can be
expressed as a gradient expansion of the initial density, with some
coefficients
$j_{0}=-\sum_{\ell=0}^{\infty}\tau_{R}^{2\ell+1}b_{2\ell+1}\partial_{x}^{2\ell+1}n_{0}(x)$.
Note that the coefficients $b_{2\ell+1}$ partially characterize the initial
state and they are dimensionless by definition.

By using the exact Green functions, it can be shown that at late times
$t\gtrsim\tau_{R}$ the system obeys a universal dispersion relation

$\displaystyle
j=-\sum_{\ell=0}^{\infty}c_{\ell}\tau_{R}^{2\ell+1}\partial_{x}^{(2\ell+1)}n\,,$
(49)

where

$c_{\ell}=(-1)^{\ell}C(\ell)+\sum_{m=1}^{\ell+1}e^{-m(t/\tau_{R})}P_{\ell,m}(t)\,,$
(50)

with $C(\ell)=(2\ell)!/(\ell!(\ell+1)!)$ being the $\ell^{\rm th}$ Catalan
number and $P_{\ell,m}(t)$ is a polynomial in $t$ of order $\ell-m+1$ whose
coefficients depend on $b_{i\leq\ell}$. It is clear from this expansion that
the dependence on the initial conditions decays exponentially for
$t\gtrsim\tau_{R}$, meaning that the system thermalizes and obeys a universal
constitutive relation. The constitutive relation can be expressed as an all-
order gradient expansion whose coefficients are given by $(-1)^{\ell}C(\ell)$.
The Catalan numbers are nothing but the Taylor coefficients of the dispersion
curve $\omega(k)$ of the diffusion mode:

$\det\left(\begin{matrix}-i\omega&ik\\\
ik&-i\omega+1/\tau_{R}\end{matrix}\right)=0\Rightarrow\omega(k)=-\frac{i}{2\tau_{R}}\left(1-\sqrt{1-4k^{2}\tau_{R}^{2}}\right),$
(51)

associated with the linear system expanded around $k=0$. The dispersion
relation has a branch singularity at $k_{*}=\pm 1/(2\tau_{R})$. Furthermore
the _large $k$_ expansion is consistent with the stability and causality
conditions given in Refs. [31, 32, 33, 43] where the Telegraph equation was
previously analyzed in these terms.

The additional eigenfrequency from (51) is gapped and approaches
$\omega(k)=-i/\tau_{R}+\mathcal{O}((\tau_{R}k)^{2})$ for $k\rightarrow 0$. The
gapped modes are responsible for the exponential decay in eq. (50). As we will
see in the next section, when the fluid is moving with velocity $v$, the
density frame constitutive relation is approached on a short timescale of
order $\tau_{R}/\gamma$, which again reflects the dynamics of gapped non-
hydrodynamic modes.

###  III.2 A random walk of massless particles: moving fluid

Let us now consider the same kinetic model of the previous section in the
background of a fluid moving uniformly with velocity $v$ in the lab frame, as
shown in fig. 2.

Figure 2: Particles moving left and right with density $N_{-}$ and $N_{+}$,
respectively, in a background fluid moving with velocity $v$.

In the lab frame, the density of left and right movers obey

$\displaystyle\partial_{t}N_{+}+c\,\partial_{x}N_{+}$ $\displaystyle=$
$\displaystyle-\Gamma_{+}N_{+}+\Gamma_{-}N_{-}\,,$ (52)
$\displaystyle\partial_{t}N_{-}-c\,\partial_{x}N_{-}$ $\displaystyle=$
$\displaystyle-\Gamma_{-}N_{-}+\Gamma_{+}N_{+}\,,$ (53)

where transition rates are

$\displaystyle\Gamma_{+}=\frac{1}{2\tau_{R}}\sqrt{\frac{1-v/c}{1+v/c}}\,,\qquad\Gamma_{-}=\frac{1}{2\tau_{R}}\sqrt{\frac{1+v/c}{1-v/c}}\,.$
(54)

The appearance of the kinematical factors $\kappa_{\pm}=\sqrt{(1\pm v/c)/(1\mp
v/c)}$ can be understood in the following way: the right movers in the local
rest frame, denoted by the subscript $R$, follow the trajectory
$x_{R}=ct_{R}$. A Lorentz boost to the lab frame coordinates, $(t,x)$, leads
to the time dilation factor $t=\gamma(1+v/c)t_{R}=\kappa_{+}t_{R}$ where
$\gamma=(1-v^{2}/c^{2})^{-1/2}$ is the Lorentz factor. Since the mean free
path time in the local rest frame is $\Gamma^{-1}=2\tau_{R}$, the mean free
path time in the lab frame is therefore $\Gamma^{-1}_{+}=\kappa_{+}2\tau_{R}$.

Similarly to the static case, by adding and rearranging equations (52) and
(53), we get

$\displaystyle\partial_{t}N+\partial_{x}J$ $\displaystyle=$ $\displaystyle
0\,,$ (55) $\displaystyle\partial_{t}J+c^{2}\partial_{x}N$ $\displaystyle=$
$\displaystyle\frac{\gamma}{\tau_{R}}\left(v\,N-J\right)\,,$ (56)

where $N=N_{+}+N_{-}$ and $J=c\,(N_{+}-N_{-})$ denote the density and the
total current respectively. Following the density frame approach, we write the
current as a sum of the advective and diffusive parts,

$J=vN+J_{D}\,.$ (57)

The diffusive part satisfies the relaxation equation

$\partial_{t}J_{D}-v\partial_{x}J_{D}+\frac{c^{2}}{\gamma^{2}}\partial_{x}N+\frac{\gamma}{\tau_{R}}J_{D}=0\,.$
(58)

For clarity we will set $c=1$ for the rest of this section.

Let us now consider an initial state given at $t=0$ in the lab frame. Based on
our findings in the static case, we expect the system to lose information
about the initial conditions for $t\gtrsim\tau_{R}/\gamma$ and obey a
universal gradient expansion

$J_{D}=-\frac{1}{\gamma^{2}}\sum_{n=1}^{\infty}c_{n}\left(\frac{\tau_{R}}{\gamma}\right)^{n}\partial_{x}^{(n)}N\,.$
(59)

Note that due to the nonzero velocity that explicitly breaks parity, we expect
both even and odd terms in the derivative expansion as opposed to the static
case, which only contains odd terms. We factored out an overall factor of
$\gamma^{-2}$ and the characteristic time scale $\tau_{R}/\gamma$ explicitly
to simplify the expressions for $c_{n}$. As in the static case, the
coefficients of the gradient expansion, $c_{n}$, can be calculated from the
dispersion relation

$\displaystyle\det\left(\begin{matrix}-i\omega+vik&ik\\\
i\gamma^{-2}k&-i\omega-vik+\gamma/\tau_{R}\end{matrix}\right)=0\,,$ (60)
$\displaystyle\omega(k)=\frac{-i}{2(\tau_{R}/\gamma)}\left(1-\sqrt{1-4k^{2}(\tau_{R}/\gamma)^{2}-4vik(\tau_{R}/\gamma)}\right)\,.$
(61)

The branch singularity in this case is in the complex plane, $k_{*}=(\pm
1-iv\gamma)/(2\tau_{R})$. Note that the radius of convergence,
$|k_{*}|=\gamma/(2\tau_{R})$, grows with $\gamma$ so that the hydrodynamic
description applies to modes of wavenumbers $k\lesssim\gamma/\ell_{\rm mfp}$
in the lab frame. Here and below we define the rest-frame mean-free-path:

$\ell_{\rm mfp}\equiv c/\Gamma\equiv 2\tau_{R}\,.$ (62)

The other branch in the dispersion relation is gapped,
$\omega(k)=-i/(\tau_{R}/\gamma)$ for $k\rightarrow 0$. This branch describes a
non-hydrodynamic mode decaying exponentially on a timescale of
$\tau_{R}/\gamma$. Notably, the decay time is not time-dilated, but instead
time-contracted by a factor of $\gamma$ relative to the static case.

The $c_{n}$ follow from Taylor expanding $\omega(k)$ in eq. (61) around $k=0$.
The even and odd terms are found as

$\displaystyle c_{2n+1}$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma^{2n}}\sum_{j=0}^{n}\frac{(-1)^{n+j}(2v\gamma)^{2j}(2n)!}{(n-j+1)(2j)!\left((n-j)!\right)^{2}}\,,$
(63) $\displaystyle=$
$\displaystyle\frac{1}{\gamma^{2n}}\frac{(-4)^{n}\Gamma(2n+3/2)}{\Gamma(n+3/2)\Gamma(n+2)}\,_{2}F_{1}\left(-1-n,-n,-2n-1/2,\gamma^{2}\right)\,,$
$\displaystyle c_{2n+2}$ $\displaystyle=$ $\displaystyle
2v\frac{1}{\gamma^{2n}}\sum_{j=0}^{n}\frac{(-1)^{n+j}(2v\gamma)^{2j}(2n+1)!}{(n-j+1)(2j+1)!\left((n-j)!\right)^{2}}\,,$
(64) $\displaystyle=$ $\displaystyle
2v\frac{1}{\gamma^{2n}}\frac{(-4)^{n}\Gamma(2n+5/2)}{\Gamma(n+5/2)\Gamma(n+2)}\,_{2}F_{1}\left(-1-n,-n,-2n-3/2,\gamma^{2}\right)\,.$

For reference we write down the first few terms below:

$\displaystyle
c_{1}=1,c_{2}=2v,c_{3}=-\left(1-5v^{2}\right),c_{4}=-2v\left(3-7v^{2}\right),c_{5}=2\left(1-14v^{2}+21v^{4}\right)\dots$
(65)

Just like the static case, the large $k$ expansion of the dispersion relation
given in eq. (61) necessarily satisfy the causality and stability conditions
of Refs. [31, 32, 33, 43].

Putting everything together, we find the gradient expansion in the lab frame
as

$\displaystyle J_{D}$ $\displaystyle=$
$\displaystyle-\frac{\tau_{R}}{\gamma^{3}}\partial_{x}N+(1-5v^{2})\frac{\tau_{R}^{3}}{\gamma^{5}}\partial_{x}^{3}N-2(1-14v^{2}+21v^{4})\frac{\tau_{R}^{5}}{\gamma^{7}}\partial_{x}^{5}N+\dots$
(66)
$\displaystyle+v\left[-2\frac{\tau_{R}^{2}}{\gamma^{4}}\partial_{x}^{2}N+2(3-7v^{2})\frac{\tau_{R}^{4}}{\gamma^{6}}\partial_{x}^{4}N+\dots\right]\,.$

Let us discuss this result. The leading term in eq. (66) agrees with the
prediction from the advection-diffusion equation in the density frame. In some
sense, the higher order corrections encode the underlying microscopic dynamics
and they restore Lorentz causality as can be seen from the large $k$
expansion, for instance. Of course Lorentz causality is violated if the
gradient expansion is truncated at any finite order.

One might ask when the physics of the covariant kinetic model is captured by
the density frame diffusion framework. We first require that the system
reaches local thermal equilibrium and therefore can be described by the
universal constitutive relation. This happens in a timescale of order
$t\sim\tau_{R}/\gamma$, which is set by the inverse frequency of the gapped
modes. This is illustrated in Fig. 3 where at $t=0$ we initialize a

Figure 3: (a) The evolution of a Gaussian drop of charge in the lab frame in a
kinetic model of Section III for a moving fluid with Lorentz factor
$\gamma=10$. The Gaussian has a lab frame width of $\sigma=L/\gamma$ where $L$
is fifty times the rest frame mean-free-path, $L=50\,\ell_{\rm mfp}$. We are
studying a very short time interval, $\Delta t=5\,\tau_{R}/\gamma$, i.e. a
time interval of order the rest frame relaxation-time _divided_ by $\gamma$.
Over this interval the drop is advected, but the diffusion of the drop is
negligible. (b) The relaxation of the lab frame diffusive current in the
kinetic model over the same short time interval as (a) starting from $J_{D}=0$
at $t=0$. The current at different times is compared to the leading order
density frame prediction, $J_{D}=-(D/\gamma^{3})\partial_{x}N$, which is
essentially time independent.

Gaussian drop of charge, $N_{0}(x)=N_{\rm max}\exp(-x^{2}/2\sigma^{2})$, and
set the initial diffusive current to zero in the kinetic model, $J_{D}=0$. For
the test case shown we took a fluid with $\gamma=10$ and set $\sigma=L/\gamma$
with $L/\ell_{\rm mfp}=50$. The left figure shows the time evolution of the
charge density in the lab frame. In the short period of time considered in the
figure $\Delta t\sim\tau_{R}/\gamma$, the charge does not have enough time to
diffuse and it is simply advected by the background flow. The right figure
shows the evolution of the lab frame diffusive current over the same time
period. Clearly the diffusive current in the kinetic model relaxes on a time
of order $\tau_{R}/\gamma$ to a steady state, and at $t=5\tau_{R}/\gamma$ the
steady state agrees with the leading derivative term of density frame gradient
expansion given in eq. (66) to a very good precision.

The relaxation timescale $\tau_{R}/\gamma$ is easily understood. Indeed, the
mean collision times of right and left movers are of the order
$\gamma\tau_{R}$ and $\tau_{R}/\gamma$ respectively. In equilibrium, where
$J=Nv$, the number of right and left movers is

$N_{+}^{\rm eq}=N\frac{1+v}{2}\,,\qquad N_{-}^{\rm eq}=N\frac{1-v}{2}\,,$ (67)

and thus for $v\rightarrow 1$ there are almost $N$ right movers and of order
$\sim N/\gamma^{2}$ left movers in the equilibrium sample. If all the
particles are initially left movers, then in a time of order
$\tau_{R}/\gamma$, these initial particles will scatter with probability one,
reaching approximate equilibrium with $N_{+}\simeq N$. Similarly, if all the
particles are initially right movers, then in a time of order
$\tau_{R}/\gamma$, $N/\gamma^{2}$ of these initial particles will scatter,
generating a yield of left movers commensurate with equilibrium, $N_{-}\sim
N/\gamma^{2}$. Thus, the typical equilibration time is of order
$\tau_{R}/\gamma$.

Of course, the applicability of the gradient expansion depends on the size of
the system as well, and we expect it to break down for smaller systems. In
Fig. 4 we show the diffusive current for the same setup as above with

Figure 4: Comparison of the first three orders of density frame gradient
expansion (DF: 1, 2, 3 given in eq. (66)) with the kinetic model, for the same
setup as Fig. 3 but for different rest frame system sizes $L$. For
$L=1\,\ell_{\rm mfp}$ the gradient expansion does not converge at all,
however, the gradient expansion becomes increasingly accurate with the
increase in system size.

$\gamma=10$ at $t=5\tau_{R}/\gamma$, but with different initial Gaussian sizes
ranging from $L=(1\ldots 8)\,\ell_{\rm mfp}$ in the rest frame. From the
figure we see that for a small system when $L$ is one mean-free-path, the
gradient expansion does not converge at all. As the system size gets larger,
the gradient expansion becomes more and more accurate. Note that the leading
term in the gradient remains qualitatively well behaved even for very small
systems. Remarkably, the convergence starts already when the system size is
only $3$ or $4$ mean-free-paths long, where the second and third order terms
in the gradient expansion capture the physics quite accurately. Notably in
this intermediate region, the effect of the second order term is clearly seen
by comparing the red (exact) and blue (first order) curves in the lower left
figure. This skewness is due to the background motion, and is absent in the
static case where parity is unbroken.

The convergence of the gradient expansion can be analyzed more rigorously. Let
us consider the time interval $\tau_{R}/\gamma\ll t\ll\tau_{R}\gamma^{3}$
where the system is in local thermal equilibrium but the initial wave packet
has not diffused yet. In this regime $N(t,x)\approx N_{0}(x)$, therefore, the
gradient expansion, eq. (59), becomes

$J_{D}\sim-\gamma^{-2}\sum_{n=1}^{\infty}c_{n}\left(\frac{\tau_{R}}{\gamma}\right)^{n}\partial_{x}^{(n)}N_{0}\sim-\gamma^{-2}\sum_{n=1}^{\infty}c_{n}\left(\frac{\tau_{R}}{\sqrt{2}\gamma\sigma}\right)^{n}e^{-\frac{x^{2}}{2\sigma^{2}}}H_{n}(x)\,,$
(68)

where $H_{n}(x)$ denotes the $n^{th}$ Hermite polynomial that follows from the
gradient expansion of the Gaussian wave-packet with the width $\sigma$. We do
not keep track of the overall magnitude of the initial density as it is
inconsequential for our argument. It is straightforward to observe that the
coefficients, $c_{n}$, grow exponentially in $n$ 555This can be seen from the
dispersion relation, eq. (61), whose Taylor coefficients around $k=0$ are
$c_{n}$.. Given that asymptotically, the Hermite polynomials grow factorially
as $H_{n}\sim\Gamma(n/2)$ (apart from the finite number of points where they
vanish), the gradient expansion is an asymptotic expansion. The effective
coupling constant of this asymptotic series is $\tau_{R}/(\sigma\gamma)$. A
well known property of an asymptotic series that grow as $g^{n}\Gamma(n/2)$ is
that it starts to diverge at order $n^{*}\sim 2/g^{2}$. According to the
optimal truncation procedure a la Poincaré, the gradient expansion can be
directly summed up to $n^{*}$ terms before the series start to diverge. This
result implies that when the system size is comparable to mean-free-path,
namely $\sigma=c\tau_{R}/\gamma$ for some constant $c$, we get $n^{*}=1$ and
the optimal truncation breaks down.

Furthermore in the ultra-relativistic limit, $v\rightarrow 1$, the gradient
expansion coefficients significantly simplify; $c_{n}=2^{n-1}+{\cal
O}(\gamma^{-1})$ and we find the asymptotic behavior of the gradient expansion
to be

$J_{D}\sim-\sum_{n}(-1)^{n}\left(\frac{2\sqrt{2}\tau_{R}}{\gamma\sigma}\right)^{n}\Gamma\left(\frac{n+1}{2}\right)e^{-\frac{(x-t)^{2}}{2\sigma}}\cos(\sqrt{2n}(x-t)-n\pi/2)\,.$
(69)

Therefore the breakdown of the optimal truncation occurs when
$\sigma=2\tau_{R}/\gamma$, confirming our heuristic argument earlier.

##  IV Stochastic dynamics

###  IV.1 Noise in the density frame

In this section we will add noise to the density frame diffusion equation and
study the stochastic dynamics of a boosted fluid. The dissipative strains in
the density frame take the form

$J_{D}^{i}=-T\sigma^{ij}\partial_{j}\hat{\mu}+\xi^{i}\,,$ (70)

where $\xi^{i}$ is the noise. For the stochastic process to equilibrate to the
probability distribution determined by the entropy of the system,

$P[N]\propto\exp\left(\mathcal{S}[N]\right)\propto\exp\left(-\int{\rm
d}^{3}x\,\frac{\beta^{0}}{2\chi^{00}}N^{2}\right)\,,$ (71)

the noise must respect the fluctuation-dissipation theorem

$\left\langle\xi^{i}(t,{\bf x})\,\xi^{j}(t^{\prime},{\bf
x}^{\prime})\right\rangle=2T\sigma^{ij}\,\delta^{3}({\bf x}-{\bf
x}^{\prime})\,\delta(t-t^{\prime})\,.$ (72)

The form of the noise matrix in the density frame can also be found by
algebraically manipulating the current in the Landau frame. In the Landau
frame, we have

$J^{\mu}=n_{{\rm\scriptscriptstyle LF}}u^{\mu}+j_{\rm\scriptscriptstyle
D,LF}^{\mu}+\xi^{\mu}_{\rm\scriptscriptstyle LF}\,,$ (73)

where the noise is orthogonal to $u^{\mu}$ and satisfies

$\left\langle\xi^{\mu}_{\rm\scriptscriptstyle
LF}\,\xi^{\nu}_{\rm\scriptscriptstyle
LF}\right\rangle=2T\sigma\,\Delta^{\mu\nu}\delta^{4}(x-y)\,.$ (74)

Rearranging the Landau frame variables into the density frame form, we find

$J^{\mu}=(N,Nv^{i}+J_{D}^{i}+\xi^{i})\,,$ (75)

where the density frame noise is

$\xi^{i}=\xi_{\rm\scriptscriptstyle LF}^{i}-v^{i}\xi_{\rm\scriptscriptstyle
LF}^{0}=(\Delta^{i}_{\phantom{\nu}\alpha}-v^{i}\Delta^{0}_{\phantom{\nu}\alpha})\,\xi^{\alpha}_{\rm\scriptscriptstyle
LF}\,.$ (76)

Computing the covariance of the density frame noise using eq. (74) yields

$\displaystyle\left\langle\xi^{i}(x)\xi^{j}(y)\right\rangle=$
$\displaystyle(\Delta^{i}_{\phantom{\nu}\alpha}-v^{i}\Delta^{0}_{\phantom{\nu}\alpha})(\Delta^{j}_{\phantom{\nu}\beta}-v^{j}\Delta^{0}_{\phantom{\nu}\beta})\,2T\sigma\Delta^{\alpha\beta}\delta^{4}(x-y)\,,$
(77) $\displaystyle=$ $\displaystyle 2T\sigma^{ij}\delta^{4}(x-y)\,.$ (78)

Thus the form of the noise in the density frame can be straightforwardly found
from the Landau frame definitions.

###  IV.2 The Metropolis algorithm for stochastic equations

Next we discuss how the dissipative stochastic dynamics can be simulated using
a Metropolis algorithm, rather than directly discretizing the Langevin
dynamics. As discussed in the introduction, the approach has the advantage
that detailed balance is maintained irrespective of the time step $\Delta t$.

To understand the method, consider the one-dimensional Brownian motion of a
particle in a potential $U(q)$. The Brownian particle evolves in phase-space
as

$\displaystyle\frac{dq}{dt}+\left\\{\mathcal{H},q\right\\}=$ $\displaystyle
0\,,$ (79a) $\displaystyle\frac{dp}{dt}+\left\\{\mathcal{H},p\right\\}=$
$\displaystyle-\eta\left(\frac{\partial\mathcal{H}}{\partial
p}\right)+\xi\,,\qquad\qquad\left\langle\xi(t)\xi(t^{\prime})\right\rangle=2T\eta\,\delta(t-t^{\prime})\,,$
(79b)

where the free energy of the particle with momentum $p$ and position $q$ is

$\mathcal{H}(q,p)=\frac{p^{2}}{2m}+U(q)\,.$ (80)

Here $\eta$ is the drag coefficient and the drag force is proportional to the
velocity, $\partial\mathcal{H}/\partial p=v$, i.e. the variable
thermodynamically conjugate to $p$. The noise is chosen so that the system
evolves to the equilibrium probability distribution $P(q,p)\propto
e^{-\beta\,\mathcal{H}(q,p)}$.

A natural way to simulate the dynamics is to use operator splitting, first
setting the right hand side of eq. (79) to zero and taking a symplectic step.
Ideally, this step should be done with a symplectic integrator which preserves
the phase-space volume. The symplectic update is followed by a dissipative
step such as the Metropolis update discussed below, which respects detailed
balance. Together, the two steps correctly evolve eq. (79) over a time $\Delta
t$.

In the Metropolis update algorithm over a time interval $\Delta t$, one makes
a proposal

$p\rightarrow p+\Delta p\,,\qquad\qquad\Delta p=\sqrt{2T\eta\Delta
t}\,\mathfrak{e}\,,$ (81)

where $\mathfrak{e}$ is a random number of variance one. Then the change in
free energy in the proposed step is

$\beta\Delta\mathcal{H}=\beta\left(\mathcal{H}(p+\Delta
p)-\mathcal{H}(p)\right)\simeq\beta\frac{\partial\mathcal{H}}{\partial
p}\,\Delta p\,.$ (82)

In a Metropolis approach if $\Delta\mathcal{H}$ is negative then the proposal
is accepted; if $\Delta\mathcal{H}$ is positive then the proposal is accepted
with probability $e^{-\beta\Delta\mathcal{H}}\simeq 1-\beta\Delta\mathcal{H}$.
Because of the asymmetry between gain and loss rates, the particle will
experience drag in addition to the noise added in eq. (81). It is
straightforward to see that the mean momentum transfer $\Delta p$ from the
Metropolis step is

$\left\langle\Delta p\right\rangle\simeq\int_{-\infty}^{\infty}{\rm
d}\mathfrak{e}P(\mathfrak{e})\left[\Theta(-\Delta\mathcal{H})+\Theta(\Delta\mathcal{H})\left(1-\beta\Delta\mathcal{H}\right)\right]\Delta
p\simeq-\eta\frac{\partial H}{\partial p}\Delta t\,,$ (83)

reproducing the mean drag in the Langevin equations of motion. A similar
computation shows that $\left\langle(\Delta p)^{2}\right\rangle=2T\eta\,\Delta
t$, indicating that the Metropolis algorithm correctly reproduces the drag and
noise of the Langevin evolution.

###  IV.3 The advection diffusion equation from the Metropolis algorithm

Now we will show how the same Metropolis algorithm can be used to simulate the
Langevin updates for the relativistic advection-diffusion equation in the
density frame. The continuum equation we would like to solve is

$\partial_{t}N+\partial_{i}(Nv^{i})+\partial_{i}(J_{D}^{i}+\xi^{i})=0\,,$ (84)

where the dissipative part is

$J_{D}^{i}=-T\sigma\left(\delta^{ij}-v^{i}v^{j}\right)\partial_{j}\hat{\mu}\,.$
(85)

For simplicity we will limit the discussion to two spatial dimensions.

As in the Brownian motion example an operator splitting approach is adopted,
we first solve the advection equation

$\partial_{t}N+\partial_{i}(Nv^{i})=0\,,$ (86)

which captures the symplectic dynamics in this case. For the advective step we
adopt the Kurganov-Tadmor (KT) central scheme using a second order spatial
discretization [44, 45, 46]. Given the stochastic nature of the simulation, we
turned off limiters such as min-mod or WENO based limiters. This choice and
possible alternatives should be reexamined in the future. For the time
integration we use a second order Total Variation Diminishing (TVD) Runge
Kutta method implemented in the PETSc library [47, 48]. We use a fixed time
step $\Delta t=0.5\,\Delta x/c$ in the numerical experiments presented below.
We note that both advection and ideal relativistic hydrodynamics have a
symplectic structure, which can be derived from Poisson brackets between the
conserved charges [49]. However, the KT scheme with the TVD time
discretization is not a symplectic integrator. It would be interesting to
explore symplectic integrators for ideal hydrodynamics when physical
dissipation (which naturally leads to TVD property) is incorporated in
subsequent steps.

After taking an advective step, we propose random transfers of charge between
the fluid cells with appropriate variances. The charge transfers are accepted
or rejected according to the statistical weight, $\exp(\Delta\mathcal{S})$.
The procedure parallels the Brownian motion example of the previous subsection
and reproduces the mean diffusive current as well as the noise. In the next
paragraphs we will explicitly list the algorithm and verify this claim.

The simulation is discretized on a two dimensional lattice with a finite
volume discretization and fixed lattice spacing. The lattice metric is

$ds^{2}=a_{x}^{2}\left(\Delta I_{x}\right)^{2}+a_{y}^{2}\left(\Delta
I_{y}\right)^{2}\,,$ (87)

with $a_{x}$ and $a_{y}$ being lattice spacing, and
$I=\left(I_{x},I_{y}\right)$ denotes the (integer) lattice coordinates. The
volume of a fluid cell is $V_{0}=\sqrt{g}=a_{x}a_{y}$ and the charge in the
$I^{\rm th}$ cell is $\mathcal{N}_{I}=V_{0}N_{I}$. The entropy of the system
takes the form

$\mathcal{S}[\mathcal{N}]=\mathcal{S}_{1}-\sum_{I}\frac{{\mathcal{N}}^{\,2}_{I}}{2T\chi
u^{0}V_{0}}\,,$ (88)

and derivatives of $\mathcal{S}$ with respect to $\mathcal{N}_{I}$ determine
the chemical potential

$\frac{\partial\mathcal{S}}{\partial\mathcal{N}_{I}}=-\hat{\mu}_{I}\,,$ (89)

where $\hat{\mu}_{I}=\mathcal{N}_{I}/T\chi u^{0}V_{0}$.

The layout of the grid is shown in Fig. 5. We imagine a stochastic current
living at the corner of the computational cells with covariance

$\left\langle\xi^{i}\xi^{j}\right\rangle=\frac{2T\sigma}{\Delta
tV_{0}}\left(\delta^{ij}-v^{i}v^{j}\right)\,,$ (90)

where $\Delta t$ is the time step.

Figure 5: Discretization of Metropolis proposals. The noise $\xi^{i}$ “lives”
on the corner of the computational cells, $A$, $B$, $C$, $D$. The charge
transfers between cells are given in eq. (95).

One way to produce this noise is to generate noises parallel and perpendicular
to the fluid velocity, $\xi_{\parallel}$ and $\xi_{\perp}$, with variances

$\displaystyle\xi_{\parallel}=$ $\displaystyle\sqrt{\frac{2T\sigma}{\Delta
tV_{0}}(1-v^{2})}\;{\mathfrak{e}}_{\parallel}\,,$ (91)
$\displaystyle\xi_{\perp}=$ $\displaystyle\sqrt{\frac{2T\sigma}{\Delta
tV_{0}}}\;{\mathfrak{e}}_{\perp}\,,$ (92)

where ${\mathfrak{e}}_{\parallel}$ and ${\mathfrak{e}}_{\perp}$ are random
numbers with zero mean and unit variance. Then a rotation gives a proposal
with the expected variance in eq. (90).

The proposed charge transfer in the $x$ and $y$ directions are

$Q^{i}=\xi^{i}A_{i}\Delta t=\xi^{i}V_{0}\Delta t/a_{i}\,.\qquad\mbox{(no
sum)}$ (93)

Here and for the rest of this section no sum is implied by repeated indices.
The variance of the proposed charge transfers is

$\left\langle Q^{i}Q^{j}\right\rangle=\frac{2T\sigma\Delta
tV_{0}}{a_{i}a_{j}}\left(\delta^{ij}-v^{i}v^{j}\right)\,,$ (94)

leading to the proposed updates for the cells $A$, $B$, $C$, and $D$ (see Fig.
5)

$\displaystyle\mathcal{N}_{A}\rightarrow$
$\displaystyle\mathcal{N}_{A}+\Delta\mathcal{N}_{A}\equiv\mathcal{N}_{A}-\frac{Q^{x}}{2}+\frac{Q^{y}}{2}\,,$
(95a) $\displaystyle\mathcal{N}_{B}\rightarrow$
$\displaystyle\mathcal{N}_{B}+\Delta\mathcal{N}_{B}\equiv\mathcal{N}_{B}+\frac{Q^{x}}{2}+\frac{Q^{y}}{2}\,,$
(95b) $\displaystyle\mathcal{N}_{C}\rightarrow$
$\displaystyle\mathcal{N}_{C}+\Delta\mathcal{N}_{C}\equiv\mathcal{N}_{C}+\frac{Q^{x}}{2}-\frac{Q^{y}}{2}\,,$
(95c) $\displaystyle\mathcal{N}_{D}\rightarrow$
$\displaystyle\mathcal{N}_{D}+\Delta\mathcal{N}_{D}\equiv\mathcal{N}_{D}-\frac{Q^{x}}{2}-\frac{Q^{y}}{2}\,.$
(95d)

Then we compute the change in the entropy for a proposal, using (89), which
reads

$\displaystyle\Delta\mathcal{S}=$
$\displaystyle\sum_{U=A,B,C,D}\mathcal{S}[\mathcal{N}_{U}+\Delta\mathcal{N}_{U}]-S[\mathcal{N}_{U}]\,,$
(96) $\displaystyle\simeq$
$\displaystyle-\frac{1}{2}(\hat{\mu}_{B}+\hat{\mu}_{C}-\hat{\mu}_{A}-\hat{\mu}_{D})\,Q^{x}-\frac{1}{2}(\hat{\mu}_{A}+\hat{\mu}_{B}-\hat{\mu}_{C}-\hat{\mu}_{D})\,Q^{y}\,.$
(97)

Formally, one can also write this as

$\Delta\mathcal{S}\simeq-\sum_{i}\partial_{i}\hat{\mu}\,a_{i}Q^{i}\,,$ (98)

where it is understood that, for instance,
$\partial_{x}\hat{\mu}\equiv(\hat{\mu}_{B}+\hat{\mu}_{C}-\hat{\mu}_{A}-\hat{\mu}_{D})/2a_{x}$.
Then the probability of accepting the proposed update is

$P_{\rm
accept}(\xi)=\theta(\Delta\mathcal{S})+\theta\left(-\Delta\mathcal{S}\right)\left(e^{\Delta\mathcal{S}}-1\right)\simeq
1+\theta(-\Delta\mathcal{S})\,\Delta\mathcal{S}\,.$ (99)

Thus, the mean charge transfer in the Metropolis step is given by

$\displaystyle\left\langle Q^{i}\right\rangle\simeq$ $\displaystyle\int{\rm
d}^{2}\xi\;P(\xi)\;Q^{i}\,\left(1+\theta\left(-\Delta\mathcal{S}\right)\Delta\mathcal{S}\right)\,,$
$\displaystyle\simeq$
$\displaystyle-\frac{1}{2}\sum_{j}\partial_{j}\hat{\mu}\,a_{j}\left\langle
Q^{i}Q^{j}\right\rangle\,,$ $\displaystyle=$ $\displaystyle-\frac{T\sigma
V_{0}\Delta t}{a_{i}}\sum_{j}(\delta^{ij}-v^{i}v^{j})\partial_{j}\hat{\mu}\,.$
(100)

The factor of a half in the second line arises because we are only integrating
over proposals where $\Delta\mathcal{S}<0$. Dividing by $A_{i}\Delta
t=V_{0}\Delta t/a_{i}$, we find the mean current from the metropolis step that
is consistent with the expected form of the diffusive current in the density
frame

$\displaystyle\left\langle J^{i}_{D}\right\rangle\simeq$
$\displaystyle-T\sigma\sum_{j}(\delta^{ij}-v^{i}v^{j})\partial_{j}\hat{\mu}\,.$
(101)

To summarize, the full procedure consists of an ideal advective step, followed
by a diffusive step. In the diffusive step we step over the lattice by two’s,
first updating the group of cells $A$, $B$, $C$, $D$ and then proceeding to
the next independent group of four cells. The updates are independent of each
other and can be done in any order. This covers one quarter of the lattice. We
then loop through the remaining three corners in a similar way to complete the
diffusive steps. In fact, to eliminate the potential bias, the order of the
four corners which are updated is randomly shuffled in each diffusive step.

There are many choices and questions here which can be studied in future work.
For instance, it is not necessary to take one diffusive step per advective
step. In fact, in order to be closer to the Langevin limit we take $400$
diffusive steps per advective step. This is quite a large number and
guarantees that the rejection probability $r$ approaches zero. It is only in
the asymptotic limit $r\rightarrow 0$ that the Metropolis updates are fully
equivalent to the Langevin simulations. In our numerical experiments
$r=0.003$. It would be helpful to explore the approach to the Langevin limit
in greater detail, which allow for better algorithms that capture the
interplay between the symplectic and dissipative dynamics.

###  IV.4 Equilibrium correlation functions

As a first test of the stochastic dynamics we will compute the correlation
function of charges advecting and diffusing in a two dimensional fluid moving
with a fixed velocity, ${\bf v}=v\,(\cos\theta,\sin\theta)$. In our test case
we treated a fluid moving at $v=0.8\,{\rm c}$ at an angle of $\theta=\pi/6$
and used a lattice of $L^{2}=128^{2}\,a^{2}$ where $a=a_{x}=a_{y}$ is the
lattice spacing.

####  IV.4.1 Physical considerations

Three dimensionful parameters in the simulation can be set to unity, setting
our units of space, time, and energy. We choose the lattice length and the
speed of light to be one, $a=c=1$. The variance of the charge in a fluid cell
$\left\langle\mathcal{N}^{\,2}\right\rangle=T\chi u^{0}a^{d}$ may also be set
to unity, where $d=2$ is the number of spatial dimensions. All physical
quantities can be expressed in terms of $a$, $c$ and $T\chi u^{0}$.

The “mean free path” of the system $\ell_{\rm mfp}$ is defined through the
diffusion coefficient $D\equiv\tfrac{1}{3}\ell_{\rm mfp}c$. The mean free path
in units of the lattice spacing $\ell_{\rm mfp}/a$ is a dimensionless
parameter, which can only be fixed through physical considerations. We are
only interested in modes where $k\ell_{\rm mfp}\ll 1$, as wave-numbers of
order $1/\ell_{\rm mfp}$ have been integrated out of the hydrodynamic
effective theory. Thus, we set $a=\ell$ which cuts off the wave numbers in the
simulation at a reasonable value. Hence, long wavelength modes on the lattice
are well described by the continuum description and the diffusion equation,
while modes of order the lattice spacing are neither resolved nor adequately
described by the diffusion equation.

In modeling the charge fluctuations we have ignored the discrete nature of the
charge carriers, neglecting shot noise. Consider a field theory with a finite
number of fields and assume that charge susceptibility is of order $T\chi\sim
e^{2}(T/\hbar c)^{d}$, where $e$ is the elementary charge. This is the case
for the electric charge susceptibility of QCD at high temperatures. The
variance of the charge within a fluid cell in units of the elementary charge
$e$ is

$\frac{a^{d}}{e^{2}}\left\langle\delta
N^{2}\right\rangle=\frac{a^{d}}{e^{2}}T\chi u^{0}\sim\left(\frac{aT}{\hbar
c}\right)^{d}\,.$ (102)

If the theory is weakly coupled $\ell_{\rm mfp}\gg(\hbar c/T)$, then this
variance is large for $a\sim\ell_{\rm mfp}$, and it is appropriate to treat
the charge and charge fluctuations using continuous variables even for short
wavelengths, $ka\sim 1$. When simulating weakly coupled fluids, errors will
still arise from the space-time discretization and from using the diffusion
equation for $k\ell_{\rm mfp}\sim 1$, but not arise from treating charge as a
continuous variable. In a strongly coupled field theories where $\ell_{\rm
mfp}\sim\hbar c/T$ and $\chi\ell^{d}_{\rm mfp}/e^{2}\sim 1$, the variance of a
fluid cell in units of $e^{2}$ is of order unity and the discretized
hydrodynamic theory does not capture the quantized charge fluctuations (or
shot noise) at the scale of the mean free path. This error is the same order
of magnitude as the discretization and modeling errors made in the weakly
coupled case. (In strongly coupled, but large $N_{c}$ field theories, the
susceptibility is of order $T\chi\sim N_{c}\,e^{2}(T/\hbar c)^{d}$ and shot
noise is always negligible.) Finally, in theories where the susceptibility is
very small $\chi\ell_{\rm mfp}^{d}/e^{2}\ll 1$, it should be possible to
include shot noise systematically into the hydrodynamic description by making
discrete Poissonian proposals for the charge transfers between fluid cells.
However, we have adopted continuous charge transfers here and leave this
regime for future work.

####  IV.4.2 Numerical results

The density-density correlation function in the simulation is

$C_{NN}(t-t^{\prime},{\bf k})\equiv\int{\rm d}^{3}x\,e^{-i{\bf k}\cdot{\bf
x}}\left\langle N(t,{\bf x})N(t^{\prime},{\bf
0})\right\rangle\equiv\frac{1}{V}\left\langle N(t,{\bf k})N(t^{\prime},-{\bf
k})\right\rangle\,.$ (103)

From the density frame equation of motion

$\partial_{t}N+\partial_{i}(Nv^{i})-\partial_{i}\left(D^{ij}\partial_{j}N\right)+\partial_{i}\xi^{i}=0\,,$
(104)

the expected correlation function can be computed straightforwardly. Indeed,
one can solve for $N(t,{\bf k})$ in terms of $\xi(t,{\bf k})$

$N(t,{\bf k})=\int^{t}_{-\infty}dt^{\prime}\,(-ik_{m}\xi^{m}(t^{\prime},{\bf
k}))\,e^{-i{\bf v}\cdot{\bf
k}\,(t-t^{\prime})}e^{-D^{ij}k_{i}k_{j}\,(t-t^{\prime})}\,.$ (105)

Squaring this expression and averaging over the noise determines the expected
form of the density frame correlation function

$C_{NN}(t-t^{\prime},{\bf k})=T\chi u^{0}\cos({\bf v}\cdot{\bf
k}\,(t-t^{\prime}))\exp(-D^{ij}k_{i}k_{j}\,|t-t^{\prime}|)\,.$ (106)

Figure 6: Comparison between the expected (106) and simulated density-density
correlation functions in the density frame as a function of time (in lattice
units) for various wave numbers. The fluid is moving with velocity $v=0.8c$ at
an angle of $30^{\circ}$ above the $x$ axis. The left plots have wave numbers
$(k_{x},k_{y})=(k_{x},0)$ for various values of $k_{x}$. Similarly, the right
plots have $(k_{x},k_{y})=(0,k_{y})$ for various values of $k_{y}$.

Figure 6 shows a comparison between the expected (106) and simulated
correlation functions for different wave numbers in the fluid. Examining eq.
(106) we see that the oscillations reflect the advection of a sinusoidal wave,
while the exponential decay is controlled by the diffusion matrix $D^{ij}$.
Naturally, the longest wavelengths (smallest wave-numbers) have the slowest
decay, and this is clearly seen from the trends in Fig. 6. The diffusion
matrix takes the form

$D^{ij}=\frac{D}{\gamma^{3}}\hat{v}^{i}\hat{v}^{j}+\frac{D}{\gamma}\left(\delta^{ij}-\hat{v}^{i}\hat{v}^{j}\right)\,,$
(107)

and thus Fourier modes which are parallel to the flow velocity will decay
slowly relative to the transverse modes, i.e. at rates $\sim
Dk^{2}/\gamma^{3}$ and $\sim Dk^{2}/\gamma$ respectively. The curves with
$k_{y}=0$ and $k_{x}$ finite (the left hand side of Fig. 6) are more aligned
with the flow than the curves with $k_{x}=0$ and $k_{y}$ finite (the right
hand side of Fig. 6), and thus the right plots show a stronger exponential
decay, confirming this expectation.

##  V Discussion

We have shown how the Metropolis algorithm can be adapted to simulate
stochastic relativistic advection diffusion equation in two dimensions. In a
companion paper we will describe how the framework can be further extended to
stochastic viscous hydrodynamics in general relativity.

The algorithm is simple:  (i) take an ideal advective step,  (ii) make a
proposal to randomly transfer the conserved charges between computational
cells, and finally  (iii) accept or reject the proposed transfers based on how
they change the entropy of the fluid in a Metropolis-Hastings accept-reject
step.  The average charge transfer reproduces the mean diffusive current.

The continuum formulation of the stochastic process is not Lorentz covariant.
But, the equations of motion are invariant under Lorentz transformations
followed by a reparametrization of the hydrodynamic fields consistent with the
derivative expansion. Indeed, to describe the stochastic dynamics we have
adopted a formulation of hydrodynamics developed to describe hydrodynamics
without boosts [34, 35, 36]. In particular we made considerable use of
“density frame” of [36].

A notable feature of the density frame is that the charge in a fluid cell
($J^{0}$ in this case) is sufficient to determine the associated chemical
potential, $\mu\equiv J^{0}/\chi u^{0}$. By contrast, in the Landau frame the
charge and the three current $J^{i}$ are both needed,
$\mu_{\rm\scriptscriptstyle LF}=-u_{\mu}J^{\mu}/\chi$. Because of this
feature, the equations of motion are first order in time, and do not need any
auxiliary variables such as the diffusive current or, in full hydrodynamics,
the strain tensor $\pi^{\mu\nu}$. The approach is numerically stable for the
diffusion equation and is expected to be stable for full hydrodynamics,
providing a practical way to simulate both stochastic and noise-averaged
relativistic fluids in heavy ion collisions. The equations of motion do not
obey Lorentz causality, since the equations do not model the high momentum
modes which lie outside of the validity of hydrodynamics. These high momentum
modes are essential for Lorentz causality [31, 32, 33, 43, 50]. In the density
frame approach these Fourier modes are simply increasingly damped as
$\sim\exp(-D^{ij}k_{i}k_{j}t)$ and do not affect the long wavelength dynamics.

In spite of these “problems”, the density frame dynamics describes well the
diffusion of charge in a highly boosted fluid with $\gamma=10$. Even in the
regime where the mean free path becomes comparable to the system size, the
leading order density frame predictions remain qualitatively well behaved. For
a specific test case described in Sec. III, we were able to work out all
higher order terms of the density frame gradient expansion, which in the
regime of validity of the hydrodynamics, systematically improve the leading
order results. A resummation of this expansion reproduces the causality and
Lorentz covariant structure of the dispersion curve of the underlying kinetic
model, even though Lorentz causality is violated at any finite order in the
gradient expansion.

The next step in this project is to use the algorithm developed here to
simulate hydrodynamics in general coordinates and to develop the Metropolis
updates into a practical tool for stochastic hydrodynamic simulations of heavy
ion collisions. Indeed, the paradigm of the current paper should work for full
hydrodynamics by implementing the following steps:  (i) first take a step with
ideal hydrodynamics;  (ii) make a proposal to randomly transfer spatial
momentum between the fluid cells;  (iii) accept or reject the proposal using
the entropy as a weight. In general coordinates the only complication is that
the momentum transfers must be parallel transported from the cell-faces to the
cell-centers before applying the accept-reject criterion. It is hoped that the
Metropolis algorithm for stochastic hydrodynamics will be robust and
effective, yielding a significant advance in the modeling of the quark-gluon
plasma created in heavy ion collisions.

###### Acknowledgements.

G.B. is supported by the National Science Foundation CAREER Award PHY-2143149.
J.B. and D.T. are supported by the U.S. Department of Energy, Office of
Science, Office of Nuclear Physics, grant No. DE-FG-02-08ER41450. R.S.
acknowledges the support of Polish NAWA Bekker program No.
BPN/BEK/2021/1/00342. Finally, we are grateful for the stimulating atmosphere
at the “The Many Faces of Relativistic Fluid Dynamics” program at the Kavli
Institute for Theoretical Physics.

##  Appendix A Green functions for the kinetic model

In this section we present the Green functions associated with the kinetic
model we analyzed in Section III. Let us consider the static case first. As
usual, given initial data at some initial time which we set to be $t=0$ the
solution to the kinetic equation at a later time $t>0$ is given by propagating
the initial data, $(n_{0}(x^{\prime}),j_{0}(x^{\prime}))$ via the retarded
Green function

$\left(\begin{matrix}n(t,x)\\\ j(t,x)\end{matrix}\right)=\int
dx^{\prime}G_{R}(t,x-x^{\prime})\left(\begin{matrix}n_{0}(x^{\prime})\\\
j_{0}(x^{\prime})\end{matrix}\right)\,,$ (108)

where the Green function $G_{R}$ is a $2\times 2$ matrix which satisfies

$\left(\begin{matrix}\partial_{t}&&\partial_{x}\\\
\partial_{x}&&2\lambda+\partial_{t}\end{matrix}\right)G_{R}(t-t^{\prime},x-x^{\prime})=\mathbbm{1}_{2\times
2}\delta(t-t^{\prime})\delta(x-x^{\prime})\,.$ (109)

Here $\lambda=1/(2\tau_{R})$ and we set $c_{s}=1$ for simplicity. By Fourier
transforming and with the help of the integral

$\int\frac{dk}{2\pi}\frac{\sin\left(\sqrt{k^{2}-\lambda^{2}}t\right)}{\sqrt{k^{2}-\lambda^{2}}}e^{-ikx}=\frac{1}{2}\Theta(t^{2}-x^{2})I_{0}(\lambda\sqrt{t^{2}-x^{2}}),$
(110)

where $\Theta$ is the Heaviside Theta function and $I_{0}$ is the modified
Bessel function of the first kind, we obtain

$\displaystyle G_{R}(t,x)=\frac{\lambda}{2}e^{-\lambda
t}\Theta(t^{2}-x^{2})\left(\begin{matrix}G_{nn}&&G_{nj}\\\
G_{jn}&&G_{jj}\end{matrix}\right)+\frac{\lambda}{2}e^{-\lambda
t}\left(\begin{matrix}\delta_{+}+\delta_{-}&&-\delta_{+}+\delta_{-}\\\
-\delta_{+}+\delta_{-}&&\delta_{+}+\delta_{-}\end{matrix}\right)\,.$ (111)

Here

$\displaystyle\delta_{\pm}$ $\displaystyle=$ $\displaystyle\delta(x\pm
t),\quad\tau=\sqrt{t^{2}-x^{2}}\,,$ (112) $\displaystyle G_{nn}(t,x)$
$\displaystyle=$ $\displaystyle\frac{\lambda}{2}e^{-\lambda
t}\left(\frac{t}{\tau}I_{1}(\lambda\tau)+I_{0}(\lambda\tau)\right)\,,$ (113)
$\displaystyle G_{nj}(t,x)$ $\displaystyle=$ $\displaystyle
G_{jn}(t,x)=\frac{\lambda}{2}e^{-\lambda
t}\left(\frac{x}{\tau}I_{1}(\lambda\tau)\right)\,,$ (114) $\displaystyle
G_{jj}(t,x)$ $\displaystyle=$ $\displaystyle\frac{\lambda}{2}e^{-\lambda
t}\left(\frac{t}{\tau}I_{1}(\lambda\tau)-I_{0}(\lambda\tau)\right)\,.$ (115)

Integrating the singular part of the Green function explicitly we obtain the

$\displaystyle n(t,x)$ $\displaystyle=$ $\displaystyle\frac{e^{-\lambda
t}}{2}\left(n_{0}(x-t)+n_{0}(x+t)+j_{0}(x-t)-j_{0}(x+t)\right)$ (116)
$\displaystyle+\int_{x-t}^{x+t}dx^{\prime}\left(G_{nn}(t,x-x^{\prime})n_{0}(x^{\prime})+G_{nj}(t,x-x^{\prime})j_{0}(x^{\prime})\right)\,,$
$\displaystyle j(t,x)$ $\displaystyle=$ $\displaystyle\frac{e^{-\lambda
t}}{2}\left(n_{0}(x-t)-n_{0}(x+t)+j_{0}(x-t)+j_{0}(x+t)\right)$ (117)
$\displaystyle+\int_{x-t}^{x+t}dx^{\prime}\left(G_{nj}(t,x-x^{\prime})n_{0}(x^{\prime})+G_{jj}(t,x-x^{\prime})j_{0}(x^{\prime})\right)\,.$

Since the system is Lorentz covariant, the Green functions for the moving
fluid can be obtained by a Lorentz boost. We consider the initial value
problem where the initial conditions are still given in the lab frame $t=0$.
The charge density and current given in Eq.(53) are therefore given by

$\displaystyle N(t,x)$ $\displaystyle=$
$\displaystyle\frac{e^{-\lambda\kappa_{-}t}}{2}\left(N_{0}(x-t)+J_{0}(x-t)\right)+\frac{e^{-\lambda\kappa_{+}t}}{2}\left(N_{0}(x+t)-J_{0}(x+t)\right)$
(118)
$\displaystyle+\gamma\int_{x-t}^{x+t}dx^{\prime}\left[\left(G_{nn}(\tilde{t}-\gamma
vx,\gamma x-\tilde{x}^{\prime})+vG_{nj}(\tilde{t}-\gamma vx,\gamma
x-\tilde{x}^{\prime})\right)N_{0}(x^{\prime})\right.$
$\displaystyle\quad\quad\quad\;\left.+\left(G_{nj}(\tilde{t}-\gamma vx,\gamma
x-\tilde{x}^{\prime})+vG_{jj}(\tilde{t}-\gamma vx,\gamma
x-\tilde{x}^{\prime})\right)J_{0}(x^{\prime})\right]\,,$ $\displaystyle
J(t,x)$ $\displaystyle=$
$\displaystyle\frac{e^{-\lambda\kappa_{-}t}}{2}\left(N_{0}(x-t)+J_{0}(x-t)\right)-\frac{e^{-\lambda\kappa_{+}t}}{2}\left(N_{0}(x+t)-J_{0}(x+t)\right)$
(119)
$\displaystyle+\gamma\int_{x-t}^{x+t}dx^{\prime}\left[\left(G_{nj}(\tilde{t}-\gamma
vx,\gamma x-\tilde{x}^{\prime})+vG_{nn}(\tilde{t}-\gamma vx,\gamma
x-\tilde{x}^{\prime})\right)N_{0}(x^{\prime})\right.$
$\displaystyle\quad\quad\quad\;\left.+\left(G_{jj}(\tilde{t}-\gamma vx,\gamma
x-\tilde{x}^{\prime})+vG_{nj}(\tilde{t}-\gamma vx,\gamma
x-\tilde{x}^{\prime})\right)J_{0}(x^{\prime})\right]\,.$

where $\kappa_{\pm}=\sqrt{(1\pm v)/(1\mp v)}$,
$\tilde{t}=\gamma(t+vx^{\prime})$ and
$\tilde{x}^{\prime}=\gamma(x^{\prime}+vt)$.

## References

  * [1] U. Heinz and R. Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part. Sci. 63 (2013) 123–151, arXiv:1301.2826 [nucl-th].
  * [2] J. E. Bernhard, J. S. Moreland, and S. A. Bass, “Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma,” Nature Phys. 15 no. 11, (2019) 1113–1117.
  * [3] G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, “Transverse momentum differential global analysis of heavy-ion collisions,” Phys. Rev. Lett. 126 (May, 2021) 202301. https://link.aps.org/doi/10.1103/PhysRevLett.126.202301.
  * [4] JETSCAPE Collaboration, D. Everett et al., “Multisystem bayesian constraints on the transport coefficients of qcd matter,” Phys. Rev. C 103 (May, 2021) 054904. https://link.aps.org/doi/10.1103/PhysRevC.103.054904.
  * [5] M. Heffernan, Quantification of the Quark-Gluon Plasma with statistical learning.  PhD thesis, McGill U., 4, 2023.  https://escholarship.mcgill.ca/concern/theses/0r967886q.
  * [6] M. Arslandok et al., “Hot QCD White Paper,” arXiv:2303.17254 [nucl-ex].
  * [7] R. D. Pisarski and F. Wilczek, “Remarks on the Chiral Phase Transition in Chromodynamics,” Phys. Rev. D 29 (1984) 338–341.
  * [8] K. Rajagopal and F. Wilczek, “Static and dynamic critical phenomena at a second order QCD phase transition,” Nucl. Phys. B 399 (1993) 395–425, arXiv:hep-ph/9210253.
  * [9] HotQCD Collaboration, H. T. Ding et al., “Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD,” Phys. Rev. Lett. 123 no. 6, (2019) 062002, arXiv:1903.04801 [hep-lat].
  * [10] O. Kaczmarek, F. Karsch, A. Lahiri, L. Mazur, and C. Schmidt, “QCD phase transition in the chiral limit,”  3, 2020.  arXiv:2003.07920 [hep-lat].
  * [11] M. A. Stephanov, “QCD Phase Diagram and the Critical Point,” Prog. Theor. Phys. Suppl. 153 (2004) 139–156, arXiv:hep-ph/0402115.
  * [12] L. Du, A. Sorensen, and M. Stephanov, “The QCD phase diagram and Beam Energy Scan physics: a theory overview,”  2, 2024.  arXiv:2402.10183 [nucl-th].
  * [13] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics.  Cambridge University Press, 4 ed., 2014.
  * [14] A. Florio, E. Grossi, A. Soloviev, and D. Teaney, “Dynamics of the $O(4)$ critical point in QCD,” Phys. Rev. D 105 no. 5, (2022) 054512, arXiv:2111.03640 [hep-lat].
  * [15] P. J. Rossky, J. D. Doll, and H. L. Friedman, “Brownian dynamics as smart Monte Carlo simulation,” The Journal of Chemical Physics 69 no. 10, (11, 1978) 4628–4633. https://doi.org/10.1063/1.436415.
  * [16] G. D. Moore, “The Sphaleron rate: Bodeker’s leading log,” Nucl. Phys. B 568 (2000) 367–404, arXiv:hep-ph/9810313.
  * [17] P. B. Arnold, “Langevin equations with multiplicative noise: Resolution of time discretization ambiguities for equilibrium systems,” Phys. Rev. E 61 (2000) 6091–6098, arXiv:hep-ph/9912208.
  * [18] P. Gao, P. Glorioso, and H. Liu, “Ghostbusters: Unitarity and Causality of Non-equilibrium Effective Field Theories,” JHEP 03 (2020) 040, arXiv:1803.10778 [hep-th].
  * [19] A. Florio, E. Grossi, and D. Teaney, “Dynamics of the $O(4)$ critical point in QCD: critical pions and diffusion in Model G,” arXiv:2306.06887 [hep-lat].
  * [20] C. Chattopadhyay, J. Ott, T. Schaefer, and V. Skokov, “Dynamic scaling of order parameter fluctuations in model B,” Phys. Rev. D 108 no. 7, (2023) 074004, arXiv:2304.07279 [nucl-th].
  * [21] L. D. Landau and E. M. Lifshitz, Fluid Mechanics: Course of Theoretical Physics, vol. 6.  Elsevier, 2013.
  * [22] W. A. Hiscock and L. Lindblom, “Generic instabilities in first-order dissipative relativistic fluid theories,” Phys. Rev. D 31 (1985) 725–733.
  * [23] L. Gavassino, “Can We Make Sense of Dissipation without Causality?,” Phys. Rev. X 12 no. 4, (2022) 041001, arXiv:2111.05254 [gr-qc].
  * [24] W. Israel, “Nonstationary irreversible thermodynamics: A Causal relativistic theory,” Annals Phys. 100 (1976) 310–331.
  * [25] W. Israel and J. M. Stewart, “Transient relativistic thermodynamics and kinetic theory,” Annals Phys. 118 (1979) 341–372.
  * [26] R. P. Geroch and L. Lindblom, “Dissipative relativistic fluid theories of divergence type,” Phys. Rev. D 41 (1990) 1855.
  * [27] H. C. Öttinger, “Relativistic and nonrelativistic description of fluids with anisotropic heat conduction,” Physica A: Statistical Mechanics and its Applications 254 no. 3, (1998) 433–450. https://www.sciencedirect.com/science/article/pii/S0378437198000454.
  * [28] F. S. Bemfica, M. M. Disconzi, and J. Noronha, “Causality and existence of solutions of relativistic viscous fluid dynamics with gravity,” Phys. Rev. D 98 no. 10, (2018) 104064, arXiv:1708.06255 [gr-qc].
  * [29] P. Kovtun, “First-order relativistic hydrodynamics is stable,” JHEP 10 (2019) 034, arXiv:1907.08191 [hep-th].
  * [30] L. Lindblom, “The Relaxation effect in dissipative relativistic fluid theories,” Annals Phys. 247 (1996) 1, arXiv:gr-qc/9508058.
  * [31] M. P. Heller, A. Serantes, M. Spaliński, and B. Withers, “Rigorous Bounds on Transport from Causality,” Phys. Rev. Lett. 130 no. 26, (2023) 261601, arXiv:2212.07434 [hep-th].
  * [32] L. Gavassino, “Bounds on transport from hydrodynamic stability,” Physics Letters B 840 (2023) 137854. https://www.sciencedirect.com/science/article/pii/S0370269323001880.
  * [33] L. Gavassino, M. M. Disconzi, and J. Noronha, “Dispersion relations alone cannot guarantee causality,” arXiv:2307.05987 [hep-th].
  * [34] I. Novak, J. Sonner, and B. Withers, “Hydrodynamics without boosts,” JHEP 07 (2020) 165, arXiv:1911.02578 [hep-th].
  * [35] J. de Boer, J. Hartong, E. Have, N. A. Obers, and W. Sybesma, “Non-Boost Invariant Fluid Dynamics,” SciPost Phys. 9 no. 2, (2020) 018, arXiv:2004.10759 [hep-th].
  * [36] J. Armas and A. Jain, “Effective field theory for hydrodynamics without boosts,” SciPost Phys. 11 no. 3, (2021) 054, arXiv:2010.15782 [hep-th].
  * [37] M. Crossley, P. Glorioso, and H. Liu, “Effective field theory of dissipative fluids,” JHEP 09 (2017) 095, arXiv:1511.03646 [hep-th].
  * [38] P. Glorioso, M. Crossley, and H. Liu, “Effective field theory of dissipative fluids (II): classical limit, dynamical KMS symmetry and entropy current,” JHEP 09 (2017) 096, arXiv:1701.07817 [hep-th].
  * [39] N. Mullins, M. Hippert, and J. Noronha, “Stochastic fluctuations in relativistic fluids: Causality, stability, and the information current,” Phys. Rev. D 108 no. 7, (2023) 076013, arXiv:2306.08635 [nucl-th].
  * [40] L. Gavassino, N. Mullins, and M. Hippert, “Consistent inclusion of fluctuations in first-order causal and stable relativistic hydrodynamics,” arXiv:2402.06776 [nucl-th].
  * [41] M. Singh, C. Shen, S. McDonald, S. Jeon, and C. Gale, “Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions,” Nucl. Phys. A 982 (2019) 319–322, arXiv:1807.05451 [nucl-th].
  * [42] D. A. Teaney, Viscous Hydrodynamics and the Quark Gluon Plasma, pp. 207–266.  World Scientific, 2010.  arXiv:0905.2433 [nucl-th].
  * [43] R. E. Hoult and P. Kovtun, “Causality and classical dispersion relations,” Phys. Rev. D 109 no. 4, (2024) 046018, arXiv:2309.11703 [hep-th].
  * [44] A. Kurganov and E. Tadmor, “New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations,” J. Comput. Phys. 160 (2000) 241–282.
  * [45] L. D. Zanna and N. Bucciantini, “An efficient shock-capturing central-type scheme for multidimensional relativistic flows. I. hydrodynamics,” Astron. Astrophys. 390 (2002) 1177, arXiv:astro-ph/0205290.
  * [46] B. Schenke, S. Jeon, and C. Gale, “(3+1)D hydrodynamic simulation of relativistic heavy-ion collisions,” Phys. Rev. C 82 (2010) 014903, arXiv:1004.1408 [hep-ph].
  * [47] S. Gottlieb, D. I. Ketcheson, and C.-W. Shu, “High order strong stability preserving time discretizations,” Journal of Scientific Computing 38 no. 3, (Mar, 2009) 251–289. https://doi.org/10.1007/s10915-008-9239-z.
  * [48] S. Balay et al., “PETSc Web page.” https://petsc.org/, 2023.  https://petsc.org/.
  * [49] I. Dzyaloshinskii and G. Volovick, “Poisson brackets in condensed matter physics,” Annals of Physics 125 no. 1, (1980) 67–97. https://www.sciencedirect.com/science/article/pii/0003491680901190.
  * [50] D.-L. Wang and S. Pu, “Stability and causality criteria in linear mode analysis: Stability means causality,” Phys. Rev. D 109 no. 3, (2024) L031504, arXiv:2309.11708 [hep-th].