# Modified equilibrium distributions for Cooper–Frye particlization

M. McNelis  Department of Physics, The Ohio State University, Columbus, OH
43210, USA  U. Heinz  Department of Physics, The Ohio State University,
Columbus, OH 43210, USA

###### Abstract

We introduce a positive definite single-particle distribution that is suitable
for describing the transition from a macroscopic hydrodynamic to a microscopic
kinetic description during the late stages of heavy-ion collisions in the
presence of moderately large viscous corrections. The modified equilibrium
distribution function can be constructed with hydrodynamic input from either
relativistic viscous fluid dynamics or anisotropic fluid dynamics. We test the
modified equilibrium distribution’s hydrodynamic output for a stationary
hadron resonance gas subject to either shear stress, bulk pressure, or baryon
diffusion current at a given freeze-out temperature and baryon chemical
potential. While it does not reproduce all components of the net baryon
current and energy-momentum tensor exactly, it significantly improves upon the
customary linearized approximations for the non-equilibrium correction $\delta
f_{n}$ which typically lead to unphysical negative distribution functions at
large particle momenta. A comparison of particle spectra and
$p_{T}$–differential elliptic flow coefficients from the Cooper–Frye formula
computed with the modified equilibrium distribution and with linearized
$\delta f_{n}$ corrections is presented, for two different (2+1)–dimensional
hypersurfaces corresponding to central and non-central Pb+Pb collisions at the
Large Hadron Collider (LHC).

Cooper–Frye formula, ultrarelativistic heavy-ion collisions, relativistic
kinetic theory, 14–moment approximation, Chapman–Enskog expansion, modified
equilibrium distribution

###### pacs:

12.38.Mh, 25.75.-q, 24.10.Nz, 52.27.Ny, 51.10.+y

##  I Introduction

Among a large variety of computational approaches describing the evolution of
ultra-relativistic heavy-ion collisions, hybrid models have had the greatest
success in describing a wide variety of hadronic observables simultaneously
Song  _et al._ (2011); Schenke  _et al._ (2012); Gale  _et al._ (2013a, b);
Shen (7 25); Ryu  _et al._ (2015); Bernhard  _et al._ (2016); Bernhard (4 19).
Hybrid models include (at least) a relativistic viscous hydrodynamics module
describing the spacetime evolution of the early hot and dense quark-gluon
plasma (QGP) stage of the fireball created in the collision, followed by a
hadron cascade “afterburner” that follows the many-particle system of hadronic
resonances created from the QGP during the process of “hadronization”
microscopically to its final “kinetic freeze-out”. While experimental data
strongly suggest that the QGP stage evolves as a near-perfect liquid Heinz and
Kolb (2002); Romatschke and Romatschke (2007); Song  _et al._ (2011); Schenke
_et al._ (2011, 2010); Denicol  _et al._ (2012); Heinz and Snellings (2013),
its transport properties cannot be measured directly, nor can they (yet) be
reliably computed from first principles. Rather, they must be inferred using
various experimental probes, the most abundant of which (and therefore most
precisely measured) are the soft-momentum hadrons produced after the quark-
gluon plasma has cooled down below its pseudo-critical temperature.

Modeling the emission of these hadrons is therefore a key component in any
hybrid model Shen  _et al._ (2016). However, the exact spatial and momentum
configurations of these emitted hadrons are not well understood since their
distribution function is governed by microscopic kinetic theory while the
preceding quark-gluon plasma phase is treated as a strongly-coupled fluid. The
dynamical process in which the quark-gluon plasma converts to hadrons, known
as hadronization, is a highly complex, unresolved problem. Instead, hybrid
models bridge the transition between these two different phases of QCD matter
with a particlization model where the strongly-coupled quark-gluon plasma is
evolved hydrodynamically (i.e. macroscopically, with minimal microscopic input
such as the equation of state (EoS) which can be obtained from lattice QCD
Borsanyi  _et al._ (2010a, b); Bazavov  _et al._ (2012); Borsanyi  _et al._
(2014); Bazavov  _et al._ (2014) and parametrized transport coefficients)
until hadronization is complete, followed by an instant conversion to weakly-
interacting hadrons on a hypersurface $\Sigma(x)$ Huovinen and Petersen
(2012). The conversion between different sets of microscopic degrees of
freedom relies on the assumption that both viscous hydrodynamics and kinetic
theory are simultaneously valid on this surface $\Sigma(x)$.

Requiring that the energy, momentum and charges of the system are conserved
during this “particlization process”, the Cooper–Frye formula gives the
particle spectrum for hadron species $n$ as Cooper and Frye (1974)

$E_{n}\frac{dN_{n}}{d^{3}p}=\frac{1}{(2\pi\hbar)^{3}}\int_{\Sigma}p\cdot
d^{3}\sigma(x)\,f_{n}(x,p)\,,$ (1)

where $f_{n}$ is its distribution function in phase-space Cooper and Frye
(1974). For a locally equilibrated hadron resonance gas this distribution
function takes the form

$f_{\mathrm{eq},n}(x,p)=\frac{g_{n}}{\exp\Bigl{[}\dfrac{p\cdot
u(x)}{T(x)}-\alpha_{n}(x)\Bigr{]}+\Theta_{n}}\,,$ (2)

where $g_{n}$ is the spin degeneracy factor, $u^{\mu}(x)$ is the fluid
velocity, $T(x)$ is the temperature, $\alpha_{n}(x)=\mu_{n}(x)/T(x)$ is the
chemical potential-to-temperature ratio of species $n$, and
$\Theta_{n}\in[1,-1]$ accounts for the quantum statistics of fermions and
bosons, respectively. In this work, we only consider the baryon chemical
potential $\mu_{B}$ and write $\alpha_{n}=b_{n}\alpha_{B}$ where $b_{n}$ is
the baryon number of species $n$.

Due to dissipative effects in the preceding hydrodynamic evolution of the
fluid, the distribution function is generally out of local-equilibrium on the
conversion surface. One writes $f_{n}(x,p)=f_{\mathrm{eq},n}(x,p)+\delta
f_{n}(x,p)$ where $\delta f_{n}(x,p)$ encodes the deviation from local
equilibrium. The temperature and chemical potential for the equilibrium part
are obtained from the energy and net-baryon densities by Landau matching.
Deviations from local equilibrium are described in the hydrodynamic stage by
dissipative corrections to the energy-momentum tensor $T^{\mu\nu}(x)$ and the
net baryon current $J_{B}^{\mu}(x)$. These so-called “dissipative flows” are
the baryon diffusion current $V_{B}^{\mu}(x)$ (a spatial vector in the local
fluid rest frame (LRF), characterized by three degrees of freedom), the shear
stress tensor $\pi^{\mu\nu}(x)$ (a symmetric traceless rank-2 tensor with only
spatial components in the LRF, characterized by five degrees of freedom) and
the bulk viscous pressure $\Pi(x)$ (adding a single scalar degree of freedom).
In kinetic theory these dissipative flows are defined as momentum moments of
the off-equilibrium part $\delta f_{n}$ of the distribution functions, summed
over all hadron species $n$. At each spacetime point $x$ their values provide
altogether nine constraints on the momentum dependences of all of the off-
equilibrium corrections $\delta f_{n}(x,p)$, $n=1,\dots,N_{R}$, but obviously
these are not sufficient to fully determine them.

In this sense, the particlization problem is ill-posed. In contrast to a full
solution of the underlying kinetic equations, the hydrodynamic output provides
only limited information about $\delta f_{n}(x,p)$. To fully specify the form
of $f_{n}(x,p)$ in the Cooper–Frye formula (1), additional guiding principles
are needed to make optimal use of the hydrodynamic constraints. One well-known
scheme is Grad’s 14–moment approximation Grad (1949) in which the momentum
dependence of $\delta f_{n}$ is expanded up to second order in the momenta,
with coefficients matched to the 14 hydrodynamic moments of $\\{f_{n}\\}$,
$T^{\mu\nu}$ and $J^{\mu}_{B}$, and assuming that all particle species $n$
have the same expansion coefficients. Another is the first-order
Chapman–Enskog expansion of the Boltzmann equation in the relaxation time
approximation (RTA), with a common relaxation time for all hadron species
Chapman  _et al._ (1990); Anderson and Witting (1974). Both schemes treat
$\delta f_{n}$ as a linear perturbation of $f_{n}(x,p)$, ignoring higher-order
terms when matching to the hydrodynamic constraints. Unfortunately, these
linear corrections can turn the total distribution functions negative at high
momenta, especially for large dissipative flows. In particular, when
particlization is performed in close proximity to the quark-hadron phase
transition, dissipative corrections related to bulk viscosity (which peaks
near the pseudocritical temperature) tend to be large, and a linearized
approach can no longer be trusted Denicol  _et al._ (2009).

In the past decade, work has been done on constructing a non-equilibrium
distribution function suitable for Cooper–Frye particlization that does not
rely on a linearized expansion scheme. In the original work by Pratt and
Torrieri Pratt and Torrieri (2010), dissipative perturbations are added to the
“auxiliary fields” $T(x)$, $u^{\mu}(x)$ and $\alpha_{B}(x)$ in the Boltzmann
factor, effectively transforming the local-equilibrium distribution (2) into a
quasi-equilibrium distribution. The resulting modified equilibrium
distribution is given by the formula

$f_{\text{eq},n}^{\text{PT}}=\frac{\mathcal{Z}_{n}g_{n}}{\exp\biggl{[}\dfrac{\sqrt{{\bm{p}^{\prime
2}{+}m_{n}^{2}}}}{T{+}\delta
T}-b_{n}(\alpha_{B}{+}\delta\alpha_{B})\biggr{]}+\Theta_{n}}\,,$ (3)

where $\delta T=\lambda_{T}\Pi$ and $\delta\alpha_{B}=\lambda_{\alpha_{B}}\Pi$
are bulk viscous corrections to the effective temperature and chemical
potential, while $p^{\,\prime}_{i}=-X_{i}\cdot p^{\,\prime}$ are modified
local-rest-frame (LRF) momentum components, with
$X^{\mu}_{i}=(X^{\mu},Y^{\mu},Z^{\mu})$ being the spatial basis vectors (i.e.
the four-vectors that reduce in the LRF to the directional unit vectors along
the $x$, $y$, $z$ axes). The fluid velocity perturbation $\delta u^{\mu}$ is
encoded in the momentum space transformation

$p_{i}=A_{ij}p^{\,\prime}_{j}\,,$ (4)

where $p_{i}=-X_{i}\cdot p$ are the usual LRF momentum components and

$A_{ij}=(1{+}\lambda_{\Pi}\Pi)\delta_{ij}+\lambda_{\pi}\pi_{ij}$ (5)

is a symmetric matrix that deforms the momentum space linearly with the bulk
viscous pressure $\Pi(x)$ and LRF shear stress tensor
$\pi_{ij}(x)=X^{\mu}_{i}X^{\nu}_{j}\pi_{\mu\nu}(x)$.111The scalar coefficients
$\lambda_{T}$, $\lambda_{\alpha_{B}}$, $\lambda_{\Pi}$ and $\lambda_{\pi}$ are
functions of the temperature and chemical potential ($T$,
$\alpha_{B}$).,222Ref. Pratt and Torrieri (2010) did not consider a nonzero
baryon diffusion current $V_{B}^{\mu}(x)$. The normalization factor

$\mathcal{Z}_{n}=\frac{1}{{\det}A}$ (6)

rescales the distribution function such that its particle density agrees with
that of a local-equilibrium distribution with temperature $T{+}\delta T$ and
chemical potential $\alpha_{B}{+}\delta\alpha_{B}$.

The modified equilibrium distribution (3) is constructed such that its
momentum moments reproduce all components of $J_{B}^{\mu}$ and $T^{\mu\nu}$ at
first order in the (small) dissipative corrections. As the dissipative flows
$\Pi$ and $\pi_{ij}$ become larger, however, the mismatch between the
hydrodynamic input and the output when re-computing them as moments of the
modified distribution function (3) increases. This raises the question whether
the modifications to the local-equilibrium distribution can minimize this
mismatch even for moderately large viscous corrections. Tests conducted in
Ref. Pratt and Torrieri (2010) for the case of vanishing bulk viscous pressure
($\Pi=0$) showed that Eq. (3) accurately reproduces the target shear stress as
well as the energy and charge densities even for moderately large shear stress
modifications. However, when repeated for moderately large bulk viscous
pressures, this time setting the shear stress to $\pi_{ij}=0$, the modified
equilibrium distribution’s hydrodynamic output strongly deviated from the
input values. This suggests that the formula (3) can be further improved.

To reduce the errors associated with bulk viscous pressure, a variant of the
Pratt–Torrieri modified equilibrium distribution has been developed and
recently implemented in a Bayesian model parameter extraction of heavy-ion
collision simulations of Pb+Pb collisions at LHC energies
($\sqrt{s_{\text{NN}}}=2.76$ and 5.02 TeV) Bernhard (4 19). A different type
of quasi-equilibrium distribution, developed much earlier but with somewhat
similar goals in mind, is the anisotropic distribution Romatschke and
Strickland (2003) on which anisotropic fluid dynamics is based Martinez and
Strickland (2010). Although its analytic form has similarities with Eq. (3),
the relation between the modification parameters and the dissipative flows is
made nonlinear in order to capture the large pressure anisotropies present in
rapidly longitudinally expanding systems Heinz  _et al._ (2016). Anisotropic
fluid dynamical simulations have been interfaced with particlization models
based on the Romatschke–Strickland anisotropic distribution Romatschke and
Strickland (2003) to describe, with considerable success, the particle spectra
of Pb+Pb collisions at LHC energies ($\sqrt{s_{\text{NN}}}=2.76$ TeV) and
Au+Au collisions at RHIC energies ($\sqrt{s_{\text{NN}}}=200$ GeV) Alqahtani
_et al._ (2017); Almaalol  _et al._ (2019). A further generalization of this
approach has been proposed in Ref. Nopoush and Strickland (2019).

In this paper we present an improved formulation of the modified equilibrium
distribution (3). The perturbations made to the exponential argument of the
local-equilibrium distribution (2) are derived systematically from the
gradient expansion of the RTA Boltzmann equation. These modifications, which
are proportional to the hydrodynamic gradients, are constructed in such a way
that the modified equilibrium distribution reduces to the first-order RTA
Chapman–Enskog expansion in the limit of small viscous corrections. We test
its ability to reproduce the net baryon current and energy-momentum tensor of
a stationary hadron resonance gas subject to a variety of dissipative flows at
chemical freeze-out, including moderately large bulk viscous pressures and
baryon diffusion currents. We further apply the same method to also modify a
Romatschke–Strickland anisotropic distribution Romatschke and Strickland
(2003) at leading order (i.e. with only nonzero diagonal elements in the
momentum deformation matrix $A_{ij}$) in such a way that it can approximately
capture the smaller off-diagonal components of the energy-momentum tensor. We
demonstrate the potential use of our modified equilibrium (anisotropic)
distribution by computing the particle spectra from a longitudinally boost-
invariant Pb+Pb collision and comparing them to those computed with the
linearized $\delta f_{n}$ corrections as well as the modified equilibrium
distribution derived in Ref. Bernhard (4 19).

The paper is structured as follows: in Sec. II we review the linearized
viscous corrections to the hadronic distributions given by the 14–moment
approximation and RTA Chapman–Enskog expansion. In Sec. III we derive the
modified equilibrium distribution from the relaxation time approximation and
test the reproduction of $J_{B}^{\mu}$ and $T^{\mu\nu}$ for a stationary
hadron resonance gas. In Sec. IV, we repeat the procedure outlined in Sec. III
but starting with a leading order anisotropic distribution. In Sec. V we
compute the particle spectra from a (2+1)–dimensional viscous hydrodynamic
simulation of Pb+Pb collisions at LHC energies ($\sqrt{s_{\text{NN}}}=2.76$
TeV) using the Cooper–Frye formula and compare the differences between the
modified distributions and linearized $\delta f_{n}$ corrections.

##  II Linearized viscous corrections

###  II.1 14–moment approximation

The 14–moment approximation is a moments expansion of the distribution
function around $f_{\text{eq},n}$ that is truncated to the 14 lowest momentum
moments ($p^{\mu}$, $p^{\mu}p^{\nu}$) Grad (1949). This approximation assumes
that the distribution function can be adequately characterized by just its
hydrodynamic moments (i.e. $\int_{p}p^{\mu}f_{n}$,
$\int_{p}p^{\mu}p^{\nu}f_{n}$) Denicol  _et al._ (2012). For a multi-component
gas with nonzero baryon chemical potential, the 14–moment approximation reads
Monnai and Hirano (2009)

$\delta
f^{14}_{n}=f_{\text{eq},n}\bar{f}_{\text{eq},n}\left(b_{n}c_{\mu}p^{\mu}+c_{\mu\nu}p^{\mu}p^{\nu}\right)\,,$
(7)

where $\bar{f}_{\text{eq},n}=1-g_{n}^{-1}\Theta_{n}f_{\text{eq},n}$. For
simplicity, we follow common practice and take the expansion coefficients
$c_{\mu}$ and $c_{\mu\nu}$ to be species-independent. To solve for the
coefficients, one rewrites Eq. (7) in irreducible form,

$\begin{split}\delta
f^{14}_{n}&=f_{\text{eq},n}\bar{f}_{\text{eq},n}\Big{(}c_{T}m_{n}^{2}+b_{n}\big{(}c_{B}(u{\,\cdot\,}p){+}c_{V}^{\langle\mu\rangle}p_{\langle\mu\rangle}\big{)}\\\
&+c_{E}(u{\,\cdot\,}p)^{2}+c_{Q}^{\langle\mu\rangle}(u{\,\cdot\,}p)p_{\langle\mu\rangle}+c_{\pi}^{\langle{\mu\nu}\rangle}p_{\langle\mu}p_{\nu\rangle}\Big{)}\,,\end{split}$
(8)

where $p^{\langle\mu\rangle}=\Delta^{\mu}_{\nu}p^{\nu}$ and
$p^{\langle\mu}p^{\nu\rangle}=\Delta^{\mu\nu}_{\alpha\beta}p^{\alpha}p^{\beta}$,
with $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ and
$\Delta^{\mu\nu}_{\alpha\beta}=\frac{1}{2}\big{(}\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}{+}\Delta^{\mu}_{\beta}\Delta^{\nu}_{\alpha}\big{)}-\frac{1}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta}$,
are purely spatial in the LRF and traceless. The irreducible coefficients in
(8) are

$\displaystyle c_{B}$ $\displaystyle=u_{\mu}c^{\mu}\,,$ (9a) $\displaystyle
c_{V}^{\langle\mu\rangle}$ $\displaystyle=\Delta^{\mu}_{\nu}c^{\nu}\,,$ (9b)
$\displaystyle c_{T}$ $\displaystyle=g_{\mu\nu}c^{\mu\nu}\,,$ (9c)
$\displaystyle c_{E}$ $\displaystyle=u_{\mu}u_{\nu}c^{{\mu\nu}}\,,$ (9d)
$\displaystyle c_{Q}^{\langle\mu\rangle}$
$\displaystyle=u_{\alpha}\Delta^{\mu}_{\beta}c^{\alpha\beta}\,,$ (9e)
$\displaystyle c_{\pi}^{\langle{\mu\nu}\rangle}$
$\displaystyle=\Delta^{\mu\nu}_{\alpha\beta}\,c^{\alpha\beta}\,.$ (9f)

One solves for these coefficients as follows, by inserting Eq. (8) into the
Landau matching conditions for the net baryon density $n_{B}$ and energy
density $\mathcal{E}$, the definition of the Landau frame
$T^{\mu\nu}u_{\nu}=\mathcal{E}u^{\mu}$ (which implies a vanishing heat current
$Q^{\mu}=0$), and the kinetic definitions for the bulk viscous pressure $\Pi$,
baryon diffusion current $V_{B}^{\mu}$ and shear stress tensor $\pi^{\mu\nu}$:

$\displaystyle\delta n_{B}$
$\displaystyle=\sum_{n}\int_{p}b_{n}\,(u{\,\cdot\,}p)\delta f_{n}=0\,,$ (10a)
$\displaystyle\delta\mathcal{E}$
$\displaystyle=\sum_{n}\int_{p}(u{\,\cdot\,}p)^{2}\delta f_{n}=0\,,$ (10b)
$\displaystyle Q^{\mu}$
$\displaystyle=\sum_{n}\int_{p}(u{\,\cdot\,}p)p^{\langle\mu\rangle}\delta
f_{n}=0\,,$ (10c) $\displaystyle\Pi$
$\displaystyle=\frac{1}{3}\sum_{n}\int_{p}(-p{\,\cdot\,}\Delta{\,\cdot\,}p)\delta
f_{n}\,,$ (10d) $\displaystyle V_{B}^{\mu}$
$\displaystyle=\sum_{n}\int_{p}b_{n}\,p^{\langle\mu\rangle}\delta f_{n}\,,$
(10e) $\displaystyle\pi^{\mu\nu}$
$\displaystyle=\sum_{n}\int_{p}p^{\langle\mu}p^{\nu\rangle}\delta f_{n}\,.$
(10f)

Here the sum runs over the number of resonances, $n=1,\dots,N_{R}$, and
$\int_{p}\equiv\int\dfrac{d^{3}p}{(2\pi\hbar)^{3}E}$. After some algebra one
obtains

$\displaystyle c_{T}$
$\displaystyle=\frac{\Pi\,\mathcal{P}}{\mathcal{A}_{21}\mathcal{P}+\mathcal{N}_{31}\mathcal{Q}+\mathcal{J}_{41}\mathcal{R}}\,,$
(11a) $\displaystyle c_{B}$
$\displaystyle=\frac{\Pi\,\mathcal{Q}}{\mathcal{A}_{21}\mathcal{P}+\mathcal{N}_{31}\mathcal{Q}+\mathcal{J}_{41}\mathcal{R}}\,,$
(11b) $\displaystyle c_{E}$
$\displaystyle=\frac{\Pi\,\mathcal{R}}{\mathcal{A}_{21}\mathcal{P}+\mathcal{N}_{31}\mathcal{Q}+\mathcal{J}_{41}\mathcal{R}}\,,$
(11c) $\displaystyle c_{V}^{\langle\mu\rangle}$
$\displaystyle=\frac{V_{B}^{\mu}\mathcal{J}_{41}}{\mathcal{N}^{2}_{31}-\mathcal{M}_{21}\mathcal{J}_{41}}\,,$
(11d) $\displaystyle c_{Q}^{\langle\mu\rangle}$
$\displaystyle=-\frac{V_{B}^{\mu}\mathcal{N}_{31}}{\mathcal{N}^{2}_{31}-\mathcal{M}_{21}\mathcal{J}_{41}}\,,$
(11e) $\displaystyle c_{\pi}^{\langle{\mu\nu}\rangle}$
$\displaystyle=\frac{\pi^{\mu\nu}}{2(\mathcal{E}{+}\mathcal{P}_{\text{eq}})T^{2}}\,,$
(11f)

where $\mathcal{P}_{\text{eq}}(\mathcal{E},n_{B})$ is the equilibrium pressure
and

$\displaystyle\mathcal{P}=$
$\displaystyle\,\mathcal{N}^{2}_{30}-\mathcal{J}_{40}\mathcal{M}_{20}\,,$
(12a) $\displaystyle\mathcal{Q}=$
$\displaystyle\,\mathcal{B}_{10}\mathcal{J}_{40}-\mathcal{A}_{20}\mathcal{N}_{30}\,,$
(12b) $\displaystyle\mathcal{R}=$
$\displaystyle\,\mathcal{A}_{20}\mathcal{M}_{20}-\mathcal{B}_{10}\mathcal{N}_{30}\,.$
(12c)

The thermal integrals $\mathcal{J}_{kq}$, $\mathcal{N}_{kq}$,
$\mathcal{M}_{kq}$, $\mathcal{A}_{kq}$, and $\mathcal{B}_{kq}$ are defined in
Appendix A. It is well known that the coefficients of the 14–moment
approximation are linearly proportional to the dissipative flows of the fluid
De Groot  _et al._ (1980). An obvious problem with truncating the $\delta
f_{n}^{14}$ correction to first order is that it can overwhelm the local-
equilibrium distribution at sufficiently high momentum, potentially turning
the total distribution function negative. This problem grows worse with larger
dissipative flows. In principle, one can systematically improve the moments
expansion by including the non-hydrodynamic moments
($\int_{p}p^{\mu}p^{\nu}p^{\lambda}f_{n}$, etc.) Denicol  _et al._ (2012).
However, using them as macroscopic input along with $J_{B}^{\mu}$ and
$T^{\mu\nu}$ would require solving a set of evolution equations for these
higher-order moments together with the fluid dynamical simulation.

###  II.2 First-order RTA Chapman–Enskog expansion

Another common approach to obtaining a linear $\delta f_{n}$ correction for
the Cooper–Frye formula is the Chapman–Enskog expansion, which is a
perturbative series of the form

$f_{n}=f_{\text{eq},n}+\sum^{\infty}_{k=1}\delta
f_{n}^{(k)}=f_{\text{eq},n}+\sum^{\infty}_{k=1}\epsilon^{k}h_{n}^{(k)}\,,$
(13)

where $\epsilon$ and $h_{n}^{(k)}$ are the expansion parameter and
coefficients Jaiswal  _et al._ (2014). This series should satisfy the RTA
Boltzmann equation Anderson and Witting (1974)

$p^{\mu}\partial_{\mu}f_{n}=-\frac{(u{\,\cdot\,}p)(f_{n}{-}f_{\text{eq},n})}{\tau_{r}}\,,$
(14)

where we take the relaxation time $\tau_{r}(x)$ to be momentum and species
independent. If the hydrodynamic gradients are small compared to the
relaxation rate $\tau_{r}^{-1}$, the expansion parameter is the directional
derivative

$\epsilon=-s^{\mu}\partial_{\mu}$ (15)

with

$s^{\mu}(x)=\dfrac{\tau_{r}(x)}{u(x){\,\cdot\,}p}\,p^{\mu}.$ (16)

Inserting the ansatz (13) into Eq. (14) one finds at first order in gradients

$\delta f^{(1)}_{n}=-s^{\mu}\partial_{\mu}f_{\text{eq},n}\,.$ (17)

Using Eq. (2) and expanding the derivative one obtains

$\begin{split}\delta
f_{n}^{(1)}&=-\tau_{r}f_{\text{eq},n}\bar{f}_{\text{eq},n}\bigg{(}b_{n}\dot{\alpha}_{B}+\frac{(u{\,\cdot\,}p)\,\dot{T}}{T^{2}}+\frac{(-p{\,\cdot\,}\Delta{\,\cdot\,}p)\theta}{3(u{\,\cdot\,}p)T}\\\
&+\frac{b_{n}p^{\langle\mu\rangle}\nabla_{\mu}\alpha_{B}}{u{\,\cdot\,}p}-\frac{p^{\langle\mu\rangle}(\dot{u}_{\mu}{-}\nabla_{\mu}{\ln}T)}{T}-\frac{\sigma_{\mu\nu}p^{\langle\mu}p^{\nu\rangle}}{(u{\,\cdot\,}p)T}\bigg{)}\,,\end{split}$
(18)

where $\theta=\partial_{\mu}u^{\mu}$ is the scalar expansion rate and
$\sigma_{\mu\nu}=\partial_{\langle\mu}u_{\nu\rangle}\equiv\Delta^{\alpha\beta}_{\mu\nu}\partial_{\beta}u_{\alpha}$
is the velocity shear tensor. We denote the LRF time derivative as
$\dot{a}=u^{\mu}\partial_{\mu}a$ and the spatial gradient in the LRF as
$\nabla^{\mu}a=\Delta^{\mu\nu}\partial_{\nu}a$. To simplify Eq. (18) in terms
of the hydrodynamic quantities one makes use of the conservation equations for
the net baryon number, energy and momentum,

$\displaystyle\partial_{\mu}J_{B}^{\mu}$ $\displaystyle=0\,,$ (19a)
$\displaystyle u_{\nu}\partial_{\mu}T^{\mu\nu}$ $\displaystyle=0\,,$ (19b)
$\displaystyle\Delta^{\mu}_{\nu}\partial_{\lambda}T^{\lambda\nu}$
$\displaystyle=0\,,$ (19c)

to eliminate the time derivatives:

$\displaystyle\dot{\alpha}_{B}$ $\displaystyle\approx\mathcal{G}\theta\,,$
(20a) $\displaystyle\dot{T}$ $\displaystyle\approx\mathcal{F}\theta\,,$ (20b)
$\displaystyle\dot{u}^{\mu}$
$\displaystyle\approx\nabla^{\mu}\ln{T}+\frac{n_{B}T}{\mathcal{E}{+}\mathcal{P}_{\text{eq}}}\nabla^{\mu}\alpha_{B}\,.$
(20c)

Here we only keep the first-order terms; the coefficients $\mathcal{G}$ and
$\mathcal{F}$ are listed in Appendix B. The spatial gradients are substituted
by dissipative flows using the Navier–Stokes relations:

$\displaystyle\Pi$ $\displaystyle\approx-\zeta\theta\,,$ (21a) $\displaystyle
V_{B}^{\mu}$ $\displaystyle\approx\kappa_{B}\nabla^{\mu}\alpha_{B}\,,$ (21b)
$\displaystyle\pi^{\mu\nu}$ $\displaystyle\approx 2\eta\sigma^{\mu\nu}\,.$
(21c)

The first-order RTA Chapman–Enskog expansion (18) then reduces to

$\begin{split}\delta
f^{\text{CE}}_{n}&=f_{\text{eq},n}\bar{f}_{\text{eq},n}\bigg{[}\frac{\Pi}{\beta_{\Pi}}\bigg{(}b_{n}\mathcal{G}+\frac{(u{\,\cdot\,}p)\mathcal{F}}{T^{2}}+\frac{(-p\cdot\Delta\cdot
p)}{3(u{\,\cdot\,}p)T}\bigg{)}\\\
&+\frac{V_{B}^{\mu}p_{\langle\mu\rangle}}{\beta_{V}}\bigg{(}\frac{n_{B}}{\mathcal{E}{+}\mathcal{P}_{\text{eq}}}-\frac{b_{n}}{u{\,\cdot\,}p}\bigg{)}+\frac{\pi_{\mu\nu}p^{\langle\mu}p^{\nu\rangle}}{2\beta_{\pi}(u{\,\cdot\,}p)T}\bigg{]}\,,\end{split}$
(22)

where $\beta_{\Pi}$, $\beta_{V}$ and $\beta_{\pi}$ are the ratios of the bulk
viscosity, baryon diffusion coefficient and shear viscosity, respectively, to
the relaxation time. It is straightforward to check that $\delta
f^{\text{CE}}_{n}$ satisfies the Landau matching and frame conditions (10a-c).
The coefficients $\beta_{\Pi}$, $\beta_{V}$ and $\beta_{\pi}$ can be extracted
by inserting Eq. (22) into Eqs. (10d-f). The resulting expressions are

$\displaystyle\beta_{\Pi}$
$\displaystyle=\mathcal{G}n_{B}T+\frac{\mathcal{F}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})}{T}+\frac{5\mathcal{J}_{32}}{3T}\,,$
(23a) $\displaystyle\beta_{V}$
$\displaystyle=\mathcal{M}_{11}-\frac{n_{B}^{2}T}{\mathcal{E}{+}\mathcal{P}_{\text{eq}}}\,,$
(23b) $\displaystyle\beta_{\pi}$ $\displaystyle=\frac{\mathcal{J}_{32}}{T}\,.$
(23c)

Like the 14–moment approximation, the first-order Chapman–Enskog expansion
reproduces $J_{B}^{\mu}$ and $T^{\mu\nu}$ exactly, but the total distribution
function can turn negative at high momentum. In this case, one could improve
the expansion scheme by adding higher-order gradient corrections, but this
approach is tedious and will not be pursued here.

##  III Modified equilibrium distribution

###  III.1 Formulation

Both the 14–moment approximation (8) and the first-order RTA Chapman–Enskog
expansion (22) solve the problem of parametrizing the momentum dependence of
$\delta f_{n}(x,p)$ in terms of the available hydrodynamic information in such
a way that the dissipative flows are exactly reproduced from the corresponding
hydrodynamic moments of $\delta f_{n}$, but they also suffer from the negative
probability problem at sufficiently high momenta. The latter is a direct
consequence of truncating the expansion underlying each of the two approaches
at first order in the dissipative flows. One way to circumvent this problem is
to manipulate the argument of the exponential function in the local-
equilibrium distribution (2) with viscous corrections, as was done in Ref.
Pratt and Torrieri (2010). By construction, this avoids the negative
probability problem but, as we will see, at the expense of not being able to
match the dissipative flows exactly to the corresponding moments of this
modified equilibrium distribution.

Our starting point is the RTA Chapman–Enskog expansion. Equation (17) is the
first-order term of the Taylor series

$f_{\text{eq},n}(x{\,-\,}s,p)\approx
f_{\text{eq},n}(x,p)-s^{\mu}\partial_{\mu}f_{\text{eq},n}(x,p)\,,$ (24)

with $s^{\mu}$ given by Eq. (16). Therefore we have as a first approximation
for the modified equilibrium distribution

$f^{\text{(mod)}}_{\text{eq},n}(x,p)=f_{\text{eq},n}(x{\,-\,}s,p)\,,$ (25)

which retains the same analytic form as $f_{\text{eq},n}$, except for a shift
in the position arguments of the auxiliary fields $T(x)$, $\alpha_{B}(x)$ and
$u^{\mu}(x)$. Note that the shift depends on the particle momentum $p^{\mu}$
at which the distribution function is evaluated. Assuming small gradients we
can evaluate these corrections to first-order:

$\displaystyle T(x{\,-\,}s)$ $\displaystyle\approx
T(x)-s^{\nu}\partial_{\nu}T(x)\,,$ (26a) $\displaystyle\alpha_{B}(x{\,-\,}s)$
$\displaystyle\approx\alpha_{B}(x)-s^{\nu}\partial_{\nu}\alpha_{B}(x)\,,$
(26b) $\displaystyle u^{\mu}(x{\,-\,}s)$ $\displaystyle\approx
u^{\mu}(x)-s^{\nu}\partial_{\nu}u^{\mu}(x)\,.$ (26c)

The perturbations can be rewritten as

$\displaystyle\delta T$
$\displaystyle=-\tau_{r}\dot{T}-\frac{\tau_{r}p^{\langle\mu\rangle}\nabla_{\mu}T}{u{\,\cdot\,}p}\,,$
(27a) $\displaystyle\delta\alpha_{B}$
$\displaystyle=-\tau_{r}\dot{\alpha}_{B}-\frac{\tau_{r}p^{\langle\mu\rangle}\nabla_{\mu}\alpha_{B}}{u{\,\cdot\,}p}\,,$
(27b) $\displaystyle\\!\\!\\!\\!\delta u{\,\cdot\,}p$
$\displaystyle=\frac{\tau_{r}\theta(-p{\,\cdot\,}\Delta{\,\cdot\,}p)}{3(u{\,\cdot\,}p)}-\tau_{r}\dot{u}_{\mu}p^{\langle\mu\rangle}-\frac{\tau_{r}\sigma_{\mu\nu}p^{\langle\mu}p^{\nu\rangle}}{u{\,\cdot\,}p}\,.$
(27c)

As before, we can eliminate the time derivatives and spatial gradients by
using the conservation laws and Navier–Stokes relations. To obtain more
precise expressions for the perturbations, we first linearize the $\delta
f_{n}$ correction in Eq. (25) as

$\delta f_{n}\approx
f_{\text{eq},n}\bar{f}_{\text{eq},n}\bigg{(}b_{n}\delta\alpha_{B}+\frac{(u{\,\cdot\,}p)\delta
T}{T^{2}}-\frac{\delta u{\,\cdot\,}p}{T}\bigg{)}\,.$ (28)

After substituting Eqs. (27), (20) and (21) in Eq. (28) we recover the RTA
Chapman–Enskog expansion (22). We now compare Eqs. (22) and (28) to read off
the perturbations:333We here drop the spatial temperature gradients arising
from this procedure in Eqs. (27a,c) since they cancel in Eq. (28).

$\displaystyle\delta T$ $\displaystyle=\frac{\Pi\mathcal{F}}{\beta_{\Pi}}\,,$
(29a) $\displaystyle\delta\alpha_{B}$
$\displaystyle=\frac{\Pi\mathcal{G}}{\beta_{\Pi}}-\frac{V_{B}^{\mu}p_{\langle\mu\rangle}}{\beta_{V}(u{\,\cdot\,}p)}\,,$
(29b) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\delta u{\,\cdot\,}p$
$\displaystyle={-}\frac{\Pi({-}p{\cdot}\Delta{\cdot}p)}{3\beta_{\Pi}(u{\,\cdot\,}p)}-\frac{V_{B}^{\mu}p_{\langle\mu\rangle}n_{B}T}{\beta_{V}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})}-\frac{\pi_{\mu\nu}p^{\langle\mu}p^{\nu\rangle}}{2\beta_{\pi}(u{\,\cdot\,}p)}.$
(29c)

These corrections are now used to modify the exponent of Eq. (25) to
obtain444Equation (30) is closely related to the “maximum entropy” (ME)
distribution recently proposed in Ref. Everett  _et al._ (2021), which
modifies the exponent of the local-equilibrium distribution with additional
Lagrange multipliers. At first order in the dissipative flows, it was shown
Everett  _et al._ (2021) that the non-equilibrium corrections due to these
Lagrange multipliers reduce to the perturbative terms in Eq. (30).

$f_{\text{eq},n}^{\text{(mod)}}=\frac{g_{n}}{\exp\Bigl{[}\dfrac{(u{+}\delta
u){\,\cdot\,}p}{T{+}\delta
T}-b_{n}(\alpha_{B}{+}\delta\alpha_{B})\Bigr{]}+\Theta_{n}}\,.$ (30)

As long as the viscous corrections are not too large (say,
$|\Pi|<3\beta_{\Pi}$) this expression is positive definite, but it does not
reproduce the components of $J_{B}^{\mu}$ and $T^{\mu\nu}$ exactly since it
contains higher-order terms beyond the first-order RTA Chapman–Enskog
expansion. These discrepancies are at least second order in gradients and
should therefore be negligible if the viscous corrections are small:

$\displaystyle|\Pi|$ $\displaystyle\ll\beta_{\Pi}\,,$ (31a)
$\displaystyle\sqrt{V_{B,\mu}V_{B}^{\mu}}$ $\displaystyle\ll\beta_{V}\,,$
(31b) $\displaystyle\sqrt{\pi_{\mu\nu}\pi^{\mu\nu}}$ $\displaystyle\ll
2\beta_{\pi}.$ (31c)

However, Eq. (30) turns out to quickly lose its usefulness for even moderate
viscous corrections which are found to result in serious violations of the
matching to $J_{B}^{\mu}$ and $T^{\mu\nu}$. As it stands, Eq. (30) is not yet
sufficient and needs additional improvements in order to mitigate the errors
arising from the higher-order terms.555One option is to replace the
perturbations in Eq. (30) by the Lagrange multipliers in Ref. Everett  _et
al._ (2021) and adjust them so that the maximum entropy distribution exactly
reproduces the hydrodynamic moments. However, an exact (numerical) calculation
of these Lagrange multipliers for moderately large viscous corrections is
still outstanding.

###  III.2 Local momentum transformation and normalization

The main source of errors in the modified equilibrium distribution (30) comes
from the momentum dependent terms in the perturbations (29b,c). In the limit
of small gradients, the deformations and shifts to the momentum space are
linearly proportional to the dissipative flows. For larger viscous
corrections, however, these effects propagate into the hydrodynamic quantities
non-linearly. In order to control these errors, one should recast the
momentum-dependent perturbations in such a way that the viscous corrections to
the momentum space are strictly linear.

The local momentum transformation (4), first introduced in Ref. Pratt and
Torrieri (2010), provides an effective way of dealing with the momentum-
dependent perturbations in Eq. (29). Following their prescription, we rewrite
Eq. (30) as

$\\!\\!\\!\\!f_{\text{eq},n}^{\text{PTM}}=\frac{\mathcal{Z}_{n}g_{n}}{\exp\biggl{[}\dfrac{\sqrt{\bm{p}^{\prime
2}{+}m_{n}^{2}}}{T{+}\beta_{\Pi}^{-1}\Pi\mathcal{F}}-b_{n}\Big{(}\alpha_{B}{+}\dfrac{\Pi\mathcal{G}}{\beta_{\Pi}}\Big{)}\biggr{]}+\Theta_{n}}.\\!$
(32)

Starting from a thermal distribution of momenta $p^{\prime}$ we construct a
map $p=M(p^{\,\prime})$ to the LRF momenta $p$ such that the deformations and
shifts are linearly proportional to the dissipative flows. Under the
constraint of reducing to the momentum-dependent perturbations in the limit of
small gradients, one finds the following transformation between $p$ and
$p^{\prime}$:

$p_{i}=A_{ij}p^{\prime}_{j}-q_{i}\sqrt{\bm{p}^{\prime
2}{+}m_{n}^{2}}+b_{n}Ta_{i}\,,$ (33)

where

$\displaystyle A_{ij}$
$\displaystyle=\biggl{(}1{+}\frac{\Pi}{3\beta_{\Pi}}\biggr{)}\,\delta_{ij}+\frac{\pi_{ij}}{2\beta_{\pi}}\,,$
(34a) $\displaystyle q_{i}$
$\displaystyle=\frac{V_{B,i}\,n_{B}T}{\beta_{V}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})}\,,$
(34b) $\displaystyle a_{i}$ $\displaystyle=\frac{V_{B,i}}{\beta_{V}}\,,$ (34c)

with $V_{B,i}=-X_{i}{\,\cdot\,}V_{B}$ being the LRF spatial components of the
baryon diffusion current. The deformation matrix $A_{ij}$ is identical to the
one in Eqs. (4,5), i.e. $\lambda_{\Pi}=1/(3\beta_{\Pi})$ and
$\lambda_{\pi}=1/(2\beta_{\pi})$. Compared to Eqs. (4,5) Pratt and Torrieri
(2010), Eqs. (33,34) are further generalized to include baryon diffusion
effects.

For the normalization factor $\mathcal{Z}_{n}$ we fix the particle density of
the modified equilibrium distribution to that given by the RTA Chapman–Enskog
expansion,

$n^{(1)}_{n}=n_{\text{eq},n}+\frac{\Pi}{\beta_{\Pi}}\bigg{(}n_{\text{eq},n}+\mathcal{N}_{10,n}\mathcal{G}+\frac{\mathcal{J}_{20,n}\mathcal{F}}{T^{2}}\bigg{)}\,,$
(35)

which contains a bulk viscous correction to the equilibrium particle density
$n_{\text{eq},n}(T,\alpha_{B})$ Pratt and Torrieri (2010); Bernhard (4 19).
This suppresses the nonlinear shear and bulk viscous corrections to the
modified particle density,

$n^{\text{PTM}}_{n}=\mathcal{Z}_{n}\,{\det}A\,\cdot\,n_{\text{eq},n}(T{+}\beta_{\Pi}^{-1}\Pi\,\mathcal{F},\alpha_{B}{+}\beta_{\Pi}^{-1}\Pi\,\mathcal{G}),$
(36)

which (if $\mathcal{Z}_{n}=1$) are another major source of errors in matching
to the hydrodynamic output. Our choice for the normalization factor is
therefore

$\mathcal{Z}_{n}=\frac{1}{{\det}A}\times\frac{n_{n}^{(1)}}{n_{\text{eq},n}(T{+}\beta_{\Pi}^{-1}\Pi\,\mathcal{F},\alpha_{B}{+}\beta_{\Pi}^{-1}\Pi\,\mathcal{G})}\,,$
(37)

which is similar to Eq. (6) except it contains an additional factor. One
should keep in mind that the linearized density $n^{(1)}_{n}$ of light
hadrons, especially pions, may turn negative if the bulk viscous pressure is
too negative.

The “PTM distribution” (32) cannot be applied in every situation. The
requirements for using $f_{\text{eq},n}^{\text{PTM}}$ in the Cooper–Frye
formula are that the determinant of the Jacobian

$\det\bigg{(}\frac{\partial p_{i}}{\partial
p^{\prime}_{j}}\bigg{)}={\det}A\,\bigg{(}1{\,-\,}\frac{q_{i}A^{-1}_{ij}p^{\prime}_{j}}{\sqrt{\bm{p}^{\prime
2}{+}m_{n}^{2}}}\bigg{)}$ (38)

is positive for any value of $p^{\prime}$, the deformation matrix $A_{ij}$ is
invertible, and the normalization factor $\mathcal{Z}_{n}$ is non-negative.
These conditions can be violated when the dissipative flows are too large.
What to do in this case will be discussed in Sec. V where we study certain
situations where the modified equilibrium distribution breaks down.

###  III.3 Reproducing the hydrodynamic quantities

As already noted, the nonlinear dependence of the modified equilibrium
distribution (32) on the dissipative flows makes it impossible to achieve an
exact matching of the hydrodynamic moments of this distribution with all
components of $T^{\mu\nu}$ and $J_{B}^{\mu}$. This is different from the
linearized parametrizations described in Sec. II. In this subsection, we study
the matching violations for hadrons emitted from a stationary time-like
freeze-out cell with the modified equilibrium distribution (32).

The system is composed of the set of hadron species included in the hadronic
afterburner code URQMD Bleicher  _et al._ (1999) ($N_{R}\sim 320$). These
hadrons are assumed to be produced from a fluid in chemical equilibrium at the
chemical freeze-out (CF) temperature $T_{\text{CF}}$ and baryon chemical
potential $\mu_{B,\text{CF}}$, which can vary across different collision
systems. The values used here (shown in the legend of Figure 1) are taken from
a statistical model fit to the hadron abundance ratios measured in Pb+Pb
collisions at LHC energies ($\sqrt{s_{\mathrm{NN}}}=5.02$ TeV) and in Au+Au
collisions at RHIC and SPS energies ($\sqrt{s_{\text{NN}}}=200$ GeV and 17.3
GeV, respectively) Andronic  _et al._ (2018). For clarity, we vary either the
input shear stress, bulk viscous pressure or baryon diffusion current while
fixing the other dissipative flows to zero.

In the first example (left column of Fig. 1) we consider the case where in a
given freeze-out cell the hydrodynamic $T^{\mu\nu}$ features a negative
pressure anisotropy
$\frac{2}{3}(\mathcal{P}_{L}{-}\mathcal{P}_{\perp})=\pi^{zz}=-2\pi^{xx}=-2\pi^{yy}$,
combined with zero bulk viscous pressure and baryon diffusion current.
Starting from equilibrium at the right edge of the plots, we decrease
$\pi^{zz}$ going left, causing the longitudinal pressure
$\mathcal{P}_{L}=\mathcal{P}_{\text{eq}}+\pi^{zz}$ to decrease and the
transverse pressure
$\mathcal{P}_{\perp}=\mathcal{P}_{\text{eq}}+\frac{1}{2}(\pi^{xx}{+}\pi^{yy})$
to increase. The energy density should not change due to the Landau matching
condition (10a). These are the hydrodynamic components of $T^{\mu\nu}$ that we
aim to reproduce, shown by the solid black curves. The colored dashed curves
are the kinetic outputs, calculated as moments of the modified equilibrium
distribution with these hydrodynamic inputs.666Due to the normalization factor
(37) the modified equilibrium distribution conserves the net baryon number
exactly ($\delta n_{B}{\,=\,}0$).,777For only shear stress inputs, the outputs
$Q^{\mu}$ and $V_{B}^{\mu}$ vanish by symmetry and are therefore not of
interest here. One sees that for small $\pi^{zz}$ the kinetic output closely
follows the initial hydrodynamic input. This is because in the limit of small
dissipative flows the modified equilibrium distribution reduces to the linear
RTA Chapman–Enskog expansion for which the kinetic output reproduces the
hydrodynamic input exactly. As $\pi^{zz}$ further decreases to larger negative
values, the kinetic output $T^{\mu\nu}$ begins to deviate from the
hydrodynamic target. In particular, for a positive definite distribution
function such as Eq. (32) the kinetic output $\mathcal{P}_{L}$ stays always
above zero even if the hydrodynamic input for $\pi^{zz}$ is so large and
negative that the total longitudinal pressure in the fluid is negative. On the
other hand, the kinetic outputs for $\mathcal{E}$ and $\mathcal{P}_{\perp}$
agree well with their hydrodynamic targets even for large pressure
anisotropies. These trends are confirmed for all three combinations
$(T_{\text{CF}},\mu_{B,\text{CF}})$ studied in the three rows of Fig. 1.
Technically, the modified equilibrium distribution breaks down completely for
$\pi^{zz}\leq-2\beta_{\pi}$ (which causes ${\det}A\leq 0$), and one should not
expect $f_{\text{eq},n}^{\text{PTM}}$ to work well for very large pressure
anisotropies. However, for moderately large
$|\pi^{zz}|\leq\frac{1}{2}\mathcal{P}_{\text{eq}}$,
$f_{\text{eq},n}^{\text{PTM}}$ reproduces the target $T^{\mu\nu}$ components
quite well, with $|\Delta\mathcal{P}_{L}|/\mathcal{P}_{L}\leq 9.0\%$,
$|\Delta\mathcal{P}_{\perp}|/\mathcal{P}_{\perp}\leq 0.5\%$ and
$|\Delta\mathcal{E}|/\mathcal{E}\leq 0.8\%$.888We also successfully replicated
the shear stress test shown in Fig. 2 of Ref. Pratt and Torrieri (2010).

Figure 1:  The reproduction of selected components of the net baryon current
$J_{B}^{\mu}$ and energy-momentum tensor $T^{\mu\nu}$ using the Cooper–Frye
prescription with the PTM distribution (32), for a stationary hadron resonance
gas subject to either shear stress (left column), bulk viscous pressure
(middle column) or a baryon diffusion current (right column). The temperature
and baryon chemical potential $(T_{\text{CF}},\mu_{B,\text{CF}})$ at chemical
freeze-out are indicated in the legend (top to bottom rows). Solid black
curves represent the hydrodynamic input (i.e. the target $J_{B}^{\mu}$ and
$T^{\mu\nu}$) while the dashed color lines are the output components of
$J_{B}^{\mu}$ and $T^{\mu\nu}$ obtained by computing the corresponding moments
of $f_{\text{eq}}^{\text{PTM}}$.

In the middle column of Fig. 1 we vary the bulk viscous pressure $\Pi$ at
vanishing shear stress and baryon diffusion. We checked that for this
hydrodynamic input the kinetic outputs for $\pi^{\mu\nu}$, $Q^{\mu}$ and
$V_{B}^{\mu}$ vanish. The plots show that, for all three choices of freeze-out
parameters, the energy matching condition $\Delta\mathcal{E}=0$ holds very
well even for large (negative) values of $\Pi$. The kinetic output for the
total pressure $\mathcal{P}_{\text{eq}}+\Pi$ tends to somewhat underpredict
the hydrodynamic target, the error staying below 10% for moderately large bulk
viscous pressures $|\Pi|\leq\frac{1}{2}\mathcal{P}_{\text{eq}}$. For more
negative input values of $\Pi$, $\Pi\lesssim-2\beta_{\Pi}$, the particle
density for pions turns negative and the modified equilibrium distribution
becomes invalid.

In the right column of Fig. 1 we test the reproduction of the baryon diffusion
current (taken arbitrarily to point in the $z$–direction) at zero shear and
bulk viscous stress. One sees that the hydrodynamic input for the ratio
$\alpha_{B}V_{B}^{z}/n_{B}$ (which is approximately the same for all three
freeze-out parameter pairs) is well reproduced by the kinetic output from the
modified equilibrium distribution. The energy matching and Landau frame
($Q^{z}=0$) conditions are also reproduced with excellent precision. However,
a nonzero input for $V_{B}^{z}\neq 0$ also generates nonzero kinetic outputs
for the bulk and shear viscous stresses even when their hydrodynamic inputs
are zero: For $V_{B}^{z}/\beta_{V}=2$ we find a 5% positive difference
$\Delta\mathcal{P}_{L}/\mathcal{P}_{\text{eq}}$ between the hydrodynamic input
and the kinetic output for the longitudinal pressure of which, after
decomposition into shear and bulk viscous contributions, 1.5% can be
attributed to an induced bulk viscous pressure $\Pi$ and 3.5% to an induced
shear stress $\pi^{zz}$, both with positive signs. For moderately large baryon
diffusion currents $|V_{B}^{z}|\leq\beta_{V}$, the hydrodynamic output errors
are $|\Delta V_{B}^{z}/V_{B}^{z}|\leq 0.5\%$,
$|\Delta\mathcal{P}_{L}|/\mathcal{P}_{\text{eq}}\leq 1.3\%$,
$|\Delta\mathcal{E}|/\mathcal{E}\leq 0.2\%$ and $|Q^{z}|/\mathcal{E}\approx
0$.

Figure 2:  The isotropic scale parameter $\lambda_{\Pi}$ (solid red) and
normalization factor $z_{\Pi}$ (solid blue) as a function of
$\Pi/\mathcal{P}_{\text{eq}}$ are computed for a hadron resonance gas at a
fixed temperature $T=150$ MeV. The dotted blue and red lines are the upper
bounds for $z_{\Pi}$ and $\Pi/\mathcal{P}_{\text{eq}}$, respectively. This
plot corresponds to Fig. 3.12 in Ref. Bernhard (4 19), where $\Delta\langle
p\rangle/\langle p\rangle_{0}=\lambda_{\Pi}$ and $\Delta n/n_{0}=z_{\Pi}-1$.

###  III.4 Pratt–Torrieri–Bernhard distribution

Another variant of Pratt and Torrieri’s idea Pratt and Torrieri (2010) was
implemented by Bernhard in Ref. Bernhard (4 19):

$f_{\text{eq},n}^{\text{PTB}}=\frac{\mathcal{Z}g_{n}}{\exp\Bigl{[}\dfrac{\sqrt{\bm{p}^{\prime
2}{+}m_{n}^{2}}}{T}\Bigr{]}+\Theta_{n}}\,.$ (39)

In this “PTB distribution” the baryon chemical potential and diffusion current
are neglected ($\alpha_{B}=V_{B}^{\mu}=0$), and the effective temperature is
not modified ($\delta T=0$). The momentum transformation rule is
$p_{i}=A_{ij}p^{\prime}_{j}$ with

$A_{ij}=(1{+}\lambda_{\Pi})\delta_{ij}+\frac{\pi_{ij}}{2\beta_{\pi}}\,.$ (40)

The shear stress modification is the same as in Eq. (34a) while the bulk
pressure term is replaced by the isotropic scale parameter $\lambda_{\Pi}$.
The normalization factor $\mathcal{Z}$ is taken as species-independent:

$\mathcal{Z}=\frac{z_{\Pi}}{{\det}A}\,,$ (41)

leaving particle abundance ratios unchanged from their chemical equilibrium
values. The parameters $\lambda_{\Pi}$ and $z_{\Pi}$ are fixed such that in
the absence of shear stress ($\pi_{ij}=0$) the energy density and bulk viscous
pressure are exactly matched Bernhard (4 19):

$\displaystyle\mathcal{E}^{\prime}(\lambda_{\Pi},z_{\Pi},T)$
$\displaystyle=\mathcal{E}\,,$ (42a)
$\displaystyle\mathcal{P}^{\prime}(\lambda_{\Pi},z_{\Pi},T)$
$\displaystyle=\mathcal{P}_{\text{eq}}+\Pi\,;$ (42b)

here $\mathcal{E}^{\prime}$ and $\mathcal{P}^{\prime}$ are the kinetic theory
output for the energy density and isotropic pressure computed from the PTB
distribution (39).

The parametrization procedure is as follows: For the distribution function
(39) to be well-defined the isotropic scale parameter must satisfy
$\lambda_{\Pi}\in(-1,\infty)$. At a given temperature $T$ one sets up a grid
in $\lambda_{\Pi}$ and computes the corresponding values for $z_{\Pi}$ and
$\Pi$ numerically by rewriting Eqs. (42a,b) as

$\displaystyle z_{\Pi}$
$\displaystyle=\frac{\mathcal{E}}{\mathcal{L}_{20}(\lambda_{\Pi},T)}\,,$ (43a)
$\displaystyle\Pi$
$\displaystyle=z_{\Pi}\mathcal{L}_{21}(\lambda_{\Pi},T)-\mathcal{P}_{\text{eq}}\,.$
(43b)

The functions $\mathcal{L}_{kq}$ are defined in Appendix A. One can then
interpolate the data with respect to $\Pi$ to construct the functions
$\lambda_{\Pi}(\Pi;T)$ and $z_{\Pi}(\Pi;T)$; an example is shown in Figure 2.
Generally, these are nonlinear functions of $\Pi$. In the limit of small bulk
viscous pressure, they linearize to

$\displaystyle z_{\Pi}$ $\displaystyle\approx 1+\delta
z_{\Pi}=1-\frac{3\Pi\mathcal{P}_{\text{eq}}}{5\beta_{\pi}\mathcal{E}-3\mathcal{P}_{\text{eq}}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})}\,,$
(44a) $\displaystyle\lambda_{\Pi}$
$\displaystyle\approx\delta\lambda_{\Pi}=\frac{\Pi\mathcal{E}}{5\beta_{\pi}\mathcal{E}-3\mathcal{P}_{\text{eq}}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})}\,.$
(44b)

Equations (43a,b) imply that $z_{\Pi}$ and $\Pi$ are bounded by

$\displaystyle 0<$ $\displaystyle\,z_{\Pi}<\frac{\mathcal{E}}{\rho}\,,$ (45a)
$\displaystyle-\mathcal{P}_{\text{eq}}<$
$\displaystyle\,\Pi<\frac{\mathcal{E}}{3}-\mathcal{P}_{\text{eq}}\,,$ (45b)

where $\rho=\sum_{n}^{N_{R}}m_{n}n_{\text{eq},n}$ is the equilibrium mass
density. When $\Pi$ lies outside the bound999Violations of the upper bound in
Eq. (45b) can occur at the transition from a conformal pre-hydrodynamic model
to non-conformal viscous hydrodynamics, where the mismatch between the
conformal and QCD equations of state gives exactly
$\Pi=\mathcal{E}/3-\mathcal{P}_{\text{eq}}(\mathcal{E})$. (45b) or
${\det}A\leq 0$, we consider the modified equilibrium distribution (39) to
have broken down. As will be discussed in Sec. V.2, we then resort to
linearizing Eq. (39) around local equilibrium, i.e. we write
$f_{\text{eq},n}^{\text{PTB}}\approx f_{\text{eq},n}+\delta f_{n}$, where

$\begin{split}\delta f_{n}&=f_{\text{eq},n}\left(\delta
z_{\Pi}-3\delta\lambda_{\Pi}\right)+\\\
&\,f_{\text{eq},n}\bar{f}_{\text{eq},n}\left(\frac{\delta\lambda_{\Pi}(-p\cdot\Delta\cdot
p)}{(u\cdot
p)T}+\frac{\pi_{\mu\nu}p^{\langle\mu}p^{\nu\rangle}}{2\beta_{\pi}(u\cdot
p)T}\right)\,.\end{split}$ (46)

The effectiveness of the PTB distribution (39) in reproducing the components
of a given energy-momentum tensor $T^{\mu\nu}$ was studied in Ref. Bernhard (4
19) to which we refer the reader for comparison with Fig. 1.

##  IV Modified anisotropic distribution

###  IV.1 Anisotropic hydrodynamics

In this section we apply the formulation developed in Sec. II to modify an
anisotropic distribution function that has already been deformed at leading
order. In anisotropic hydrodynamics, the distribution function is decomposed
into a momentum-anisotropic leading order term and a residual correction Bazow
_et al._ (2014); Molnar  _et al._ (2016); McNelis  _et al._ (2018) as follows:

$f_{n}(x,p)=f_{a,n}(x,p)+\delta\tilde{f}_{n}(x,p)\,.$ (47)

Here the anisotropic leading order distribution $f_{a,n}$ is taken as101010For
simplicity, we set in this Section the effective baryon chemical potential to
$\tilde{\mu}_{B}=0$.

$f_{a,n}(x,p)=\frac{g_{n}}{\exp\biggl{[}\dfrac{\sqrt{p\cdot\Omega(x)\cdot
p+m_{n}^{2}}}{\Lambda(x)}\biggr{]}+\Theta_{n}}\,.$ (48)

Similar to the modified equilibrium distribution (32), it has an effective
temperature $\Lambda$ and a momentum transformation encoded in the ellipsoidal
tensor

$\Omega_{\mu\nu}=\frac{Z_{\mu}Z_{\nu}}{\alpha_{L}^{2}}-\frac{\Xi_{\mu\nu}}{\alpha_{\perp}^{2}}\,.$
(49)

Here $Z^{\mu}{\,\equiv\,}X_{3}^{\mu}$ is the longitudinal basis vector and
$\Xi^{\mu\nu}=g^{\mu\nu}{-}u^{\mu}u^{\nu}{+}Z^{\mu}Z^{\nu}$ the associated
transverse spatial projector Molnar  _et al._ (2016). The momentum anisotropy
parameters $\alpha_{L}$ and $\alpha_{\perp}$ deform the longitudinal and
transverse momentum space, respectively. The anisotropic distribution (48) can
be rewritten more conveniently as

$f_{a,n}=\frac{g_{n}}{\exp\Bigl{[}\dfrac{\sqrt{\bm{p}^{\prime
2}{+}m_{n}^{2}}}{\Lambda}\Bigr{]}+\Theta_{n}}\,,$ (50)

where $p_{i}=A_{ij}p^{\prime}_{j}$ with

$A_{ij}=\alpha_{\perp}\delta_{ij}+(\alpha_{L}{-}\alpha_{\perp})Z_{i}Z_{j}\,,$
(51)

and $Z_{i}=-X_{i}{\,\cdot\,}Z=(0,0,1)$ being the LRF components of the
longitudinal basis vector. Together, the anisotropic parameters $\Lambda$,
$\alpha_{L}$ and $\alpha_{\perp}$ are adjusted such that $f_{a,n}$ completely
captures the energy density $\mathcal{E}$ and the longitudinal and transverse
pressures, $\mathcal{P}_{L}$ and $\mathcal{P}_{\perp}$, in the energy-momentum
tensor Molnar  _et al._ (2016); McNelis  _et al._ (2018)

$\\!\\!\\!\\!T^{\mu\nu}=\mathcal{E}u^{\mu}u^{\nu}{+}\mathcal{P}_{L}Z^{\mu}Z^{\nu}{-}\mathcal{P}_{\perp}\Xi^{\mu\nu}+2W_{\perp
z}^{(\mu}Z^{\nu)}+\pi_{\perp}^{\mu\nu}.\\!$ (52)

The residual correction $\delta\tilde{f}_{n}$ does not contribute to
$\mathcal{E}$, $\mathcal{P}_{L}$ and $\mathcal{P}_{\perp}$ but accounts for
the longitudinal momentum diffusion current $W_{\perp z}^{\mu}$ and the
transverse shear stress tensor $\pi_{\perp}^{\mu\nu}$:

$\displaystyle W_{\perp z}^{\mu}$
$\displaystyle=-\Xi^{\mu}_{\alpha}Z_{\nu}T^{\alpha\nu}\,,$ (53a)
$\displaystyle\pi_{\perp}^{\mu\nu}$
$\displaystyle=\Xi^{\mu\nu}_{\alpha\beta}T^{\alpha\beta}\,.$ (53b)

Here
$\Xi^{\mu\nu}_{\alpha\beta}=\frac{1}{2}\big{(}\Xi^{\mu}_{\alpha}\Xi^{\nu}_{\beta}+\Xi^{\nu}_{\beta}\Xi^{\mu}_{\alpha}-\Xi^{\mu\nu}\Xi_{\alpha\beta}\big{)}$
is the double-transverse traceless projector. An expression for
$\delta\tilde{f}_{n}$ can be derived using a linearization approach, such as
the Chapman–Enskog expansion. After substituting Eq. (47) in the RTA Boltzmann
equation (14), the first-order expression for $\delta\tilde{f}_{n}$ is

$\delta\tilde{f}_{n}=-s^{\mu}\partial_{\mu}f_{a,n}-(f_{a,n}{-}f_{\text{eq},n})\,,$
(54)

where the second-order term $-s^{\mu}\partial_{\mu}\delta\tilde{f}_{n}$ was
neglected. The first term in Eq. (54) can be expanded into

$-s^{\mu}\partial_{\mu}f_{a,n}=f_{a,n}\bar{f}_{a,n}\left(\dfrac{E_{a,n}\delta\Lambda}{\Lambda^{2}}-\frac{p\cdot\delta\Omega\cdot
p}{2E_{a,n}\Lambda}\right)\,,$ (55)

where $\bar{f}_{a,n}{\,=\,}1{\,-\,}g_{n}^{-1}\Theta_{n}f_{a,n}$,
$E_{a,n}{\,=\,}\sqrt{p{\,\cdot\,}\Omega{\,\cdot\,}p{+}m_{n}^{2}}$,

$\displaystyle\frac{p\cdot\delta\Omega\cdot p}{2}$ $\displaystyle=$
$\displaystyle\frac{(\alpha_{L}^{2}{-}\alpha_{\perp}^{2})(-Z\cdot p)(\delta
Z\cdot p)}{\alpha_{\perp}^{2}\alpha_{L}^{2}}+\frac{(u{\,\cdot\,}p)(\delta
u{\,\cdot\,}p)}{\alpha_{\perp}^{2}}$ (56)
$\displaystyle-\frac{\delta\alpha_{L}(-Z\cdot
p)^{2}}{\alpha_{L}^{3}}-\frac{\delta\alpha_{\perp}(-p\cdot\Xi\cdot
p)}{\alpha_{\perp}^{3}}\,,$

and the perturbations $\delta u^{\mu}$ and $\delta Z^{\mu}$ are111111We find
that the $\delta\tilde{f}_{n}$ terms $\propto$ ($\delta\Lambda$,
$\delta\alpha_{L}$, $\delta\alpha_{\perp}$) do not contribute to $W_{\perp
z}^{\mu}$ and $\pi_{\perp}^{\mu\nu}$ so we can effectively set these
perturbations to zero.

$\displaystyle\delta u{\,\cdot\,}p$
$\displaystyle=\frac{-\tau_{r}}{u{\,\cdot\,}p}\Big{[}(u\cdot p)\big{(}(-Z\cdot
p)Z^{\mu}{+}p^{\\{\mu\\}}\big{)}\dot{u}_{\mu}$ (57a) $\displaystyle-(-Z\cdot
p)^{2}\theta_{L}-\frac{1}{2}(-p\cdot\Xi\cdot p)\theta_{\perp}$
$\displaystyle+p^{\\{\mu}p^{\nu\\}}\sigma_{\perp,{\mu\nu}}+(-Z\cdot
p)p_{\\{\mu\\}}\big{(}Z_{\nu}\nabla_{\perp}^{\mu}u^{\nu}{-}D_{z}u^{\mu}\big{)}\Big{]}\,,$
$\displaystyle\delta Z\cdot p$
$\displaystyle=\frac{-\tau_{r}}{u{\,\cdot\,}p}\Big{[}{-}(u\cdot
p)^{2}Z_{\mu}\dot{u}^{\mu}+(u{\,\cdot\,}p)(-Z\cdot p)\theta_{L}$ (57b)
$\displaystyle+(u{\,\cdot\,}p)(-Z\cdot
p)\theta_{L}+p^{\\{\mu}p^{\nu\\}}\sigma_{z,{\mu\nu}}$ $\displaystyle-(-Z\cdot
p)p_{\\{\mu\\}}D_{z}Z^{\mu}-\frac{1}{2}(-p\cdot\Xi\cdot
p)\nabla_{\perp\mu}Z^{\mu}$ $\displaystyle+(u\cdot
p)p_{\\{\mu\\}}\big{(}\dot{Z}^{\mu}{-}Z_{\nu}\nabla_{\perp}^{\mu}u^{\nu}\big{)}\Big{]}\,,$

where $p^{\\{\mu\\}}\equiv\Xi^{\mu}_{\nu}\,p^{\nu}$,
$p^{\\{\mu}p^{\nu\\}}\equiv\Xi^{\mu\nu}_{\alpha\beta}\,p^{\alpha}p^{\beta}$,
$\theta_{L}=Z_{\mu}D_{z}u^{\mu}$ is the longitudinal expansion rate,
$\theta_{\perp}=\nabla_{\perp,\mu}u^{\mu}$ is the transverse expansion rate,
$D_{z}=-Z{\,\cdot\,}\partial$ is the LRF longitudinal derivative,
$\nabla_{\perp}^{\mu}=\Xi^{\mu\nu}\partial_{\nu}$ is the transverse spatial
gradient,
$\sigma_{\perp,{\mu\nu}}=\Xi^{\alpha\beta}_{\mu\nu}\partial_{\alpha}u_{\beta}$
is the transverse velocity shear tensor and
$\sigma_{z,{\mu\nu}}=\Xi^{\alpha\beta}_{\mu\nu}\partial_{\alpha}Z_{\beta}$
Molnar  _et al._ (2016). For the sake of simplicity, we only consider the
terms with nonzero contributions to $W_{\perp z}^{\mu}$ and
$\pi_{\perp}^{\mu\nu}$, whose kinetic theory definitions are

$\displaystyle W_{\perp z}^{\mu}$ $\displaystyle=\sum_{n}\int_{p}(-Z\cdot
p)p^{\\{\mu\\}}\delta\tilde{f}_{n}\,,$ (58a)
$\displaystyle\pi_{\perp}^{\mu\nu}$
$\displaystyle=\sum_{n}\int_{p}p^{\\{\mu}p^{\nu\\}}\delta\tilde{f}_{n}\,.$
(58b)

By symmetry, the term $(f_{a,n}{-}f_{\text{eq},n})$ in Eq. (54) and
corrections $\propto(\delta\Lambda,\delta\alpha_{L},\delta\alpha_{\perp})$
have zero contributions, so we can effectively eliminate them. The only
nonzero contributions from $\delta u{\,\cdot\,}p$ and $\delta Z{\,\cdot\,}p$
in Eqs. (57a,b) are

$\displaystyle\delta u{\,\cdot\,}p\rightarrow$
$\displaystyle-\frac{\tau_{r}(-Z\cdot
p)p_{\\{\mu\\}}\left(Z_{\nu}\nabla_{\perp}^{\mu}u^{\nu}{-}D_{z}u^{\mu}\right)}{u{\,\cdot\,}p}$
(59a)
$\displaystyle-\frac{\tau_{r}p_{\\{\mu}p_{\nu\\}}\sigma_{\perp}^{\mu\nu}}{u{\,\cdot\,}p}\,,$
$\displaystyle\delta Z\cdot p\rightarrow$
$\displaystyle-\tau_{r}p_{\\{\mu\\}}\big{(}\dot{Z}^{\mu}{-}Z_{\nu}\nabla_{\perp}^{\mu}u^{\nu}\big{)}\,.$
(59b)

The $\delta\tilde{f}_{n}$ correction then reduces to

$\begin{split}\delta\tilde{f}_{n}&=\frac{\tau_{r}f_{a,n}\bar{f}_{a,n}}{E_{a,n}\Lambda}\left[\frac{p_{\\{\mu}p_{\nu\\}}\sigma_{\perp}^{\mu\nu}}{\alpha_{\perp}^{2}}{-}\frac{({-}Z{\,\cdot\,}P)p_{\\{\mu\\}}\kappa_{\perp}^{\mu}}{\alpha_{\perp}\alpha_{L}}\right]\,,\end{split}$
(60)

where

$\kappa_{\perp}^{\mu}=\frac{(\alpha_{\perp}^{2}{-}\alpha_{L}^{2})\Xi^{\mu}_{\nu}\dot{Z}^{\nu}{-}\alpha_{\perp}^{2}Z_{\nu}\nabla_{\perp}^{\mu}u^{\nu}{+}\alpha_{L}^{2}\Xi^{\mu}_{\nu}D_{z}u^{\nu}}{\alpha_{\perp}\alpha_{L}}\,.$
(61)

The gradients $\sigma_{\perp}^{\mu\nu}$ and $\kappa_{\perp}^{\mu}$ are
proportional to the “Navier–Stokes” values for $\pi_{\perp}^{\mu\nu}$ and
$W_{\perp z}^{\mu}$, respectively. After inserting Eq. (60) in Eq. (58) one
obtains

$\pi_{\perp}^{\mu\nu}=2\tau_{r}\beta_{\pi}^{\perp}\sigma_{\perp}^{\mu\nu}\,,\qquad
W_{\perp z}^{\mu}=\tau_{r}\beta_{W}^{\perp}\kappa_{\perp}^{\mu}\,,$ (62)

where

$\beta_{\pi}^{\perp}=\frac{\mathcal{J}_{402-1}}{\alpha_{\perp}^{2}\Lambda}\,,\qquad\beta_{W}^{\perp}=\frac{\mathcal{J}_{421-1}}{\alpha_{L}\alpha_{\perp}\Lambda},$
(63)

and the integrals $\mathcal{J}_{krqs}$ are defined in Appendix A. Thus, the
first-order anisotropic RTA Chapman–Enskog expansion is

$\\!\\!\\!\\!\delta\tilde{f}_{n}=f_{a,n}\bar{f}_{a,n}\left(\frac{p_{\\{\mu}p_{\nu\\}}\pi_{\perp}^{\mu\nu}}{2\beta_{\pi}^{\perp}E_{a,n}\alpha_{\perp}^{2}\Lambda}{-}\frac{(-Z\cdot
p)p_{\\{\mu\\}}W_{\perp
z}^{\mu}}{\beta_{W}^{\perp}E_{a,n}\alpha_{\perp}\alpha_{L}\Lambda}\right)\\!.\\!\\!$
(64)

One can show that Eq. (64) satisfies the anisotropic matching conditions for
$\mathcal{E}$, $\mathcal{P}_{L}$ and $\mathcal{P}_{\perp}$ and the Landau
frame constraints:

$\displaystyle\delta\tilde{\mathcal{E}}$
$\displaystyle=\sum_{n}\int_{p}(u\cdot p)^{2}\delta\tilde{f}_{n}=0\,,$ (65a)
$\displaystyle\delta\tilde{\mathcal{P}}_{L}$
$\displaystyle=\sum_{n}\int_{p}(-Z\cdot p)^{2}\delta\tilde{f}_{n}=0\,,$ (65b)
$\displaystyle\delta\tilde{\mathcal{P}}_{\perp}$
$\displaystyle=\frac{1}{2}\sum_{n}\int_{p}(-p\cdot\Xi\cdot
p)\delta\tilde{f}_{n}=0\,,$ (65c) $\displaystyle Q_{L}$
$\displaystyle=\sum_{n}\int_{p}(u\cdot p)(-Z\cdot p)\delta\tilde{f}_{n}=0\,,$
(65d) $\displaystyle Q^{\mu}_{\perp}$ $\displaystyle=\sum_{n}\int_{p}(u\cdot
p)p^{\\{\mu\\}}\delta\tilde{f}_{n}=0\,,$ (65e)

where $Q_{L}={\,-\,}Z{\,\cdot\,}Q$ is the LRF longitudinal heat flow and
$Q_{\perp}^{\mu}=\Xi^{\mu}_{\nu}Q^{\nu}$ is the transverse heat current.

###  IV.2 Modifying the anisotropic distribution

Figure 3:  The reproduction of selected components of the energy-momentum
tensor $T^{\mu\nu}$, using either the PTMA distribution (68) (dashed color) or
the PTM distribution (32) (dotted color), for a stationary hadron resonance
gas at temperature $T=150$ MeV with zero net baryon density ($\alpha_{B}=0$).
The system is assumed to have a fixed pressure anisotropy
$\pi^{zz}=-\frac{2}{5}\mathcal{P}_{\text{eq}}$ and bulk viscous pressure
$\Pi=-\frac{2}{5}\mathcal{P}_{\text{eq}}$ and is further subjected to either
transverse shear stress (left panel) or longitudinal momentum diffusion (right
panel). The solid black lines show the components of the hydrodynamic input
$T^{\mu\nu}$.

With the anisotropic RTA Chapman–Enskog expansion (64) at hand we can proceed
to modify the leading order anisotropic distribution (50). First, we insert a
perturbation $\delta\Omega$ in $f_{a,n}$, keeping $\delta\Lambda=0$:

$f_{a,n}^{\text{(mod)}}=\frac{g_{n}}{\exp\biggl{[}\dfrac{\sqrt{p{\,\cdot\,}(\Omega{+}\delta\Omega){\,\cdot\,}p{\,+\,}m_{n}^{2}}}{\Lambda}\biggr{]}+\Theta_{n}}\,.$
(66)

Linearizing this modified anisotropic distribution and comparing it to Eq.
(64) one finds

$p\cdot\delta\Omega\cdot
p=-\frac{p_{\\{\mu}p_{\nu\\}}\pi_{\perp}^{\mu\nu}}{\beta_{\pi}^{\perp}\alpha_{\perp}^{2}}+\frac{2(-Z\cdot
p)p_{\\{\mu\\}}W_{\perp
z}^{\mu}}{\beta_{W}^{\perp}\alpha_{\perp}\alpha_{L}}\,.$ (67)

Next we rewrite Eq. (66) as121212The generalization of Eq. (68) for nonzero
baryon chemical potential and diffusion is non-trivial and will be left to
future work.

$f^{\text{PTMA}}_{a,n}=\frac{\mathcal{Z}g_{n}}{\exp\biggl{[}\dfrac{\sqrt{\bm{p}^{\prime\prime
2}{+}m_{n}^{2}}}{\Lambda}\biggr{]}+\Theta_{n}}\,,$ (68)

where the normalization $\mathcal{Z}$ will be discussed further below and

$p_{i}=B_{ij}p^{\prime\prime}_{j}=C_{im}A_{mj}p^{\prime\prime}_{j}\,.$ (69)

By comparing with Eq. (51) one sees that the residual shear transformation
$C_{im}$ further deforms the anisotropic momentum space. The additional
deformations, which are linearly proportional to $\pi_{\perp}^{\mu\nu}$ and
$W_{\perp z}^{\mu}$, are assumed to be smaller than the anisotropy parameters
$\alpha_{L}$ and $\alpha_{\perp}$. The matrix $C_{im}$ is constructed such
that the total deformation matrix $B_{ij}=C_{im}A_{mj}$ is symmetric and
reproduces the perturbation $\delta\Omega$ in Eq. (66) for small residual
shear stresses. After some algebra one finds

$C_{im}=\delta_{im}+\frac{\pi_{\perp,im}}{2\beta_{\pi}^{\perp}}+\frac{\alpha_{\perp}W_{\perp
z,i}Z_{m}+\alpha_{L}W_{\perp
z,m}Z_{i}}{\beta_{W}^{\perp}(\alpha_{\perp}{+}\alpha_{L})}\,.$ (70)

Here $\pi_{\perp,im}=X_{i}{\,\cdot\,}\pi_{\perp}{\cdot\,}X_{m}$ and $W_{\perp
z,i}=-X_{i}{\,\cdot\,}W_{\perp z}$ are the LRF residual shear stress
components. Although $C_{im}$ in Eq. (70) is not symmetric, $B_{ij}$ is:

$B_{ij}=A_{ij}+\frac{\alpha_{\perp}\pi_{\perp,ij}}{2\beta_{\pi}^{\perp}}+\frac{\alpha_{\perp}\alpha_{L}(W_{\perp
z,i}Z_{j}+W_{\perp
z,j}Z_{i})}{\beta_{W}^{\perp}(\alpha_{\perp}{+}\alpha_{L})}.$ (71)

Finally we renormalize the particle density to the one given by the leading
order anisotropic distribution:

$n_{a,n}={\det}A\cdot n_{\text{eq},n}(\Lambda,0)\,$ (72)

where $n_{a,n}$ is the anisotropic particle density and
${\det}A=\alpha_{\perp}^{2}\alpha_{L}$ (there is no contribution from the
anisotropic RTA Chapman–Enskog correction $\delta\tilde{f}_{n}$). The modified
particle density from Eq. (68) is

$n_{a,n}^{\text{PTMA}}=\mathcal{Z}\,{\det}C\cdot n_{a,n}\,,$ (73)

which leads us to the normalization factor

$\mathcal{Z}=\frac{1}{{\det}C}\,.$ (74)

In Figure 3, we test the reproduction of the energy-momentum tensor by the
“PTMA distribution” (68), similar to the test done in Fig. 1. We consider a
hadron resonance gas at temperature $T=150$ MeV and baryon chemical potential
$\alpha_{B}=0$ with an input $T^{\mu\nu}$ featuring the two fixed dissipative
flows $\pi^{zz}=\Pi=-\frac{2}{5}\mathcal{P}_{\text{eq}}$ (or, equivalently,
$\mathcal{P}_{L}=\frac{1}{5}\mathcal{P}_{\text{eq}}$ and
$\mathcal{P}_{\perp}=\frac{4}{5}\mathcal{P}_{\text{eq}}$). These viscous
pressures are captured by the leading order distribution $f_{a,n}$ by
numerically adjusting (“Landau matching”) the anisotropy parameters $\Lambda$,
$\alpha_{L}$ and $\alpha_{\perp}$ accordingly.

In the left panel of Fig. 3 we explore the reliability of the modified
anisotropic distribution in reproducing a nonzero transverse shear stress
$\pi_{\perp}^{xx}=-\pi_{\perp}^{yy}$. Moving from the right edge of the plot
leftward we decrease $\pi_{\perp}^{xx}$, causing the pressures
$\mathcal{P}_{x}=\mathcal{P}_{\perp}+\pi_{\perp}^{xx}$ and
$\mathcal{P}_{y}=\mathcal{P}_{\perp}+\pi_{\perp}^{yy}$ to decrease and
increase, respectively. According to the anisotropic matching conditions (65),
the energy density and longitudinal pressure should remain constant as we do
so. The hydrodynamic input (solid black lines) and kinetic output $T^{\mu\nu}$
components (colored dashed lines) from this test are very similar to those
shown in the left panels of Fig. 1, with the kinetic output for
$\mathcal{P}_{x}$ (blue dashed) approaching zero for large negative
$\pi_{\perp}^{xx}$ while the output for $\mathcal{P}_{y}$ (green dashed)
overestimates the hydrodynamic target value but to a lesser degree. The
kinetic outputs for $\mathcal{E}$ and $\mathcal{P}_{L}$ are in very good
agreement with their hydrodynamic target values. The modified anisotropic
distribution breaks down for $\pi_{\perp}^{xx}<-2\beta_{\pi}^{\perp}$ (or
$\det C<0$) but for moderately large values
$|\pi_{\perp}^{xx}|\leq\frac{1}{6}(\mathcal{P}_{L}{+}2\mathcal{P}_{\perp})$
the discrepancies between the input and output values are
$|\Delta\mathcal{P}_{x}|/\mathcal{P}_{x}\leq 4.0\%$,
$|\Delta\mathcal{P}_{y}|/\mathcal{P}_{y}\leq 1.5\%$,
$|\Delta\mathcal{P}_{L}|/\mathcal{P}_{L}\leq 0.17\%$ and
$|\Delta\mathcal{E}|/\mathcal{E}\leq 0.47\%$. We also checked that for this
hydrodynamic input the off-diagonal components of the kinetic output
$T^{\mu\nu}$ are zero.

For the second test (right panel of Fig. 3) we start from the same system but
add a longitudinal momentum diffusion component $W_{\perp z}^{x}$ of
increasing magnitude. The energy density and pressure components
$\mathcal{P}_{x}$, $\mathcal{P}_{y}$ and $\mathcal{P}_{L}$ should not change
as we do so. Overall, we find that the kinetic outputs for $\mathcal{E}$,
$\mathcal{P}_{y}$ and $W_{\perp z}^{x}$ are in good agreement with this
expectation. However, the kinetic outputs for the pressure components
$\mathcal{P}_{x}$ and $\mathcal{P}_{L}$ are seen to exhibit stronger
sensitivity to the errors caused by mismatched nonlinear terms in $W_{\perp
z}^{x}$. Still, for moderate values of $|W_{\perp
z}^{x}|\leq\frac{1}{6}(\mathcal{P}_{L}{+}2\mathcal{P}_{\perp})$, the output
errors stay below $|\Delta W_{\perp,z}^{x}/W_{\perp,z}^{x}|\leq 1.6\%$,
$|\Delta\mathcal{P}_{x}|/\mathcal{P}_{x}\leq 4.0\%$,
$|\Delta\mathcal{P}_{y}|/\mathcal{P}_{y}\leq 0.55\%$,
$|\Delta\mathcal{P}_{L}|/\mathcal{P}_{L}\leq 21\%$ and
$|\Delta\mathcal{E}|/\mathcal{E}\leq 0.84\%$

Finally, we repeat the previous tests with the PTM modified equilibrium
distribution (32) and plot the corresponding kinetic output $T^{\mu\nu}$ as
dotted lines in Fig. 3 for comparison. One observes that the PTM distribution
follows essentially the same trends as the PTMA distribution, but it does not
fully capture the pressure anisotropy and bulk viscous pressure at zero
residual shear stress, similar to what we saw in Fig. 1. On the other hand, by
imposing the generalized Landau matching conditions (65b,c) the leading order
anisotropic distribution precisely reproduces the longitudinal and transverse
pressures in that limit.

Figure 4:  The $\tau{-}x$ slice $(y=\eta_{s}=0)$ of the shear inverse Reynolds
number
$\text{Re}^{-1}_{\pi}=\sqrt{\pi^{\mu\nu}\pi_{\mu\nu}}/(\mathcal{P}_{\text{eq}}\sqrt{3})$
(left panels) and bulk inverse Reynolds number
$\text{Re}^{-1}_{\Pi}=|\Pi|/\mathcal{P}_{\text{eq}}$ (right panels) showing
the strength of the $\delta f_{n}$ corrections on a particlization
hypersurface of constant temperature $T_{\text{sw}}=150$ MeV (white contour)
from a (2+1)–d central Pb+Pb collision with smooth T${}_{\text{\sc R}}$ENTo
initial conditions. We vary the peak temperature of the specific bulk
viscosity $\zeta/\mathcal{S}$ to either $T_{\zeta}=220$ MeV (top row) or
$T_{\zeta}=160$ MeV (bottom row).

##  V Continuous particle spectra

In this section we compute the continuous momentum spectra of identified
hadrons $(\pi^{+},K^{+},p)$ using the Cooper–Frye formula (1). To generate the
hypersurfaces, we run the VAH code McNelis  _et al._ (2021a) to evolve central
and non-central Pb+Pb collisions using standard viscous hydrodynamics with
smooth initial conditions. For simplicity, we only consider the central slice
($\eta_{s}=0$) in the transverse plane and assume longitudinal boost-
invariance to extend the solution in the spacetime rapidity direction. For the
different hadron phase-space distribution models discussed in Secs. II – IV,
we will compare the azimuthally averaged transverse momentum spectra

$\frac{dN_{n}}{2\pi
p_{T}dp_{T}dy_{p}}=\int_{0}^{2\pi}\,\frac{d\phi_{p}}{2\pi}\frac{dN_{n}}{p_{T}dp_{T}d\phi_{p}dy_{p}}$
(75)

and the $p_{T}$–differential elliptic flow coefficient

$v^{(2)}_{2}(p_{T})=\frac{\int_{0}^{2\pi}d\phi_{p}\cos(2\phi_{p})\,\dfrac{dN_{n}}{p_{T}dp_{T}d\phi_{p}dy_{p}}}{\int_{0}^{2\pi}d\phi_{p}\,\dfrac{dN_{n}}{p_{T}dp_{T}d\phi_{p}dy_{p}}}$
(76)

at mid-rapidity ($y_{p}=0$). Specifically, we compare results obtained with
the 14–moment approximation (8), the RTA Chapman–Enskog expansion (22), the
PTM and PTB modified equilibrium distributions (32) and (39), and the PTMA
modified anisotropic distribution (68).

###  V.1 Setup

We evolve an azimuthally symmetric, event-averaged T${}_{\text{\sc R}}$ENTo
transverse energy density profile with (2+1)–dimensional second-order viscous
hydrodynamics Moreland  _et al._ (2015); McNelis  _et al._ (2021a). We start
the simulation at the longitudinal proper time $\tau_{0}=0.5$ fm/$c$ with an
initial central temperature of $T_{0,\text{center}}=400$ MeV. The spatial
components of the fluid velocity $u^{\mu}$ and bulk viscous pressure $\Pi$ are
initialized to zero. The shear stress tensor is initialized as

$\pi^{\mu\nu}=\frac{1}{3}\big{(}\mathcal{P}_{L}{-}\mathcal{P}_{\perp}\big{)}\left(\Delta^{\mu\nu}+3Z^{\mu}Z^{\nu}\right)\,,$
(77)

where the initial pressure anisotropy is set to
$\mathcal{P}_{L}-\mathcal{P}_{\perp}=\frac{12}{11}\mathcal{P}_{\text{eq}}$ and
the initial longitudinal basis vector is $Z^{\mu}=(0,0,0,\tau_{0}^{-1})$. For
the equilibrium pressure $\mathcal{P}_{\text{eq}}(\mathcal{E})$, we use the
lattice QCD equation of state from the HotQCD collaboration Bazavov  _et al._
(2014). The baryon chemical potential $\alpha_{B}$ and baryon diffusion
current $V_{B}^{\mu}$ are fixed to zero for the entire simulation.131313The
current version of the relativistic hydrodynamic code VAH McNelis  _et al._
(2021a) does not propagate the net baryon density and baryon diffusion
current.

The shear and bulk viscosities are modeled using the temperature dependent
parametrizations from the JETSCAPE collaboration Everett  _et al._ (2020a, b):

$\displaystyle(\eta/\mathcal{S})(T)=$
$\displaystyle\,(\eta/\mathcal{S})_{\text{kink}}\,+\,a_{\text{low}}\,(T{-}T_{\eta})\Theta(T_{\eta}{-}T)$
(78a)
$\displaystyle\,+\,a_{\text{high}}(T{-}T_{\eta})\,\Theta(T{-}T_{\eta})\,,$
$\displaystyle(\zeta/\mathcal{S})(T)=$
$\displaystyle\,\frac{(\zeta/\mathcal{S})_{\text{max}}\,\Lambda_{\zeta}(T)^{2}}{\Lambda_{\zeta}(T)^{2}+(T{-}T_{\zeta})^{2}}\,,$
(78b)

where $\Theta(x)$ is the Heaviside step function and

$\Lambda_{\zeta}(T)=w_{\zeta}\left(1+\lambda_{\zeta}\,\mathrm{sgn}(T{-}T_{\zeta})\right)\,.$
(79)

In this work we fix the viscosity parameters to
$(\eta/\mathcal{S})_{\text{kink}}=0.093$, $T_{\eta}=0.22$ GeV,
$a_{\text{low}}=-0.77$ GeV-1, $a_{\text{high}}=0.21$ GeV-1,
$(\zeta/\mathcal{S})_{\text{max}}=0.1$, $w_{\zeta}=0.05$ GeV and
$\lambda_{\zeta}=0$. However, for exploration purposes we vary the temperature
$T_{\zeta}$ at which $\zeta/\mathcal{S}$ peaks, setting it to either $220$ MeV
or $160$ MeV.

From the hydrodynamic simulation we generate an isothermal particlization
hypersurface of temperature $T_{\text{sw}}=150$ MeV, using the freeze-out
finder code CORNELIUS Huovinen and Petersen (2012). For our central Pb+Pb
collision, the longitudinally boost-invariant hypersurface at spacetime
rapidity $\eta_{s}=0$ contains about $6.8{\,\times\,}10^{4}$ freeze-out cells.
Figure 4 shows a $\tau{-}x$ slice of the particlization surface at
$y=\eta_{s}=0$, as well as the shear and bulk inverse Reynolds numbers to
gauge the strength of the shear and bulk $\delta f_{n}$ corrections.141414We
note the setup used in this Section is different from the one in an earlier
study to test the iS3D particle sampler (see Sec. 3 in Ref. McNelis  _et al._
(2021b)). There we found that the viscous corrections were so large that the
PTMA distribution was unusable for a significant part of the hypersurface. The
hypersurface shown in Fig. 4 was computed with parameter settings that allow
for a meaningful comparison between the PTM and PTMA distributions.  Finally,
we evaluate the Cooper–Frye formula with the code iS3D McNelis  _et al._
(2021b), which was developed from the particlization module iSS Shen  _et al._
(2016). The code gives the user the option to use one of five $\delta f_{n}$
corrections described in Secs. II – IV. To compute the longitudinally boost-
invariant particle spectra, the hypersurface volume needs to be extended to
the $\eta_{s}$ dimension and centered around the momentum rapidity $y_{p}$. We
perform the numerical integration of the Cooper–Frye formula along the
$\eta_{s}$–direction using Gauss–Legendre integration on a 48-point grid
$(y_{p}{\,-\,}\eta_{s})_{j}$ with integration weights $\omega_{j}$ given by

$\displaystyle(y_{p}{\,-\,}\eta_{s})_{j}$
$\displaystyle=\sinh^{-1}\left(\frac{x_{j}}{1{\,-\,}x_{j}^{2}}\right)\,,$
(80a) $\displaystyle\omega_{j}$
$\displaystyle=w_{j}\frac{1{\,+\,}x_{j}^{2}}{|1{\,-\,}x_{j}^{2}|\sqrt{1{\,-\,}x_{j}^{2}{\,+\,}x_{j}^{4}}}\,,$
(80b)

where $x_{j}$ and $w_{j}$ are the Gauss–Legendre roots and weights,
respectively.

###  V.2 Breakdown of the modified distributions and technical issues

As stated in the previous sections, the modified equilibrium distributions
(32) and (39) break down in freeze-out cells with large viscous corrections.
For the first scenario where the specific bulk viscosity’s peak temperature is
$T_{\zeta}=220$ MeV (top row in Fig. 4), both the shear and bulk viscous
corrections are small enough to use the modified equilibrium distribution for
all freeze-out cells. For $T_{\zeta}=160$ MeV (bottom row in Fig. 4) the bulk
viscosity peaks much closer to the particlization hypersurface. Not only are
the bulk viscous corrections much larger in this case, but the shear
corrections are also enhanced due to shear-bulk coupling effects in the
hydrodynamic simulation Denicol  _et al._ (2014). Together they cause the
modified equilibrium distribution to break down for about 4400 freeze-out
cells at $r{\,\sim\,}10$ fm and $0.5$ fm/$c{\,<\,}\tau{\,<\,}1.3$ fm/$c$.
There are several options to handle the particle spectra contributions from
these freeze-out cells: (i) ignore such cells entirely, (ii) use for them the
local-equilibrium distribution $f_{n}=f_{\text{eq},n}$, which is positive
definite but neglects the viscous components of $T^{\mu\nu}$, or (iii)
linearize in such cells the modified equilibrium distribution
$f_{\text{eq},n}^{(\text{mod})}\approx f_{\text{eq},n}+\delta f_{n}$, which
may turn negative at high momenta but captures all components of $T^{\mu\nu}$.
Here we choose the third option, i.e. for the PTM distribution we switch to
the Chapman–Enskog expansion (22) when ${\det}A<10^{-5}$ or when the
normalization factor $\mathcal{Z}_{n}$ of any hadron (usually the lightest
pion $\pi^{0}$) turns negative. For the PTB distribution, we use $\delta
f_{n}$ from Eq. (46) when ${\det}A<10^{-5}$ or
$-\mathcal{P}_{\text{eq}}<\Pi<\mathcal{E}/3-\mathcal{P}_{\text{eq}}$. In this
study, the impact of these manipulations on the total particle spectra is
negligible since the number of “bad” freeze-out cells is small. However, if a
significant fraction of the hypersurface requires one to use something other
than the modified equilibrium distribution, then Eqs. (32) and (39) should
probably not be used for particlization.

Figure 5:  The azimuthally-averaged transverse momentum spectra of
$(\pi^{+},K^{+},p)$ (top panels) computed with the Cooper–Frye formula for the
(2+1)–d central Pb+Pb collision described in Sec. V.1. The peak temperature of
$\zeta/\mathcal{S}$ is set to $T_{\zeta}=220$ MeV. We compare the shear (red),
bulk (blue) and combined shear + bulk (purple) $\delta f_{n}$ corrections of
the 14–moment approximation (dotted color), RTA Chapman–Enskog expansion
(dashed color), PTM distribution (solid color) and PTB distribution (dot-
dashed color) relative to the ideal spectra with $\delta f_{n}=0$ (solid
black). The bottom panels show the ratio of the particle spectra with $\delta
f_{n}$ corrections to the ideal spectra.

Even if the modified equilibrium distribution can be used, a very small
${\det}A$ value can cause its width along the rapidity direction
$y_{p}{\,-\,}\eta_{s}$ to be extremely narrow. This leads to numerical errors
in the particle yields if the spacetime rapidity grid
$(y_{p}{\,-\,}\eta_{s})_{j}$ is not fine enough. To resolve this technical
issue, we rescale the spatial grid $(y_{p}{\,-\,}\eta_{s})_{j}$ by the
distribution’s rapidity width $\delta y_{p}$. One can estimate $\delta y_{p}$
from a diagonal deformation matrix
$A_{ij}{\,=\,}\mathop{\mathrm{diag}}(1{\,+\,}\bar{\Pi}{\,-\,}\frac{1}{2}\bar{\pi}_{zz},1{\,+\,}\bar{\Pi}{\,-\,}\frac{1}{2}\bar{\pi}_{zz},1{\,+\,}\bar{\Pi}{\,+\,}\bar{\pi}_{zz})$:

$\delta y_{p}\sim 1+\bar{\Pi}+\bar{\pi}_{zz}\,,$ (81)

where $\bar{\Pi}=\Pi/(3\beta_{\Pi})$ (or $\lambda_{\Pi}$) and
$\bar{\pi}_{zz}=\pi_{zz}/(2\beta_{\pi})$. Relative to the transverse momentum
space, the rapidity distribution becomes very narrow when
$\bar{\pi}_{zz}\approx-(1+\bar{\Pi})$. In this limit, the rapidity width is
proportional to

$\delta y_{p}\propto\frac{{\det}A}{({\det}A_{\Pi})^{2/3}}\,,$ (82)

where ${\det}A_{\Pi}=(1{\,+\,}\bar{\Pi})^{3}$. Therefore, for each
longitudinally boost-invariant freeze-out cell $d^{2}\sigma_{\mu,i}$ we
rescale the spacetime rapidity grid points and weights (80) as

$\displaystyle(y_{p}{\,-\,}\eta_{s})_{i,j}$
$\displaystyle=\frac{{\det}A_{i}}{({\det}A_{\Pi,i})^{2/3}}\times(y_{p}{\,-\,}\eta_{s})_{j}\,,$
(83a) $\displaystyle\omega_{i,j}$
$\displaystyle=\frac{{\det}A_{i}}{({\det}A_{\Pi,i})^{2/3}}\times\omega_{j}\,.$
(83b)

This rescaling trick is found to work well even for small values of
${\det}A\sim 10^{-5}$. For (3+1)–d hypersurfaces, however, this method cannot
be used because the freeze-out finder fixes the freeze-out cells’ spacetime
rapidity. Instead, we switch to a linearized $\delta f_{n}$ correction if the
rapidity width is too small (e.g. ${\det}A/({\det}A_{\Pi})^{2/3}<0.01$).

The modified anisotropic distribution (68) can also break down, usually
because the longitudinal pressure $\mathcal{P}_{L}$ turns negative during the
viscous hydrodynamic simulation. When this happens one cannot construct a
solution for the momentum deformation parameter $\alpha_{L}$ in the leading
order anisotropic distribution (48). For the hypersurface with large bulk
viscous corrections we find $\mathcal{P}_{L}{\,<\,}0$ in about 7800 freeze-out
cells during the period $0.5$ fm/$c<\tau<2.3$ fm/$c$ (see footnote 14). For
these freeze-out cells, we here simply replace Eq. (68) by the local-
equilibrium distribution (2). In practice it would be more appropriate to use
the PTMA distribution on hypersurfaces constructed from anisotropic
hydrodynamic simulations in which the occurrence of negative longitudinal
pressures is largely avoided McNelis  _et al._ (2021a).

Figure 6:  Same as Fig. 5 but with a $\zeta/\mathcal{S}$ peak temperature of
$T_{\zeta}=160$ MeV.

###  V.3 Central collisions

Figure 5 shows the continuous transverse momentum spectra (75) of ($\pi^{+}$,
$K^{+}$, $p$), without resonance decay contributions or hadronic rescattering,
for our central Pb+Pb collision, using a specific bulk viscosity peak
temperature of $T_{\zeta}=220$ MeV. We study the shear and bulk viscous
corrections to the ideal spectra (i.e. $\delta f_{n}=0$) computed with the
four models for $\delta f_{n}$ discussed in Sec. II.

As seen in the figure, the shear stress (although small) slightly flattens the
$p_{T}$ spectra (red curves), without affecting the total yields, since it
slows the longitudinal expansion and pushes the particles outward in the
transverse direction. The bulk viscous pressure has the opposite effect by
counteracting the scalar expansion rate and reducing the average pressure,
thereby softening the $p_{T}$ spectra (blue curves). One observes the
“shoulder” at low values of $p_{T}$ that is typical for thermal flow spectra
Schnedermann  _et al._ (1993); there the bulk corrections are more pronounced
in protons than in pions and kaons. Furthermore, the bulk viscous correction
decreases the total pion and kaon yields while increasing the proton yield.
The PTB distribution is the only exception to these trends: for non-zero bulk
viscous pressure it both raises the shoulder and increases the yields of all
particles, which sets it apart from the other $\delta f_{n}$ corrections.
Finally, the purple curves show the combined effects of the shear and bulk
viscous corrections to the spectra. Overall, the slope of the total $p_{T}$
spectra is steeper than the ideal one since here the bulk viscous pressure is
larger than the shear stress for most of the hypersurface.

We also compare the $\delta f_{n}$ corrections of the 14–moment approximation,
RTA Chapman–Enskog expansion, PTM distribution and PTB distribution. We see
that the momentum dependence of the viscous correction varies considerably
between models: the ratio of the total particle spectra to the ideal one at
high $p_{T}$ is approximately quadratic for the 14–moment approximation and
linear for the RTA Chapman–Enskog expansion, as illustrated by the bottom
panels in Fig. 5. The PTM spectra are almost identical to the ones computed
with the Chapman–Enskog expansion because the shear stress and bulk viscous
pressures are small; differences between the two approaches emerge only at
large values of $p_{T}$. Overall, there are no significant differences among
the first three $\delta f_{n}$ models across most of the $p_{T}$ spectrum in
Fig. 5. The PTB distribution, on the other hand, has a moderate excess of low
$p_{T}$ pions and kaons relative to the other $\delta f_{n}$ corrections.

Figure 7: (Color online) The $p_{T}$–differential elliptic flow coefficient of
$(\pi^{+},K^{+},p)$ for the (2+1)–d non-central Pb+Pb collision (top panels).
The peak temperature of $\zeta/\mathcal{S}$ is set to $T_{\zeta}=220$ MeV.
Similar to Fig. 5, we plot the shear and bulk viscous corrections of each
$\delta f_{n}$ model to the ideal $v_{2}(p_{T})$ (solid black). The bottom
panels show the differences between $v_{2}(p_{T})$ with $\delta f_{n}$
corrections and the ideal $v_{2}(p_{T})$.

In Figure 6 we show how the central Pb+Pb collision spectra from Fig. 5 change
when we move the peak of the specific bulk viscosity from $T_{\zeta}=220$ MeV
to 160 MeV. This obviously increases the magnitude of the bulk viscous
pressure on the particlization hypersurface at $T_{\mathrm{sw}}=150$ MeV. The
shear-bulk coupling effect in the hydrodynamic simulation then also increases
the strength of the shear viscous corrections on the hypersurface, resulting
in much flatter $p_{T}$ slopes than in the previous case. Surprisingly, the
shear viscous corrections from the four $\delta f_{n}$ models are still very
close to each other (for PTM and PTB this is expected because these
distributions use the same shear stress modification.) Instead, the largest
differences are found in their bulk viscous corrections. Because the bulk
viscous pressure on the hypersurface is quite large (lower right panel in Fig.
4), the linearized bulk viscous correction in the 14–moment approximation
causes the pion spectra to turn negative already at intermediate momentum
$p_{T}\sim 1.8$ GeV while the kaon and proton spectra flip sign for $p_{T}>2$
GeV and 2.4 GeV, respectively. The bulk viscous correction in the RTA
Chapman–Enskog expansion is less severe than in the 14–moment approximation,
but the pion and kaon spectra still turn negative for $p_{T}\gtrsim 2$ GeV
(the proton spectra remain positive up to $p_{T}=3$ GeV). In contrast, the PTM
spectra are positive definite by construction even for moderately large bulk
viscous pressures, maintaining an exponential tail at high values of
$p_{T}$.151515Without regulation, the linearized $\delta f_{n}$ corrections
from the “bad” freeze-out cells (see Sec. V.2) eventually turn the PTM and PTB
spectra negative but this occurs outside the experimental range of interest
for soft hadron emission $p_{T}<3$ GeV. One also observes a slight excess of
protons at low $p_{T}$ (this effect is much less visible for kaons). The PTB
spectra also remain positive but have significantly steeper slopes than the
PTM spectra; the shoulder enhancements at low $p_{T}$ are also much larger,
especially for pions.

After combining the shear and bulk viscous corrections, the $p_{T}$ spectra
become positive except for the one computed with the 14–moment approximation,
which has a strong quadratic momentum dependence from its bulk viscous
correction.

Figure 8:  Same as Fig. 7 but with a $\zeta/\mathcal{S}$ peak temperature of
$T_{\zeta}=160$ MeV.

###  V.4 Non-central collisions

Next we repeat the same hydrodynamic simulations for a nonzero impact
parameter $b=5$ fm, to study the effects of the different viscous corrections
$\delta f_{n}$ on the $p_{T}$–differential elliptic flow coefficient, shown in
Fig. 7 for $\pi^{+}$, $K^{+}$, and $p$. Here we start with $T_{\zeta}=220$
MeV, resulting in relatively weak viscous stresses on the particlization
hypersurface. The shear viscous corrections are seen to decrease the
differential elliptic flow $v_{2}(p_{T})$, counteracting the effects of
anisotropic transverse flow. The bulk viscous pressure, on the other hand,
tends to increase $v_{2}(p_{T})$ by suppressing the radial flow and making the
$p_{T}$ spectra (75) steeper. Individually, the small shear and bulk viscous
corrections to the ideal $v_{2}(p_{T})$ (defined by setting $\delta
f_{n}{\,=\,}0$) are roughly linear in $p_{T}$, except for the bulk $\delta
f_{n}$ correction of the 14–moment approximation, which causes $v_{2}(p_{T})$
to diverge at high $p_{T}$ when the spectrum (which enters the denominator of
$v_{n}(p_{T})$) passes through zero. Overall, there is a net suppression on
the differential elliptic flow since it is more sensitive to the shear viscous
corrections at the space-like edges of the hypersurface, whose fluid cells
have undergone the strongest transverse acceleration. Interestingly, the only
exception is the PTB elliptic flow, whose shear and bulk viscous corrections
nearly cancel each other. Similar to the spectra discussed in the previous
subsection, the four $\delta f_{n}$ models produce very similar viscous
corrections to the elliptic flow when $T_{\zeta}=220$ MeV. Clear distinctions
between these models appear only for large transverse momenta $p_{T}>2$ GeV.

When $T_{\zeta}$ is lowered to 160 MeV (see Fig. 8), the much larger bulk
viscous pressure on the space-like part of the hypersurface results in the
divergence of the linearized bulk viscous corrections to $v_{2}(p_{T})$ at
intermediate values of $p_{T}{\,\sim\,}1.5{-}2.5$ GeV, for all three particle
species considered, reflecting the corresponding sign change of their
azimuthally averaged $p_{T}$ spectra (75). The shear viscous corrections help
offset this effect by flattening the slope of the $p_{T}$ spectra, but the
pion and kaon elliptic flows still diverge, albeit now at slightly higher
$p_{T}$ (as does the proton $v_{2}(p_{T})$ from the 14–moment approximation).
For realistic event-by-event simulations where the hadrons are Monte Carlo
sampled from the hypersurface prior to the afterburner phase, these
divergences are removed by enforcing the regulation $f_{\text{eq},n}+\delta
f_{n}\geq 0$ in the Cooper–Frye formula Shen  _et al._ (2016); McNelis  _et
al._ (2021b).

In the PTM distribution, on the other hand, the bulk viscous modifications
prevent the $p_{T}$ spectra from turning negative, and the resulting elliptic
flow curves are well behaved even at high $p_{T}$. After including the shear
stress modification, the PTM differential elliptic flows decrease but remain
above the ideal curves since the bulk viscous pressure overwhelms the shear
stress on the hypersurface. The PTB elliptic flow coefficients also stay
finite at high $p_{T}$ but are significantly larger than those computed with
the PTM distribution.

Figure 9:  The azimuthally-averaged transverse momentum spectra of
($\pi^{+},K^{+},p)$ (top panels) for the central Pb+Pb collision with
$T_{\zeta}=160$ MeV. We compare the full PTMA distribution (solid light blue)
to the PTM distribution with shear and bulk modifications (solid purple), the
leading order anisotropic (LOA) distribution (48) without residual shear
corrections (dotted light blue), and the ideal (local thermal equilibrium)
distribution (solid black). The bottom panels show the ratio of the particle
spectra with viscous corrections to the ideal spectra.

Figure 10:  The $p_{T}$–differential elliptic flow of ($\pi^{+},K^{+},p)$ for
the non-central Pb+Pb collision with $T_{\zeta}=160$ MeV. Similar to Fig. 10,
we compare the PTMA distribution (solid light blue) to the PTM distribution
(solid purple) and leading order anisotropic (LOA) distribution (dotted light
blue).

###  V.5 Modified anisotropic distribution

We close this section by looking at the particle $p_{T}$ spectra and
differential elliptic flows computed with the PTMA distribution (68). Figure
10 shows the resulting transverse momentum spectra for the same central Pb+Pb
collision as in Fig. 6. Compared to the PTM spectra, which here include both
shear and bulk viscous corrections, the PTMA kaon and proton distributions
have a slightly higher mean transverse momentum, indicated by corresponding
shifts in the slope and shoulder, and the pion yield slightly increases. We
know from Fig. 1 that the PTM distribution underpredicts the input isotropic
pressure $\mathcal{P}_{\text{eq}}+\Pi$ because the bulk viscous modifications
to the effective temperature and isotropic momentum space in Eq. (32) do not
perfectly reproduce the input bulk viscous pressure. The anisotropic
parameters ($\Lambda,\alpha_{\perp},\alpha_{L}$) in the PTMA distribution are
optimized to correct for these errors, outputting a slightly larger isotropic
pressure than the PTM distribution. One also notices that the PTMA spectra are
virtually identical to the ones computed with leading order anisotropic
distribution (48), which excludes the residual shear corrections in Eq. (71).
Due to the approximate azimuthal symmetry of the hypersurface, we expect the
transverse shear stress $\pi_{\perp}^{\mu\nu}$ ($W_{\perp z}^{\mu}=0$ by
longitudinal boost-invariance) to have little to no impact on the azimuthally-
averaged particle spectra.

Figure 10 shows the elliptic flow coefficients from the PTMA distribution for
the same non-central Pb+Pb collision as in Fig. 8. Although the PTM and PTMA
models yield very similar results, the smaller bulk viscous modification in
the PTMA distribution slightly brings down the $v_{2}(p_{T})$ curves. Relative
to the leading order anisotropic distribution, we see that the transverse
shear corrections play a larger role in non-central collisions, damping the
anisotropic flow.

##  VI Conclusions

In this work we proposed and studied a positive definite hadronic distribution
for Cooper–Frye particlization that improves upon the Pratt–Torrieri
distribution introduced in Ref. Pratt and Torrieri (2010). At an intermediate
stage of the derivation, we used the RTA Boltzmann equation to identify the
non-equilibrium corrections to the Boltzmann factor. This ensures that in the
limit of small dissipative flows our “modified equilibrium” distribution
function reduces exactly to the first-order RTA Chapman–Enskog expansion. Even
with these dissipative modifications, the distribution function remains
positive definite for arbitrarily large momenta, but there is a trade-off:
different from the first-order RTA Chapman–Enskog expansion, the modified
equilibrium distribution no longer matches the input energy-momentum tensor
and net baryon current exactly. The mismatch becomes significant already for
moderately large viscous corrections.

To minimize these errors, in Sec. III.2 we slightly restructured the
dissipative perturbations by linearizing the viscous corrections to the
momentum scales of the local-equilibrium distribution as well as their
contributions to the corresponding particle yields. The resulting PTM modified
equilibrium distribution bears close resemblance to the Pratt–Torrieri
distribution Pratt and Torrieri (2010) but it can better reproduce the input
hydrodynamic quantities in freeze-out cells subject to moderately large bulk
viscous pressures. Although the output of the energy-momentum tensor
calculated with the modified equilibrium distribution does not perfectly match
the input $T^{\mu\nu}$, one can further improve it by applying the same
technique to modify not the isotropic local-equilibrium distribution but
instead the leading order anisotropic distribution (which accounts for some
large dissipative effects non-perturbatively). At this point this PTMA
distribution has only been constructed for systems with zero net baryon
density; its generalization to nonzero baryon chemical potential will be left
to future work.

We also compared the PTM distribution to the linearized 14–moment
approximation and the first-order RTA Chapman–Enskog expansion, by using the
Cooper–Frye formula to compute the momentum spectra and $p_{T}$–differential
elliptic flow coefficients of hadrons emitted from a particlization
hypersurface. For small viscous corrections the azimuthally averaged
transverse momentum spectra and $p_{T}$–differential elliptic flows generated
by all of these $\delta f_{n}$ models are very similar. When the dissipative
flows on the hypersurface (in our case mostly the bulk viscous pressure)
become moderately large, the transverse momentum spectra of the linearized
$\delta f_{n}$ corrections turn negative at intermediate $p_{T}$ values while
those of the PTM and PTMA distributions remain positive at all $p_{T}$, by
construction. This simultaneously prevents the $p_{T}$–differential elliptic
flow coefficients from developing singularities caused by zero crossings of
the $p_{T}$ spectra, a problem that plagues the linearized viscous corrections
whenever the hypersurface features large bulk viscous pressures.

In practical applications the hydrodynamic evolution of the fluid inside the
particlization hypersurface must be complemented by a microscopic kinetic
evolution outside that surface, typically simulated by a hadronic cascade that
allows for the decay of unstable resonances and their regeneration via
hadronic rescattering (such as e.g. URQMD Bleicher  _et al._ (1999) or SMASH
Weil  _et al._ (2016)). Initialization of this cascade requires the sampling
of hadron positions and momenta from the Cooper–Frye formula using Monte Carlo
techniques Shen (7 25), rather than computing continuous particle spectra by
numerically performing the Cooper–Frye integral as we have done in this work.
Monte Carlo sampling requires a positive definite hadron distribution function
$f_{n}=f_{\text{eq},n}+\delta f_{n}$. Positivity during the particle sampling
process can be ensured by appropriate regularization prescriptions Shen (7
25); McNelis  _et al._ (2021b). Combined with fluctuations from finite
sampling statistics, it then becomes difficult to distinguish the modified
equilibrium distributions from the regulated linearized $\delta f_{n}$
corrections at large momenta where few particles are emitted in any given
collision event. The differences that matter in practice are those that are
visible at soft momenta $p_{T}\lesssim 1$ GeV/$c$.

The Monte Carlo sampling of hadrons from the Cooper–Frye formula using one of
the five $\delta f_{n}$ models discussed in this work has already been
implemented in the code iS3D McNelis  _et al._ (2021b). In addition to the PTM
and PTMA distributions, it may be insightful to consider other $\delta f_{n}$
correction candidates, in order to fully quantify the theoretical
uncertainties associated with the particlization process in heavy-ion
collisions. A recent model-to-data analysis that utilized our code to reduce
model selection bias in the phenomenological constraints on the quark-gluon
plasma’s transport properties was reported in Refs. Everett  _et al._ (2020a,
b). The only $\delta f_{n}$ model discussed here that was not used in that
analysis is the PTMA modified anisotropic distribution, which is the latest
addition to the iS3D particle sampler.

An interesting recent development Everett  _et al._ (2021) is the proposal of
a new type of positive definite particle distribution that introduces
dissipative corrections to the local thermal equilibrium distribution by
maximizing the entropy while using additional Lagrange multipliers to keep all
ten components of the energy-momentum tensor $T^{\mu\nu}$ (including the bulk
and shear viscous stresses) fixed exactly at their hydrodynamic values. The
conceptual advantage of this approach is that it does not rely on any
uncontrolled assumptions about the microscopic dynamics of the distribution
function but uses only the available macroscopic information encoded in the
hydrodynamic output for $T^{\mu\nu}$. We note that this so-called “maximum
entropy” distribution Everett  _et al._ (2021) shares many features with our
PTM and PTMA modified equilibrium distributions. It would be interesting to
see whether a deeper connection exists between these two approaches, and to
compare their particle spectra in event-by-event simulations of heavy-ion
collisions with fluctuating initial conditions. Such studies will become
feasible once a full numerical calculation of the Lagrange multipliers in the
maximum entropy distribution for arbitrary energy-momentum tensors becomes
available.

###### Acknowledgements.

The authors would like to thank Scott Pratt, Jonah Bernhard, Derek Everett and
Chandrodoy Chattopadhyay for useful discussions. This work was supported by
the National Science Foundation (NSF) within the framework of the JETSCAPE
Collaboration under Award No. ACI-1550223. Additional partial support by the
U.S. Department of Energy (DOE), Office of Science, Office for Nuclear Physics
under Award No. DE-SC0004286 and within the framework of the BEST and JET
Collaborations is also acknowledged.

##  Appendix A Moments of the distribution function

We define the thermal, isotropic and anisotropic integrals that appear in this
paper. For each species $n$ let

$J_{kq,n}\equiv\int_{p}\frac{(u{\,\cdot\,}p)^{k-2q}({-}p{\,\cdot\,}\Delta{\,\cdot\,}p)^{q}}{(2q{+}1)!!}f_{\text{eq},n}\bar{f}_{\text{eq},n}\,.$
(84)

The thermal integrals over the local-equilibrium distribution are then given
by

$\displaystyle\mathcal{J}_{kq}$ $\displaystyle=\sum_{n}J_{kq,n}\,,$ (85a)
$\displaystyle\mathcal{N}_{kq}$ $\displaystyle=\sum_{n}b_{n}J_{kq,n}\,,$ (85b)
$\displaystyle\mathcal{M}_{kq}$ $\displaystyle=\sum_{n}b^{2}_{n}J_{kq,n}\,,$
(85c) $\displaystyle\mathcal{A}_{kq}$
$\displaystyle=\sum_{n}m_{n}^{2}J_{kq,n}\,,$ (85d)
$\displaystyle\mathcal{B}_{kq}$
$\displaystyle=\sum_{n}b_{n}m_{n}^{2}J_{kq,n}\,.$ (85e)

The isotropic integrals over the PTB distribution (39) without shear stress
modifications (i.e. for $\pi_{ij}=0$) are

$\mathcal{L}_{kq}=\sum_{n}\int_{p}\frac{(u{\,\cdot\,}p)^{k-2q}({-}p{\,\cdot\,}\Delta{\,\cdot\,}p)^{q}}{(2q{+}1)!!}f_{\lambda,n}\,,$
(86)

where

$f_{\lambda,n}=\frac{g_{n}(1{+}\lambda_{\Pi})^{-3}}{\exp\left[\dfrac{1}{T}\sqrt{m_{n}^{2}-\dfrac{p{\,\cdot\,}\Delta{\,\cdot\,}p}{(1{+}\lambda_{\Pi})^{2}}}\right]+\Theta_{n}}\,.$
(87)

In particular, the modified energy density and isotropic pressure are given by
$\mathcal{E}^{\prime}(\lambda_{\Pi},z_{\Pi},T)=z_{\Pi}\mathcal{L}_{20}(\lambda_{\Pi},T)$
and
$\mathcal{P}^{\prime}(\lambda_{\Pi},z_{\Pi},T)=z_{\Pi}\mathcal{L}_{21}(\lambda_{\Pi},T)$.

Finally, the anisotropic integrals are given by

$\begin{split}\mathcal{J}_{krqs}=&\sum_{n}\frac{1}{(2q)!!}\int_{p}(u{\,\cdot\,}p)^{k{-}r{-}2q}(-Z{\,\cdot\,}p)^{r}\\\
&\ \ \times({-}p{\,\cdot\,}\Xi{\,\cdot\,}p)^{q}(p\cdot\Omega\cdot
p)^{s/2}f_{a,n}\bar{f}_{a,n}\,.\end{split}$ (88)

##  Appendix B Conservation laws in first-order approximation

Here we derive the expression for the time derivatives (20) by making a first-
order approximation to the conversation laws (19). The conservation equations
for the net baryon number and energy up to first order in gradients are

$\dot{n}_{B}=-n_{B}\theta,\qquad\dot{\mathcal{E}}=-(\mathcal{E}{+}\mathcal{P}_{\text{eq}})\theta\,.$
(89)

Taking the time derivative of the kinetic definitions for the net baryon and
energy densities,

$\\!\\!n_{B}=\sum_{n}b_{n}\int_{p}(u{\,\cdot\,}p)f_{\text{eq},n},\quad\mathcal{E}=\sum_{n}\int_{p}(u{\,\cdot\,}p)^{2}f_{\text{eq},n},\\!\\!$
(90)

one obtains

$\dot{\alpha}_{B}=\mathcal{G}(T,\alpha_{B})\theta\,,\qquad\dot{T}=\mathcal{F}(T,\alpha_{B})\theta\,.$
(91)

The coefficients $\mathcal{G}$ and $\mathcal{F}$ appear in the RTA
Chapman–Enskog expansion (22):

$\displaystyle\mathcal{G}$
$\displaystyle=T\left(\frac{(\mathcal{E}{+}\mathcal{P}_{\text{eq}})\mathcal{N}_{20}-n_{B}\mathcal{J}_{30}}{\mathcal{J}_{30}\mathcal{M}_{10}-\mathcal{N}_{20}^{2}}\right)\,,$
(92a) $\displaystyle\mathcal{F}$
$\displaystyle=T^{2}\left(\frac{n_{B}\mathcal{N}_{20}-(\mathcal{E}{+}\mathcal{P}_{\text{eq}})\mathcal{M}_{10}}{\mathcal{J}_{30}\mathcal{M}_{10}-\mathcal{N}_{20}^{2}}\right)\,.$
(92b)

The evolution equation for the fluid velocity up to first-order in gradients
is given by

$\dot{u}^{\mu}=\frac{\nabla^{\mu}\mathcal{P}_{\text{eq}}}{\mathcal{E}{+}\mathcal{P}_{\text{eq}}}\,.$
(93)

After computing the spatial gradient of the equilibrium pressure
$\mathcal{P}_{\text{eq}}=\frac{1}{3}\sum_{n}\int_{p}(-p{\,\cdot\,}\Delta{\,\cdot\,}p)f_{\text{eq},n}$,
Eq. (93) can be rewritten as

$\dot{u}^{\mu}=\nabla^{\mu}{\ln}T+\frac{n_{B}T}{\mathcal{E}{+}\mathcal{P}_{\text{eq}}}\nabla^{\mu}\alpha_{B}\,,$
(94)

where we made use of the identities Jaiswal  _et al._ (2014)
$\mathcal{J}_{31}=(\mathcal{E}{+}\mathcal{P}_{\text{eq}})T$ and
$\mathcal{N}_{21}=n_{B}T$.

## References

  * Song  _et al._ (2011) H. Song, S. A. Bass, U. Heinz, T. Hirano, and C. Shen, Phys. Rev. Lett. 106, 192301 (2011), [Erratum: Phys. Rev. Lett. 109, 139904 (2012)], arXiv:1011.2783 [nucl-th] .
  * Schenke  _et al._ (2012) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C85, 024901 (2012), arXiv:1109.6289 [hep-ph] .
  * Gale  _et al._ (2013a) C. Gale, S. Jeon, B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev. Lett. 110, 012302 (2013a), arXiv:1209.6330 [nucl-th] .
  * Gale  _et al._ (2013b) C. Gale, S. Jeon, and B. Schenke, Int. J. Mod. Phys. A28, 1340011 (2013b), arXiv:1301.5893 [nucl-th] .
  * Shen (7 25) C. Shen,  _The standard model for relativistic heavy-ion collisions and electromagnetic tomography_ , Ph.D. thesis, Ohio State U. (2014-07-25).
  * Ryu  _et al._ (2015) S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * Bernhard  _et al._ (2016) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Bernhard (4 19) J. E. Bernhard,  _Bayesian parameter estimation for relativistic heavy-ion collisions_ , Ph.D. thesis, Duke University (2018-04-19), arXiv:1804.06469 [nucl-th] .
  * Heinz and Kolb (2002) U. Heinz and P. F. Kolb, Nucl. Phys. A 702, 269 (2002), arXiv:hep-ph/0111075 .
  * Romatschke and Romatschke (2007) P. Romatschke and U. Romatschke, Phys. Rev. Lett. 99, 172301 (2007), arXiv:0706.1522 [nucl-th] .
  * Schenke  _et al._ (2011) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Schenke  _et al._ (2010) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Denicol  _et al._ (2012) G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Phys. Rev. D85, 114047 (2012), [Erratum: Phys. Rev.D91,no.3,039902(2015)], arXiv:1202.4551 [nucl-th] .
  * Heinz and Snellings (2013) U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
  * Shen  _et al._ (2016) C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, Comput. Phys. Commun. 199, 61 (2016), arXiv:1409.8164 [nucl-th] .
  * Borsanyi  _et al._ (2010a) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo (Wuppertal-Budapest), JHEP 09, 073 (2010a), arXiv:1005.3508 [hep-lat] .
  * Borsanyi  _et al._ (2010b) S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, JHEP 11, 077 (2010b), arXiv:1007.2580 [hep-lat] .
  * Bazavov  _et al._ (2012) A. Bazavov _et al._ , Phys. Rev. D 85, 054503 (2012), arXiv:1111.1710 [hep-lat] .
  * Borsanyi  _et al._ (2014) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Phys. Lett. B 730, 99 (2014), arXiv:1309.5258 [hep-lat] .
  * Bazavov  _et al._ (2014) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Huovinen and Petersen (2012) P. Huovinen and H. Petersen, Eur. Phys. J. A48, 171 (2012), arXiv:1206.3371 [nucl-th] .
  * Cooper and Frye (1974) F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
  * Grad (1949) H. Grad, Communications on Pure and Applied Mathematics 2, 331 (1949).
  * Chapman  _et al._ (1990) S. Chapman, T. Cowling, D. Burnett, and C. Cercignani,  _The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases_ , Cambridge Mathematical Library (Cambridge University Press, 1990).
  * Anderson and Witting (1974) J. Anderson and H. Witting, Physica 74, 466 (1974).
  * Denicol  _et al._ (2009) G. S. Denicol, T. Kodama, T. Koide, and P. Mota, Phys. Rev. C80, 064901 (2009), arXiv:0903.3595 [hep-ph] .
  * Pratt and Torrieri (2010) S. Pratt and G. Torrieri, Phys. Rev. C82, 044901 (2010), arXiv:1003.0413 [nucl-th] .
  * Romatschke and Strickland (2003) P. Romatschke and M. Strickland, Phys. Rev. D68, 036004 (2003), arXiv:hep-ph/0304092 [hep-ph] .
  * Martinez and Strickland (2010) M. Martinez and M. Strickland, Nucl. Phys. A848, 183 (2010), arXiv:1007.0889 [nucl-th] .
  * Heinz  _et al._ (2016) U. Heinz, D. Bazow, G. Denicol, M. Martinez, M. Nopoush, J. Noronha, R. Ryblewski, and M. Strickland, Nucl. Part. Phys. Proc. 276-278, 193 (2016), arXiv:1509.05818 [nucl-th] .
  * Alqahtani  _et al._ (2017) M. Alqahtani, M. Nopoush, R. Ryblewski, and M. Strickland, Phys. Rev. C96, 044910 (2017), arXiv:1705.10191 [nucl-th] .
  * Almaalol  _et al._ (2019) D. Almaalol, M. Alqahtani, and M. Strickland, Phys. Rev. C99, 044902 (2019), arXiv:1807.04337 [nucl-th] .
  * Nopoush and Strickland (2019) M. Nopoush and M. Strickland, Phys. Rev. C 100, 014904 (2019), arXiv:1902.03303 [nucl-th] .
  * Monnai and Hirano (2009) A. Monnai and T. Hirano, Phys. Rev. C80, 054906 (2009), arXiv:0903.4436 [nucl-th] .
  * De Groot  _et al._ (1980) S. R. De Groot, W. A. Van Leeuwen, and C. G. Van Weert,  _Relativistic Kinetic Theory. Principles and Applications_ (North Holland / Elsevier, 1980).
  * Jaiswal  _et al._ (2014) A. Jaiswal, R. Ryblewski, and M. Strickland, Phys. Rev. C90, 044908 (2014), arXiv:1407.7231 [hep-ph] .
  * Everett  _et al._ (2021) D. Everett, C. Chattopadhyay, and U. Heinz, (2021), arXiv:2101.01130 [hep-ph] .
  * Bleicher  _et al._ (1999) M. Bleicher _et al._ , J. Phys. G25, 1859 (1999), arXiv:hep-ph/9909407 [hep-ph] .
  * Andronic  _et al._ (2018) A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, Nature 561, 321 (2018), arXiv:1710.09425 [nucl-th] .
  * Bazow  _et al._ (2014) D. Bazow, U. Heinz, and M. Strickland, Phys. Rev. C90, 054910 (2014), arXiv:1311.6720 [nucl-th] .
  * Molnar  _et al._ (2016) E. Molnar, H. Niemi, and D. H. Rischke, Phys. Rev. D93, 114025 (2016), arXiv:1602.00573 [nucl-th] .
  * McNelis  _et al._ (2018) M. McNelis, D. Bazow, and U. Heinz, Phys. Rev. C97, 054912 (2018), arXiv:1803.01810 [nucl-th] .
  * McNelis  _et al._ (2021a) M. McNelis, D. Bazow, and U. Heinz, (2021a), arXiv:2101.02827 [nucl-th] .
  * Moreland  _et al._ (2015) J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Everett  _et al._ (2020a) D. Everett _et al._ (JETSCAPE), (2020a), arXiv:2011.01430 [hep-ph] .
  * Everett  _et al._ (2020b) D. Everett _et al._ (JETSCAPE), (2020b), arXiv:2010.03928 [hep-ph] .
  * McNelis  _et al._ (2021b) M. McNelis, D. Everett, and U. Heinz, Comput. Phys. Commun. 258, 107604 (2021b), arXiv:1912.08271 [nucl-th] .
  * Denicol  _et al._ (2014) G. S. Denicol, W. Florkowski, R. Ryblewski, and M. Strickland, Phys. Rev. C90, 044905 (2014), arXiv:1407.4767 [hep-ph] .
  * Schnedermann  _et al._ (1993) E. Schnedermann, J. Sollfrank, and U. W. Heinz, Phys. Rev. C 48, 2462 (1993), arXiv:nucl-th/9307020 .
  * Weil  _et al._ (2016) J. Weil _et al._ , Phys. Rev. C94, 054905 (2016), arXiv:1606.06642 [nucl-th] .