# Far-off-equilibrium expansion trajectories in the QCD phase diagram

Chandrodoy Chattopadhyay  Department of Physics, North Carolina State
University, Raleigh, NC 27695, USA  Ulrich Heinz  Department of Physics, The
Ohio State University, Columbus, OH 43210, USA  Thomas Schäfer  Department of
Physics, North Carolina State University, Raleigh, NC 27695, USA

###### Abstract

We consider the hydrodynamic evolution of a quark-gluon gas with non-zero
quark masses and net baryon number in its phase diagram. For far-off-
equilibrium initial conditions the expansion trajectories appear to violate
simple rules based on the second law of thermodynamics that were previously
established for ideal or weakly dissipative fluids. For Bjorken flow we
present a detailed analysis within kinetic theory that provides a full
microscopic understanding of these macroscopic phenomena and establishes their
thermodynamic consistency. We point out that, for certain far-off-equilibrium
initial conditions, the well-known phenomenon of “viscous heating” turns into
“viscous cooling” where, driven by dissipative effects, the temperature
decreases faster than in adiabatic expansion.

##  I Introduction

One of the primary goals of relativistic heavy-ion collisions is to study the
phase diagram of Quantum Chromodynamics (QCD). In a simplified version which
ignores additional dimensions Stephanov:2006zvm  associated with strangeness
and isospin (which are conserved by the strong interactions), this phase
diagram reduces to a plane spanned by the temperature $T$ and the baryon
chemical potential $\mu_{B}$. Based on first principles lattice QCD
calculations at vanishing $\mu_{B}$ Aoki:2006we ; Bazavov:2009zn ;
Borsanyi:2010cj ; Bazavov:2011nk  and theoretical modeling at large $\mu_{B}$
Stephanov:2006zvm , it is conjectured that there exists a first-order phase
transition line in this plane which terminates at a critical point at finite
$(T_{\mathrm{cr}},\mu_{B,\mathrm{cr}})$ Stephanov:1998dy ; Stephanov:1999zu ,
turning into a continuous crossover at smaller $\mu_{B}$ values. The search
for the QCD critical point has been an area of intense research over the last
two decades Busza:2018rrf ; Bzdak:2019pkr .

A key feature common to equilibrated systems near a critical point is the
enhancement of thermodynamic fluctuations accompanied by diverging correlation
length of these fluctuations hohenberg1977theory . Accordingly, the QCD
critical point search program using heavy-ion collisions is driven by
exploration of final-state observables that are sensitive to the correlation
length. One such class of observables are higher-order cumulants of net-proton
fluctuations, which are predicted to exhibit non-monotonic behavior as the
freeze-out region of collision fireballs is systematically varied in the
$(T,\mu_{B})$ plane from low to high $\mu_{B}$ Stephanov:1998dy ;
Stephanov:1999zu ; Hatta:2003wn ; Stephanov:2008qz ; Stephanov:2011pb ;
Luo:2017faz . This strategy has been pursued in the two-stage Beam Energy Scan
(BES) program at RHIC, by varying the collision energy $\sqrt{s_{NN}}$ from
$200$ GeV down to $7.7$ GeV (for BES-I) and $3.0$ GeV (for BES-II)
Bzdak:2019pkr ; STAR:2017sal .

However, this simple concept is complicated along several fronts by the
realities of relativistic heavy-ion collisions Nahrgang:2011mg ;
Nahrgang:2011mv ; Mukherjee:2015swa ; Herold:2016uvv ; Stephanov:2017ghc ;
Hirano:2018diu ; Singh:2018dpk ; Nahrgang:2018afz ; Yin:2018ejt ; An:2021wof .
The main complication arises from the explosive expansion dynamics of the
fireballs formed in the collisions which drive the fireball medium out of
local thermal equilibrium. At the lower BES energies additional considerations
make thermalization even more elusive Gupta:2022phu . Equally important is the
observation Steinheimer:2007iy ; Du:2021zqz  (reviewed in Shen:2020gef ;
An:2021wof ) that different fluid cells inside the expanding medium follow
different expansion trajectories through the phase diagram such that the
measured final-state particles represent an average over contributions with
different chemical compositions at freeze-out. This has led to revived
interest in fireball expansion trajectories through the QCD phase diagram.

Such expansion trajectories were first studied over 35 years ago in the
context of ideal fluid dynamics where both entropy and baryon number are
conserved and the expansion therefore occurs at fixed entropy per baryon,
$S/A{\,=\,}s/n{\,=\,}$const. (where $S$ is the total entropy, $A$ is the
nuclear mass (baryon) number, $s=S/V$ is the entropy density, and $n=A/V$ is
the net baryon density) Heinz:1987ca ; Lee:1992hn ; Cho:1993mv ; Hung:1997du .
From these studies it is known that a thermalized gas of massless quarks,
antiquarks and gluons expands isentropically along lines of constant
$\mu_{B}/T$, ending at the origin $(T_{\infty}=0,\mu_{B\infty}=0)$, whereas
for a thermalized gas of hadron resonances the adiabatic expansion
trajectories end at $(T_{\infty}=0,\,\mu_{B\infty}=m_{N})$ (where $m_{N}$ is
the mass of the nucleon, the lightest baryon number carrying hadron). It was
also understood that, at constant temperature $T>0$, an increase in specific
entropy $s/n$ leads in local thermal equilibrium to a decrease of the chemical
potential $\mu_{B}$, irrespective of the mass of the constituents. It came
therefore to the present authors as a surprise when in a recent study
Dore:2020jye  of dissipative expansion trajectories, where viscous heating
leads to an increase of entropy with time, for some initial conditions in the
quark-gluon plasma phase the chemical potential $\mu_{B}$ did not initially
decrease, but rather grow with decreasing temperature. Clearly, along such
trajectories the equilibrium value of the specific entropy $s/n$ decreased
initially Dore:2022qyz . What causes this decrease, and how can this be
compatible with the second law of thermodynamics asserting that the total
entropy can never decrease (Boltzmann’s “H-theorem”)? The authors of Ref.
Dore:2022qyz  pointed to non-equilibrium entropy corrections as the likely
culprit for this phenomenon. In this paper we show that this conjecture was
correct and explain in detail all of the possible manifestations of
dissipative effects in the heavy-ion fireball expansion trajectories through
the QCD phase diagram.

We point out that the above phenomenon of decreasing equilibrium entropy is
associated with more rapid cooling along the expansion trajectory than
expected for adiabatic expansion. One usually associates dissipative effects
with the phenomenon of “viscous heating” — a form of internal friction that
causes the medium to cool down more slowly in dissipative evolution than in
adiabatic expansion. To distinguish the far-off-equilibrium dynamics
discovered numerically by Dore et al. Dore:2020jye ; Dore:2022qyz  from this
expected behavior we introduce the concept of “viscous cooling” where the
dissipative system cools more rapidly than the ideal fluid.

To allow for a largely analytic discussion we follow Ref. Dore:2020jye  and
assume 1-dimensional Bjorken expansion, but our insights are generic and
generalize to arbitrary expansion geometries that can only be studied
numerically. Moreover, Dore et al. employed a Lattice QCD based equation of
state that included the effects of criticality Parotto:2018pwx . Since a
realistic equation of state will not be crucial for our discussion, we
simplify the system further and consider a weakly interacting gas of massive
quarks and gluons at finite chemical potential. Such a gas does not have a
phase transition or a critical point, but it has non-vanishing transport
coefficients arising from the weak interactions among its constituents. We
calculate these transport coefficients from kinetic theory self-consistently
with our model assumptions.

Dividing the entropy (or any other macroscopic hydrodynamic quantity) into
equilibrium and non-equilibrium contributions requires a definition of the
temperature $T$ and chemical potential $\mu_{B}$ associated with the
equilibrium part. This matching of equilibrium parameters to the non-
equilibrium state is ambiguous. An important aspect of the work presented here
is a comparison between macroscopic and microscopic descriptions; for the
former we use second-order hydrodynamics, whereas for the latter we implement
the relativistic Boltzmann equation with a collision kernel in Relaxation Time
Approximation (“RTA Boltzmann”). To ensure conservation of energy-momentum and
baryon number by the RTA collision kernel Romatschke:2011qp ; Denicol:2012cn ;
Florkowski:2013lya ; Jaiswal:2013npa ; Romatschke:2017acs 111This works for a
momentum independent relaxation time in the RTA collision kernel as used here.
For a momentum-dependent relaxation time the RTA collision term must be
generalized to remain compatible with these conservation laws under Landau
matching Teaney:2013gca ; Rocha:2021zcw ; Rocha:2021lze ; Dash:2021ibx . we
employ the Landau matching conditions DeGroot:1980dk
$e=e_{\mathrm{eq}}(T,\mu_{B})$ and $n=n_{\mathrm{eq}}(T,\mu_{B})$, i.e. we
assign $T$ and $\mu_{B}$ such that the associated equilibrium values reproduce
the actual energy and net baryon densities of the fluid cell in its non-
equilibrium state.

This paper is organized as follows: In Sec. II we describe the non-conformal
hydrodynamic equations studied in this work. Since we study a multi-component
quantum gas at finite baryon chemical potential $\mu_{B}$, its transport
coefficients differ from those used in Refs. Dore:2020jye  for a single-
component Boltzmann gas without conserved charge Denicol:2014vaa ; these
transport coefficients are derived in App. A. Section III briefly summarizes
the characteristics of the isentropic expansion trajectories in ideal fluid
dynamics for conformal and non-conformal systems. This is followed by a
discussion of dissipative expansion trajectories for conformal systems in Sec.
IV and for non-conformal systems in Sec. V. In each case two different
subsections describe and analyze the hydrodynamic and kinetic theory solutions
for these trajectories; for the non-conformal case the two are compared in
App. C. We explain how, for certain types of far-off-equilibrium initial
conditions, the system can expand following trajectories along which the
equilibrium part of the entropy initially decreases. For the conformal case,
the splitting of the entropy current into equilibrium and non-equilibrium
parts is described in App. B. Our conclusions are presented in Sec. VI.

##  II Dissipative hydrodynamics with Bjorken flow

Bjorken expansion is defined by the flow velocity field
$v^{x}=v^{y}=0,v^{z}=z/t$ Bjorken:1982qr . Symmetries dictate that all
macroscopic quantities are independent of transverse position $(x,y)$ and
space-time rapidity $\eta_{s}\equiv\tanh^{-1}(z/t)$, and thus depend only on
Milne proper time $\tau\equiv\sqrt{t^{2}-z^{2}}$. Together with invariance
under longitudinal reflections ($\eta_{s}\to-\eta_{s}$), these symmetries
force the dissipative part of the conserved baryon current to vanish, and the
shear stress tensor reduces to a single independent component
$\pi(\tau)\equiv\pi^{\eta_{s}}_{\eta_{s}}(\tau)$. The hydrodynamic evolution
equations reduce to Dore:2020jye ; Jaiswal:2014isa

$\displaystyle\frac{de}{d\tau}$
$\displaystyle=-\frac{1}{\tau}\left(e+P+\Pi-\pi\right),$ (1)
$\displaystyle\frac{dn}{d\tau}$ $\displaystyle=-\frac{n}{\tau},$ (2)
$\displaystyle\frac{d\Pi}{d\tau}$
$\displaystyle=-\frac{\Pi}{\tau_{R}}-\frac{\beta_{\Pi}}{\tau}-\delta_{\Pi\Pi}\,\frac{\Pi}{\tau}+\lambda_{\Pi\pi}\,\frac{\pi}{\tau},$
(3) $\displaystyle\frac{d\pi}{d\tau}$
$\displaystyle=-\frac{\pi}{\tau_{R}}+\frac{4}{3}\,\frac{\beta_{\pi}}{\tau}-\left(\frac{1}{3}\,\tau_{\pi\pi}+\delta_{\pi\pi}\right)\frac{\pi}{\tau}+\frac{2}{3}\,\lambda_{\pi\Pi}\,\frac{\Pi}{\tau},$
(4)

where $e$ and $n$ are the energy and net quark density, $P$ is the thermal
equilibrium pressure, given by the equation of state $P(e,n)=P(T,\mu)$, and
$\Pi$ and $\pi$ are the bulk and shear viscous stresses. $\tau_{R}$ is the
microscopic relaxation time in an underlying kinetic description to be
discussed below, and
$\beta_{\pi},\,\beta_{\Pi},\,\delta_{\pi\pi},\,\delta_{\Pi\Pi},\,\tau_{\pi\pi},\,\lambda_{\pi\Pi},\,\lambda_{\Pi\pi}$
are associated transport coefficients, derived in Appendix A using the
Chapman-Enskog expansion Jaiswal:2014isa .

In Bjorken flow the system volume grows linearly with Milne time $\tau$. Eq.
(2) thus reflects conservation of the comoving net quark number per unit
transverse area, $n\tau$, owing to the lack of net quark diffusion out of a
fluid cell. Equation (1) shows that the same is not true for the energy
density: As long as the effective longitudinal pressure $P{+}\Pi{-}\pi$ is
positive, the comoving energy density per unit area, $e\tau$, decreases with
$\tau$ due to work done by the pressure which converts some of the thermal
energy into kinetic motion energy associated with the collective longitudinal
expansion.

The evolution equations (3,4) for the bulk and shear viscous stresses are
derived by iteratively solving (up to second-order in velocity gradients) the
Boltzmann equation in relaxation-time approximation, followed by coarse-
graining and partial resummation of the gradient terms Jaiswal:2014isa ;
Jaiswal:2013npa . We parametrize the relaxation time $\tau_{R}$ in terms of
the local temperature $T$ as $T\tau_{R}=5C$ where $C$ is a unitless
constant.222This is strictly justified only in the massless limit where the
temperature is the only energy scale in the problem but we checked that
corrections in powers of $m/T$ are small enough to not affect any of the
qualitative features of our new results.,333For a conformal system at
vanishing chemical potential, the relation $\tau_{R}\beta_{\pi}\equiv\eta$
(where $\eta$ is the shear viscosity), along with the parametrization
$\tau_{R}=5C/T$, implies $C=\eta/s$. However, for $\mu\neq 0$, the same pair
of relations implies $C=\eta/s\times(1+(\mu/T)(n/s))^{-1}$.

As we will be interested in obtaining the hydrodynamic trajectories in the
$(T,\mu)$-plane we shall directly solve Eqs. (1)–(4) in terms of the Lagrange
parameters, $x_{a}\equiv\\{T,\mu\\}$, instead of the associated densities
$X_{a}=\\{e,n\\}$. Making use of the Landau matching conditions (55) (see App.
A) we express $X_{a}$ in terms of $x_{a}$ using the equilibrium distributions

$\displaystyle f^{q}_{\mathrm{eq}}$ $\displaystyle=\frac{1}{\exp(\beta
E_{p}-\alpha)+1},\,\,f^{\bar{q}}_{\mathrm{eq}}=\frac{1}{\exp(\beta
E_{p}+\alpha)+1},$ $\displaystyle f^{g}_{\mathrm{eq}}$
$\displaystyle=\frac{1}{\exp(\beta E_{p})-1}.$ (5)

Here $\beta$ is the inverse temperature, $\alpha\equiv\beta\mu$, and
$E_{p}=\sqrt{\bm{p}^{2}+m^{2}}$ is the energy of a particle of momentum
$\bm{p}$ and mass $m$ in the fluid rest frame.

For massless quarks and gluons with degeneracy factors $g_{q}$ and $g_{g}$,
respectively, the desired relations are simply

$\displaystyle e_{\mathrm{eq}}(T,\mu)$
$\displaystyle=T^{4}\Bigg{[}\frac{(4g_{g}+7g_{q})\pi^{2}}{120}+\frac{g_{q}}{4}\left(\frac{\mu}{T}\right)^{2}+\frac{g_{q}}{8\pi^{2}}\left(\frac{\mu}{T}\right)^{4}\Bigg{]}$
$\displaystyle=3P_{\mathrm{eq}}(T,\mu),$ (6) $\displaystyle
n_{\mathrm{eq}}(T,\mu)$
$\displaystyle=T^{3}\left[\frac{g_{q}}{6}\left(\frac{\mu}{T}\right)+\frac{g_{q}}{6\pi^{2}}\left(\frac{\mu}{T}\right)^{3}\right],$
(7)

which we use to compute $dx_{a}=M_{a}^{\,b}\,dX_{b}$ and convert Eqs. (1,2)
into evolution equations for $T$ and $\mu$. We use $g_{q}=2\times N_{c}\times
N_{f}=12$ (with $N_{c}{\,=\,}3$ colors and $N_{f}{\,=\,}2$ flavors) for quarks
and antiquarks and $g_{g}=2\times(N_{c}^{2}-1)=16$ for gluons.

Giving the (anti-)quarks a mass $m\neq 0$ changes these relations to

$\displaystyle e_{\mathrm{eq}}(T,\mu)$
$\displaystyle=\frac{T^{4}}{2\pi^{2}}\,\int_{0}^{\infty}du\,u^{2}\,\Bigg{\\{}\sqrt{u^{2}+\frac{m^{2}}{T^{2}}}\,g_{q}\Bigg{[}\tilde{f}^{q}_{\mathrm{eq}}\left(u,\frac{m}{T},\frac{\mu}{T}\right)$
$\displaystyle+\tilde{f}^{\bar{q}}_{\mathrm{eq}}\left(u,\frac{m}{T},\frac{\mu}{T}\right)\Bigg{]}+u\,g_{g}\,\tilde{f}^{g}_{\mathrm{eq}}(u)\Bigg{\\}},$
(8) $\displaystyle P_{\mathrm{eq}}(T,\mu)$
$\displaystyle=\frac{T^{4}}{6\pi^{2}}\,\int_{0}^{\infty}du\,u^{3}\,\biggl{\\{}\frac{u}{\sqrt{u^{2}+\frac{m^{2}}{T^{2}}}}\,g_{q}\Bigg{[}\tilde{f}^{q}_{\mathrm{eq}}\left(u,\frac{m}{T},\frac{\mu}{T}\right)$
$\displaystyle+\tilde{f}^{\bar{q}}_{\mathrm{eq}}\left(u,\frac{m}{T},\frac{\mu}{T}\right)\Bigg{]}+\,g_{g}\,\tilde{f}^{g}_{\mathrm{eq}}(u)\biggr{\\}},$
(9) $\displaystyle n_{\mathrm{eq}}(T,\mu)$
$\displaystyle=\frac{T^{3}}{2\pi^{2}}\,\int_{0}^{\infty}du\,u^{2}\,g_{q}\Bigg{[}\tilde{f}^{q}_{\mathrm{eq}}\left(u,\frac{m}{T},\frac{\mu}{T}\right)$
$\displaystyle-\tilde{f}^{\bar{q}}_{\mathrm{eq}}\left(u,\frac{m}{T},\frac{\mu}{T}\right)\Bigg{]},$
(10)

where
$\tilde{f}^{i}_{\mathrm{eq}}(u,z,\alpha)=\bigl{[}\exp\bigl{(}\sqrt{u^{2}{+}z^{2}}-\alpha_{i}\bigr{)}-\theta_{i}\bigr{]}^{-1}$
with $\alpha_{q}=-\alpha_{\bar{q}}=\alpha,\ \alpha_{g}=0$ as well as
$\theta_{q}=\theta_{\bar{q}}=-1,\ \theta_{g}=1$.

##  III Expansion trajectories in ideal hydrodynamics

We first solve the hydrodynamic equations in the absence of dissipation, by
setting the shear and bulk viscous stresses to zero throughout the evolution:
$\pi=\Pi=0$. The ideal expansion trajectories are shown in Fig. 1, for the
massless (conformal) case in panel (a) and for a quark-gluon gas with massive
($m=1$ GeV) quarks and antiquarks in panel (b). Throughout this work, initial
conditions are imposed at initial Milne time $\tau_{0}=0.1$ fm/$c$.

Figure 1: Ideal hydrodynamic expansion trajectories in the $(T,\mu)$ plane for
a conformal ($m=0$, panel (a)) and a non-conformal ($m_{q}=m_{\bar{q}}=1$ GeV,
$m_{g}=0$, panel (b)) quark-gluon gas with non-zero conserved net baryon
charge, for identical sets of initial values $(T_{0},\mu_{0})$ in (a) and (b).
The corresponding conserved entropy/baryon charge values $s/n$ are indicated.

In panel (a) the initial temperature and chemical potentials $(T_{0},\mu_{0})$
are chosen to yield the integer $s/n$ ratios indicated.444Note that
$n=\frac{1}{3}n_{B}$ is the net quark density such that
$s/n=\textstyle{\frac{1}{3}}(S/A)$ where $S$ and $A$ are the total entropy and
net baryon number of the system. For the non-conformal case, the same initial
values $(T_{0},\mu_{0})$ correspond to the larger $s/n$ ratios listed in panel
(b). In ideal fluids both entropy and net quark number of a fluid element in
its local rest frame are conserved during evolution, and the decrease of the
associated comoving densities arises entirely from volume expansion. This
follows directly from the thermodynamic relations $de=T\,ds+\mu\,dn$ and
$s=(e+P-\mu n)/T$ which allow to re-cast Eqs. (1), (2) as

$\displaystyle\frac{ds}{d\tau}=-\frac{s}{\tau},\qquad\frac{dn}{d\tau}=-\frac{n}{\tau},$
(11)

with the straightforward solution $s/n=\mathrm{const}$. Since for Bjorken flow
Milne coordinates serve as a comoving coordinate system, $s/n$ remains
constant as a function of Milne time $\tau$ which can be used as a line
parameter along each of the trajectories shown in Fig. 1.

For the conformal case shown in panel (a), with massless constituents yielding
the EoS (6,7), lines of constant $s/n$ correspond to lines of constant
$\mu/T$. This follows from conformality: in a system devoid of any
dimensionful microscopic scales, $s$ and $n$ (which share the same units) must
be related by

$s(T,\mu)=g(T,\mu)\,n(T,\mu)$ (12)

where $g$ can depend on $T$ and $\mu$ only through the dimensionless ratio
$\alpha=\mu/T$. Thus $s/n=g(\alpha)=\mathrm{constant}$ implies constancy of
$\alpha$. As a result, all the trajectories in Fig. 1a are straight lines
passing through the origin.555For each trajectory, its mirror image with
respect to the temperature axis is itself a solution with the opposite sign of
$s/n$, reflecting the oddness of the function $g(\alpha)$.

The analogous expansion trajectories for a non-conformal quark-gluon gas are
shown in Fig. 1b. In this case, since the function $g=s/n$ now depends on two
dimensionless ratios $\mu/T$ and $m/T$, they are no longer straight lines. To
illustrate the effect we have chosen a rather large quark mass of 1 GeV such
that the quark mass effects are clearly visible at all temperatures shown.
While at high temperatures $T/m\sim 1$ the expansion trajectories still
exhibit approximately linear behavior as observed in the massless case, they
drastically change shape at lower temperatures $T/m\sim 0.2-0.3$ where they
veer sharply to the right. This is forced by the existence of a non-vanishing
net baryon charge which, due to Fermi statistics, requires $\mu\to m=1$ GeV in
the zero temperature limit $T\to 0$.666This also explains the massless case
where $\mu\to m=0$ as $T\to 0$. The rightward bend of the trajectories near
$T/m\sim 0.2-0.3$ reflects the transition from a high-temperature state where
quarks, antiquarks and gluons contribute democratically to the thermodynamic
quantities, to a low-temperature state where antiquark and gluon contributions
to $n$ are exponentially suppressed relative to those from quarks by combined
mass and fugacity effects.

The late-time slopes of the various trajectories can be obtained from
thermodynamic arguments: At low temperatures, quarks dominate over antiquarks
and gluons, and their rest mass dwarfs their kinetic energy, such that $P\ll
e$. As a result, work done by the pressure can be neglected and the energy
density falls off $\propto 1/\tau$ while the ratio between energy and entropy
density, $e/s$, stays approximately constant. Using this in the thermodynamic
relation $Ts=e+P-\mu n$ yields

$\displaystyle T\approx-\frac{1}{(s/n)}\,\mu+{\cal C}$ (13)

where ${\cal C}=e/s\approx\,\mathrm{constant}$. Thus, at late-times the
isentropic trajectories are straight lines with negative slopes whose
magnitude is inversely proportional to $s/n$, as borne out in Fig. 1b.

We note for later use that it follows from Fig. 1 that in ideal fluids at
fixed temperature $g=s/n$ is a monotonic function of $\mu$, i.e. $\mu$
decreases when the specific entropy $s/n$ increases.

##  IV Dissipative expansion trajectories for conformal systems with non-zero
baryon chemical potential

We now proceed to discuss the expansion trajactories for dissipative systems.
We begin in this section with conformal systems without bulk viscous pressure,
$\Pi=0$. Non-conformal systems with non-vanishing values for both the shear
and bulk viscous stresses will be studied in the following Section V. In both
sections, we will begin with a discussion of macroscopic hydrodynamic
phenomena, followed by an in-depth microscopic analysis using kinetic theory.

###  IV.1 Conformal Bjorken hydrodynamics

Figure 2: (a) Expansion trajectories of a conformal gas in the phase diagram
using dissipative hydrodynamics. Solid lines show adiabatic expansion for
comparison. (b) Scaled time evolution of the normalised shear stress tensor.

We turn on dissipation by letting the shear stress tensor evolve via Eq. (4),
with $\Pi$ set to zero. For illustration the constant $C$ appearing in
$\tau_{R}=5C/T$ is set to $C=10/4\pi$, and we consider two sets of initial
temperatures and chemical potentials at initial time $\tau_{0}=0.1$ fm/$c$:
$(T_{0},\mu_{0})=(0.8,0.5)$ GeV and $(0.7,0.83)$ GeV, shown as blue and green
curves in Fig. 2. For both sets we further consider two different values for
the initial normalised shear stress: $(\pi/4p)_{0}=0.2$ and
$(\pi/4p)_{0}=-0.3$. The resulting dissipative expansion trajectories are
shown as dotted and dashed lines in Fig. 2a, together with solid lines for
ideal expansion starting from the same point $(T_{0},\mu_{0})$. Fig. 2b shows
the evolution of the normalized shear stress $\pi/(e{+}P)$ as a function of
the scaled time $\tau/\tau_{R}$. As the initial temperatures of the green and
blue sets of curves are slightly different, their starting values of
$\tau/\tau_{R}=\tau T/(5C)$ differ correspondingly. The solid red line in Fig.
2b shows the Navier-Stokes solution $\pi_{NS}=4\eta/3\tau$ which can also be
written as $\pi_{NS}/(e{+}P)=4\tau_{R}/(15\tau)$. It is a universal function
of the scaled time and provides a late-time attractor for the evolution of
$\pi/(e{+}P)$ Heller:2015dha ; Romatschke:2017vte ; Romatschke:2017acs ;
Jaiswal:2019cju ; Kurkela:2019set .

Looking at panel (a) we notice immediately that, in spite of the production of
additional entropy by dissipation, only the dotted lines (where the initial
shear stress $\pi_{0}$ is positive) exhibit the well-known phenomenon of
viscous heating where the dissipative fluid cools more slowly than the ideal
one. The dashed lines, corresponding to negative initial shear stress, instead
deviate from the ideal expansion trajectories initially towards the right
(i.e. towards larger $\mu$ or lower $T$, depending on your point of view),
i.e. they initially exhibit viscous cooling where the dissipative fluid cools
more rapidly than the ideal one. Only at later times do the dashed lines cross
over to the left, towards smaller $\mu$.777Qualitatively similar features were
first observed by the authors of Dore:2020jye  (see also Dore:2022qyz ). In
their work, which also includes bulk viscous effects, some dissipative
expansion trajectories even started moving towards the right of the starting
point $(T_{0},\mu_{0})$, not just right of the ideal isentrope with
$(s/n)(\tau)=(s/n)_{0}$. We will return to this observation further below.

Figure 3: Evolution of the total specific entropy $s_{\mathrm{tot}}/n$
computed from Eq. (16), along trajectories of corresponding line style and
color shown in Fig. 2a. The brown curves are obtained using the equilibrium
definition.

At first sight this feature appears to contradict the argument from the
preceding subsection that points with larger $(s/n)(\tau)>(s/n)_{0}$ should be
located in the phase diagram to the left of the isentropic expansion
trajectory with $(s/n)_{\mathrm{eq}}=(s/n)_{0}$. When (incorrectly)
interpreted within the ideal thermodynamic framework used in that subsection,
the specific entropy $s/n$ initially seems to decrease along the dashed
trajectories! The flaw in that argument is that it does not recognize the fact
that the dissipative stresses describe deviations of the microscopic phase-
space distributions from local thermal equilibrium which manifest themselves
in non-equilibrium corrections to the entropy density. Initial conditions with
non-vanishing shear stress (or, for that matter, any other dissipative flows)
have different $(s/n)_{0}$ than an ideal fluid prepared with the same initial
temperature and chemical potential $(T_{0},\mu_{0})$. An initial condition
with a large negative non-equilibrium correction to the entropy density can
produce additional total entropy by dissipative heating (thus satisfying the
second law of thermodynamics) while, at the same time, reducing the
equilibrium entropy, by transferring negative entropy from the non-equilibrium
to the equilibrium part as the system thermalizes and moves closer to local
equilibrium.

A little manipulation of Eqs. (1,2) in the presence of dissipation leads (for
$\Pi=0$) to

$\frac{d(s_{\mathrm{eq}}\tau)}{d\tau}=\frac{\pi}{T}.$ (14)

Here $s_{\mathrm{eq}}$ is the equilibrium entropy density, related to $e,\,p,$
and $n$ via the fundamental relation which holds in thermal equilibrium.
Clearly, the equilibrium entropy per transverse area of a fluid cell can
decrease with time whenever the shear stress $\pi$ is negative (i.e. far from
its positive first-order Navier-Stokes limit for Bjorken flow). This is
precisely the case for the dashed curves in Fig. 2a. Since in these cases it
takes some time for $\pi$ to turn positive (see Fig. 2b), the growth with time
of the total specific entropy by viscous heating is reflected also in its
equilibrium contribution only at later times.

For a conformal gas of quarks and gluons at finite chemical potential, the
second-order out-of-equilibrium entropy four-current using kinetic theory
is888This is derived in Appendix B, along with the coefficients $c_{nn}$ and
$c_{n\pi}$.

$\displaystyle S^{\mu}$
$\displaystyle=s_{\mathrm{eq}}\,u^{\mu}-\alpha\,n^{\mu}-\frac{\beta}{4\beta_{\pi}}\,u^{\mu}\,\pi^{\alpha\beta}\,\pi_{\alpha\beta}+c_{nn}\,u^{\mu}\,n^{\alpha}\,n_{\alpha}$
$\displaystyle\quad+\,c_{n\pi}\,\pi^{\mu\alpha}\,n_{\alpha},$ (15)

where $\pi^{\mu\nu}$ is the shear stress tensor and $n^{\mu}$ is the net-quark
diffusion current. Due to the vanishing of baryon diffusion in Bjorken flow,
the second-order entropy density simplifies to

$s_{\mathrm{tot}}=s_{\mathrm{eq}}-\frac{3\beta}{8\beta_{\pi}}\pi^{2},$ (16)

where the second term is a (negative!) non-equilibrium correction. The entropy
density is largest in local thermal equilibrium and reduced by dissipative
corrections: $s_{\mathrm{tot}}<s_{\mathrm{eq}}$. The second law of
thermodynamics applies only to the total specific entropy $s_{\mathrm{tot}}/n$
but not to the equilibrium part $s_{\mathrm{eq}}/n$ in isolation. This is
illustrated in Fig. 3 where only $s_{\mathrm{tot}}/n$ is seen to always grow
monotonically with time while for the dashed lines (corresponding to initially
negative shear stress) the equilibrium part $s_{\mathrm{eq}}/n$ initially
decreases with time. As the system approaches local thermal equilibrium at
late times, dissipative effects become small and the total and equilibrium
specific entropies converge and saturate.

We next turn to solving kinetic theory for Bjorken flow in order to find out
whether, and to what extent, the features observed in hydrodynamics also
manifest themselves in the underlying microscopic theory from which CE
hydrodynamics is obtained.

###  IV.2 Kinetic theory for a system of massless quarks and gluons

Let us consider a system of massless quarks, anti-quarks, and gluons whose
evolution is governed by the Boltzmann equation with a relaxation-type
collision kernel, sharing a common relaxation time $\tau_{R}$:

$p^{\mu}\partial_{\mu}f^{i}(x,p)=-\frac{u\cdot
p}{\tau_{R}}\left(f^{i}-f^{i}_{\mathrm{eq}}\right).$ (17)

Their equilibrium distributions are given by Eqs. (II), setting $m=0$. Energy-
momentum and net-baryon number conservation by the RTA collision kernel are
ensured (see, e.g., Romatschke:2011qp ; Romatschke:2017acs ; Rocha:2021zcw ;
Rocha:2021lze ; Dash:2021ibx ) if the macroscopic parameters $\beta=1/T$,
$\alpha=\mu/T$, and fluid four-velocity $u^{\mu}$ are defined by the Landau
matching conditions:

$\displaystyle u_{\mu}u_{\nu}T^{\mu\nu}$ $\displaystyle=\int
dP\,(p{\,\cdot\,}u)^{2}\,\Bigl{[}g_{q}\left(f^{q}{+}f^{\bar{q}}\right)+g_{g}f^{g}\Bigr{]}$
$\displaystyle=e_{\mathrm{eq}}(T,\mu),$ (18) $\displaystyle u_{\mu}N^{\mu}$
$\displaystyle=\int
dP\,(p{\,\cdot\,}u)\,g_{q}\left(f^{q}{-}f^{\bar{q}}\right)=n_{\mathrm{eq}}(T,\mu),$
(19)

with the right hands sides given by Eqs. (6-7). The Lorentz invariant
integration measure is defined as $dP\equiv d^{3}p/[(2\pi)^{3}E_{p}]$, with
$E_{p}=p$ here.

For a system undergoing Bjorken expansion, Eq. (17) reduces in Milne
coordinates to the ordinary differential equation

$\displaystyle\frac{df^{i}}{d\tau}=-\frac{f^{i}-f^{i}_{\mathrm{eq}}}{\tau_{R}},$
(20)

where the superscript $i$ denotes the species considered. The formal solution
of the above equation is Baym:1984np ; Florkowski:2013lya

$\displaystyle f^{i}(\tau;p_{T},w)$
$\displaystyle=D(\tau,\tau_{0})f^{i}_{\mathrm{in}}(\tau;p_{T},w)$
$\displaystyle+\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}}{\tau_{R}(\tau^{\prime})}\,D(\tau,\tau^{\prime})\,f^{i}_{\mathrm{eq}}(\tau^{\prime};p_{T},w),$
(21)

where $w\equiv p_{\eta}$ is the longitudinal momentum in Milne coordinates and
the damping function $D$ is defined by

$D(\tau_{2},\tau_{1})=\exp\biggl{(}-\int_{\tau_{1}}^{\tau_{2}}\frac{d\tau^{\prime}}{\tau_{R}(\tau^{\prime})}\biggr{)}.$
(22)

For the initial momentum distributions we choose a Romatschke-Strickland type
parametrization Romatschke:2003ms :999We note that this parametrization does
not allow for a full exploration of all possible values for the bulk viscous
pressure $\Pi$ allowed in kinetic theory Jaiswal:2021uvv . Very large negative
values of $\Pi$ require an additional fugacity factor that allows for
oversaturation of at least one of the three particle species considered here.
However, for reasons discussed below, we are not primarily interested in this
work in large negative $\Pi$ values – the hydrodynamic expansion trajectories
develop their most striking features when the initial value for $\Pi$ is
positive.

$\displaystyle f_{\mathrm{in}}^{q}$
$\displaystyle=\frac{1}{\exp\left(\sqrt{p_{T}^{2}+(1+\xi_{0})w^{2}/\tau_{0}^{2}+m^{2}}/\Lambda_{0}-\frac{\nu_{0}}{\Lambda_{0}}\right)+1},$
(23) $\displaystyle f_{\mathrm{in}}^{\bar{q}}$
$\displaystyle=\frac{1}{\exp\left(\sqrt{p_{T}^{2}+(1+\xi_{0})w^{2}/\tau_{0}^{2}+m^{2}}/\Lambda_{0}+\frac{\nu_{0}}{\Lambda_{0}}\right)+1},$
(24) $\displaystyle f_{\mathrm{in}}^{g}$
$\displaystyle=\frac{1}{\exp\left(\sqrt{p_{T}^{2}+(1+\xi_{0})w^{2}/\tau_{0}^{2}}/\Lambda_{0}\right)-1},$
(25)

where in this subsection we set the quark mass to $m=0$. For simplicity we
have taken a common anisotropy parameter $\xi_{0}$ and momentum scale
$\Lambda_{0}$ for all three species, and the same parameter $\nu_{0}$ for
quarks and anti-quarks characterizing a non-zero initial net-quark density. To
solve for temperature and chemical potential we use the Landau matching
conditions Florkowski:2013lya

$\displaystyle e_{\mathrm{eq}}(T(\tau),\mu(\tau))=D(\tau,\tau_{0}){\cal
E}(\tau)$ (26)
$\displaystyle\qquad\quad+\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}}{\tau_{R}(\tau^{\prime})}\,D(\tau,\tau^{\prime})\,{\cal
H}_{e}\Bigl{(}\frac{\tau^{\prime}}{\tau}\Bigr{)}\,e_{\mathrm{eq}}\bigl{(}T(\tau^{\prime}),\mu(\tau^{\prime})\bigr{)},$
$\displaystyle
n_{\mathrm{eq}}(T(\tau),\mu(\tau))=\frac{\tau_{0}}{\tau}\,n_{\mathrm{in}}.$
(27)

Noting the absence of net-quark diffusion in Bjorken flow we directly used the
non-diffusive solution for the net-quark density in the number matching
condition. The ${\cal H}_{e}$ function is defined as

$\displaystyle{\cal
H}_{e}(x)\equiv\frac{1}{2}\left(x^{2}+\frac{\tanh^{-1}\left(\sqrt{1-\frac{1}{x^{2}}}\right)}{\sqrt{1-\frac{1}{x^{2}}}}\right).$
(28)

The time-dependent function ${\cal E}$ arises from the momentum integral of
the initial distributions:

$\displaystyle{\cal E}(\tau)={\cal
H}_{e}\left(\frac{\tau_{0}}{\tau\sqrt{1+\xi_{0}}}\right)\,e_{\mathrm{eq}}(\Lambda_{0},\nu_{0}).$
(29)

The initial net-quark density is given by
$n_{\mathrm{in}}=n_{\mathrm{eq}}(\Lambda_{0},\nu_{0})/\sqrt{1+\xi_{0}}$.

The integral equations (26,27) are solved iteratively for the temperature and
chemical potential by numerical quadrature. After finding the solutions
$(T(\tau),\mu(\tau))$ we use them to compute the shear stress tensor component
$\pi\equiv\pi^{\eta}_{\eta}$:

$\displaystyle\pi(\tau)=D(\tau,\tau_{0})\,{\cal C}(\tau)$
$\displaystyle+\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}}{\tau_{R}(\tau^{\prime})}\,D(\tau,\tau^{\prime})\,{\cal
H}_{\pi}\Bigl{(}\frac{\tau^{\prime}}{\tau}\Bigr{)}\,e_{\mathrm{eq}}\bigl{(}T(\tau^{\prime}),\mu(\tau^{\prime})\bigr{)},$
(30)

where we defined

$\displaystyle{\cal H}_{\pi}(x)$
$\displaystyle\equiv\frac{1}{6\left(1-x^{2}\right)}\Bigg{[}x^{2}\,\left(1+2\,x^{2}\right)$
$\displaystyle+\left(1-4x^{2}\right)\frac{\tanh^{-1}\left(\sqrt{1-\frac{1}{x^{2}}}\right)}{\sqrt{1-\frac{1}{x^{2}}}}\Bigg{]}$
(31)

as well as

$\displaystyle{\cal C}(\tau)={\cal
H}_{\pi}\left(\frac{\tau_{0}}{\tau\sqrt{1+\xi_{0}}}\right)e_{\mathrm{eq}}(\Lambda_{0},\nu_{0}).$
(32)

We are particularly interested in computing the evolution of the out-of-
equilibrium entropy density in kinetic theory. For this, we start from the
following definition of the entropy current of a gas of quarks, anti-quarks,
and gluons:

$\displaystyle S^{\mu}(x)=-\sum_{i=1}^{3}g_{i}\int
dP\,p^{\mu}\,\phi_{i}[f^{i}].$ (33)

Here $i=1,2,3$ labels quarks, antiquarks, and gluons, respectively. The
function $\phi_{i}[f^{i}]$ is defined by

$\displaystyle\phi_{i}[f^{i}]=f^{i}\ln\left(f^{i}\right)-\frac{1+\theta_{i}\,f^{i}}{\theta_{i}}\ln\left(1+\theta_{i}\,f^{i}\right).$
(34)

It distinguishes between statistics of the constituents via the parameter
$\theta_{i}$: $\theta_{1}=\theta_{2}=-1$ (Fermi-Dirac) and $\theta_{3}=1$
(Bose-Einstein). For Bjorken flow, the entropy density $s_{\mathrm{tot}}\equiv
u_{\mu}S^{\mu}$ is given by

$\displaystyle
s_{\mathrm{tot}}(\tau)=-\sum_{i=1}^{3}g_{i}\int\frac{dp_{T}\,dw\,p_{T}}{4\pi^{2}\tau}\,\phi_{i}[f^{i}(\tau;p_{T},w)]$
(35)

where $f^{i}(\tau;p_{T},w)$ is the solution given in Eq. (IV.2).

Figure 4: Comparison of the expansion trajectories (panel (a)) and the
evolution of the inverse Reynolds number (panel (b)) between exact kinetic
theory (magenta and black curves) and their second-order conformal
hydrodynamic approximation (blue and green curves).

In Fig. 4 we compare, for identical initial conditions, the dissipative
expansion trajectories through the phase diagram (panel (a)) and the evolution
of the shear inverse Reynolds number $\pi/(e{+}P)$ (panel (b)) obtained from
the exact solution (26,27) of conformal kinetic theory (magenta and black
curves) with those computed with second-order conformal hydrodynamics (blue
and green curves, same as those shown in Fig. 2). This comparison tests the
accuracy of the hydrodynamic approximation to the underlying microscopic
dynamics described by the RTA Boltzmann equation. In particular, the
phenomenon of viscous cooling exhibited by the dashed-line trajectories in
panel (a) is seen to be a robust feature of the exact solution to the
underlying microscopic kinetic theory and not an artefact of the macroscopic
hydrodynamic approximation. The corresponding initial parameters
$(\Lambda_{0},\nu_{0},\xi_{0})$ in the initial distribution functions
(23)–(25) are listed in Table 1.

Color |  $\Lambda_{0}$ (GeV) |  $\nu_{0}$ (GeV) | $\xi_{0}$
---|---|---|---
Magenta dotted | $5.671$ | $4.142$ | $9.258$
Magenta dashed | $2.651$ | $1.997$ | $-0.885$
Black dotted | $4.856$ | $6.892$ | $9.258$
Black dashed | $2.257$ | $3.326$ | $-0.885$

Table 1: Initial parameters $(\Lambda_{0},\nu_{0},\xi_{0})$ used to generate
the kinetic theory trajectories in Fig. 4. Due to conformal invariance,
$\bar{\pi}_{0}$ depends only on $\xi_{0}$. Hence, the $\xi_{0}$ values for the
dotted ($\bar{\pi}_{0}=0.2$) and dashed ($\bar{\pi}_{0}=-0.3$) pairs of curves
are each identical.

In panel (a) the magenta and black curves for kinetic theory lie slightly to
the right of the corresponding blue and green hydrodynamic trajectories,
suggesting that their specific entropies are somewhat smaller in kinetic
theory than in the hydrodynamic approximation. This is further substantiated
in Fig. 5 below. Consistent with this observation, the normalised shear stress
shown in panel (b) is somewhat smaller in kinetic theory than for
hydrodynamics which is seen to over-predict the deviation from equilibrium,
especially if the initial shear stress is large and positive (dotted curves).
Up to these differences, hydrodynamics shows excellent agreement with kinetic
theory for conformal Bjorken evolution.

Figure 5: (a) Comparison of the out-of-equilibrium specific entropy evolution
in conformal RTA Boltzmann theory and its hydrodynamic approximation. Line
styles and colors are identical to Fig. 4. (b) Comparison of the kinetic
evolution of $s_{\mathrm{tot}}/n$ with its equilibrium part
$s_{\mathrm{eq}}/n$ (red lines) along the magenta and black trajectories shown
in Fig. 4a.

As noted earlier in a different context Chattopadhyay:2018apf , differences
between hydrodynamics and kinetic theory become more readily apparent in the
evolution of the out-of-equilibrium specific entropy. Its time evolution is
shown in Fig. 5a. Similar to Fig. 4, the black and magenta curves in Fig. 5a
denote kinetic theory trajectories, whereas the green and blue curves
correspond to hydrodynamics. The hydrodynamic solutions are identical to those
shown in Fig. 3. All kinetic theory solutions lie below the hydrodynamic ones,
consistent with the observations made in Figs. 4a and 4b. For example, for
$\bar{\pi}_{0}=0.2$ (dotted lines) the initial value of the total specific
entropy, $s_{\mathrm{tot}}/n$, computed in kinetic theory is about $7\%$
percent less than the one obtained within hydrodynamics. For the dotted black
and magenta curves, kinetic theory yields
$(s_{\mathrm{tot}}/n)_{0}=s_{\mathrm{in}}(\Lambda_{0},\nu_{0})/n_{\mathrm{in}}(\Lambda_{0},\nu_{0})\approx
5.92$ and $11.23$, respectively, whereas the dotted green and blue
hydrodynamic curves correspond to $(s_{\mathrm{tot}}/n)_{0}\approx 6.39$ and
$11.98$. This initial difference between dotted kinetic theory and dotted
hydrodynamic results persists throughout the evolution, resulting in an error
of approximately 10 percent at late times.

For the other initial condition, $\bar{\pi}_{0}=-0.3$ (dashed lines), the
kinetic theory and hydrodynamic trajectories start from approximately the same
initial value for $s_{\mathrm{tot}}/n$. However, the hydrodynamic trajectories
are characterised by larger entropy production and therefore differ at late
times from the microscopic theory prediction by $\sim 7\%$. It is fascinating
to note the sensitivity of entropy production to the initial conditions of the
fluid: in spite of starting from nearly equal initial values of
$s_{\mathrm{tot}}/n$, the dotted and dashed kinetic theory solutions
(corresponding to $\bar{\pi}_{0}=0.2$ and $\bar{\pi}_{0}=-0.3$, respectively)
differ appreciably in the late time values of $s_{\mathrm{tot}}/n$. For the
dashed black and magenta trajectories in Fig. 4a, the late-time $s/n$ is
approximately twice its initial value, whereas for the dotted trajectories the
viscous heating factor is more than 3 times. The reason for this is the
following: the total amount of specific entropy $\Delta(s/n)$ that can be
generated during the fluid’s equilibration is

$\displaystyle\Delta\left(\frac{s}{n}\right)=-\left(\frac{\delta
s_{\mathrm{neq}}}{n}\right)_{0}+\frac{1}{(n\tau)_{0}}\int_{\tau_{0}}^{\infty}\frac{d\tau}{T(\tau)}\,\pi(\tau),$
(36)

where $\delta s_{\mathrm{neq}}=s_{\mathrm{tot}}-s_{\mathrm{eq}}$. Although
both dotted and dashed kinetic theory trajectories have nearly equal $(\delta
s_{\mathrm{neq}}/n)_{0}$ initially, the dotted lines are characterised by
positive shear throughout, leading to a larger value of the integral in Eq.
(36) than for the dashed ones where the integrand has both negative and
positive contributions.

In Fig. 5b we compare the evolution of the equilibrium specific entropy
($s_{\mathrm{eq}}/n$, red curves) with the out-of-equilibrium
$s_{\mathrm{tot}}/n$ (magenta and black), both computed using kinetic theory.
The magenta and black curves are the same as those in panel (a). The observed
features are qualitatively similar to Fig. 3 obtained within hydrodynamics.
For initialisation $\bar{\pi}_{0}=0.2$ (dotted lines), $s_{\mathrm{eq}}/n$
increases monotonically — consistent with the dotted magenta and green phase
trajectories in Fig. 4a lying to the left of the isentropic straight line
trajectories shown in Fig. 2. In contrast, the dashed red curves for
$\bar{\pi}_{0}=-0.3$ show $s_{\mathrm{eq}}/n$ evolving non-monotonically: the
equilibrium part of the specific entropy initially decreases before starting
to increase after $\tau/\tau_{R}\approx 0.3$. This reflects a transfer of
large negative specific entropy from the non-equilibrium to the equilibrium
sector at early times. During this stage the local momentum distribution forms
a prolate ellipsoid whose deformation is continuously decreasing via free-
streaming which red-shifts the longitudinal momenta. Once the momentum
distribution reaches local isotropy, the sign of the shear stress $\pi$ flips
and, according to Eq. (14), so does the direction of entropy flow between the
non-equilibrium and equilibrium sectors. In Figure 4a the initial decrease of
$s_{\mathrm{eq}}/n$ manifests itself by causing the dashed magenta and black
curves to first move to the right of the straight-line ideal expansion
trajectory before veering left at late times, eventually crossing the ideal
line (see Fig. 2) due to the overall production of entropy by viscous heating
and the approach to local equilibrium at late times.

This transfer of entropy from the non-equilibrium to the equilibrium sector is
a novel effect which to our knowledge has not been previously described. The
second law of thermodynamics demands (i) $s/n\leq s_{\mathrm{eq}}/n$ and (ii)
$d(s/n)/d\tau\geq 0$, where $s=s_{\mathrm{tot}}$ stands for the total specific
entropy (i.e. the sum of its equilibrium contribution
$\bigl{(}s_{\mathrm{eq}}/n\bigr{)}(T,\mu)$ and its non-equilibrium correction
which depends on the dissipative flows). As time evolves, entropy can be
transferred from the (negative) non-equilibrium correction to the equilibrium
part. Therefore $d(s_{\mathrm{eq}}/n)/d\tau$ does not need to be positive
definite.

Comparison of the dashed and dotted curves in Fig. 4b shows that the rate
(sign and magnitude) at which entropy flows from the non-equilibrium to the
equilibrium sector depends on the values of the dissipative fluxes. Within the
hydrodynamic framework this is expressed by Eq. (14). However, while the
dashed lines in Fig. 4a move to the right of the ideal expansion trajectories,
they do not move towards larger chemical potential as has been observed for
some expansion trajectories shown in Refs. Dore:2020jye ; Dore:2022qyz . In
the following subsection we will therefore discuss the criteria that must be
satisfied for obtaining dissipative expansion trajectories that start out by
moving towards larger $\mu$. We will see that this very hard to achieve in
conformal systems but much easier in non-conformal systems. This will set the
stage for discussion of non-conformal expansion trajectories in Section V.

###  IV.3 Criteria for generating trajectories with increasing chemical
potential

Let us explore the necessary conditions required to generate Bjorken expansion
trajectories where the chemical potential increases with Milne time. The
discussion will be general, including both shear and bulk viscous stresses. We
start with the thermodynamic relations among differentials of the Lagrange
parameters $x_{a}=\\{T,\mu\\}$ and the densities $X_{a}=\\{e,n\\}$:

$\begin{pmatrix}dT\\\
d\mu\end{pmatrix}=\frac{1}{\mathcal{A}}\begin{pmatrix}a_{4}&-a_{2}\\\
-a_{3}&a_{1}\end{pmatrix}\begin{pmatrix}de\\\ dn\end{pmatrix},$ (37)

with the determinant $\mathcal{A}\equiv(a_{1}a_{4}-a_{3}a_{2})$ and the
positive definite thermodynamic response functions

$\displaystyle a_{1}$ $\displaystyle=\left(\frac{\partial e}{\partial
T}\right)_{\mu},\ a_{2}=\left(\frac{\partial e}{\partial\mu}\right)_{T},$
$\displaystyle a_{3}$ $\displaystyle=\left(\frac{\partial n}{\partial
T}\right)_{\mu},\ a_{4}=\left(\frac{\partial n}{\partial\mu}\right)_{T}.$ (38)

The determinant $\mathcal{A}$ is also positive since a thermodynamic
transformation $(de{\,>\,}0,\,dn{\,=\,}0)$ must lead to
$dT{\,>\,}0$.101010Similarly, a transformation where $de=0$ and $dn>0$ must
yield $d\mu>0$. We spot-checked both expectations numerically. Thus, the
criterion for an increase in chemical potential, $d\mu/d\tau\geq 0$, is

$\displaystyle-a_{3}\,\frac{de}{d\tau}+a_{1}\,\frac{dn}{d\tau}\geq 0\
\implies\ a_{3}\frac{e+P-\phi}{\tau}-a_{1}\frac{n}{\tau}\geq 0,$ (39)

where we used the equations of motion (1, 2) and introduced the shorthand
$\phi\equiv\pi-\Pi$. Accordingly, the amount of normalised dissipative stress,
$\bar{\phi}\equiv\phi/(e{+}P)$, required for generating $d\mu>0$ trajectories
is

$\displaystyle\bar{\phi}\leq
1-\frac{a_{1}}{a_{3}}\,\frac{n}{(e{+}P)}\equiv\psi.$ (40)

This obviously prefers negative shear stress values, $\pi<0$, and positive
bulk viscous pressure, $\Pi>0$. Both are “unnatural” for Bjorken expansion
where the Navier-Stokes values for the dissipative fluxes have the opposite
signs.

Figure 6: The solid blue line shows the upper bound on the normalised shear
stress that ensures $d\mu/d\tau>0$ in conformal Bjorken dynamics. The dashed
black line, $\psi_{c}=-0.5$, represents the lower bound on $\bar{\pi}$ allowed
by kinetic theory.

Let us first explore the implications of this criterion for the simpler case
of a conformal system $(\Pi{\,=\,}0)$ of massless quarks, anti-quarks, and
gluons at non-zero net baryon density.111111Note that the upper limit of
$\bar{\phi}$ depends on the equation of state of the system. Accordingly, the
bound changes somewhat when replacing quantum statistics (which we use here)
by classical Boltzmann statistics, or changing the degeneracy factors for
(anti-)quarks or gluons. In this case the response functions $a_{1}$ and
$a_{3}$ simplify to

$\displaystyle a_{1}=3\,s_{\mathrm{eq}},\quad a_{3}=\frac{g_{q}}{3}\,T\,\mu,$
(41)

such that the condition (40) yields

$\displaystyle\bar{\pi}\leq
1-3\,\Bigl{(}\frac{1}{2}+\frac{\alpha^{2}}{2\pi^{2}}\Bigr{)}\Big{/}\Bigl{(}1+\alpha\frac{n}{s_{\mathrm{eq}}}\Bigr{)}\equiv\psi_{c}(\alpha).$
(42)

A plot of $\psi_{c}(\alpha)$ is shown in Fig. 6. For the trajectory to move
towards the right in the phase diagram, $d\mu>0$, $\bar{\pi}$ must lie below
the blue curve. In kinetic theory there is a lower bound on the normalised
shear, $\bar{\pi}\geq-0.5$ (arising from the condition of non-negative
effective transverse pressure, $P_{T}=P+\pi/2\geq 0$), indicated by the black
dashed line. In conformal kinetic theory the condition $d\mu>0$ thus restricts
$\bar{\pi}$ to a very thin sliver of parameter space bounded by the solid blue
and dashed black lines which covers only a small fraction of the overall
parameter space $\bar{\pi}\in[-0.5,0.25]$. Moreover, even if the system is
initialized within this sliver at early times, it quickly moves out of it as
the normalised shear stress rapidly approaches the free-streaming attractor
$\bar{\pi}\approx 0.25$. Thus, in conformal kinetic theory it is nearly
impossible to find trajectories characterised by substantial increase in
chemical potential.

Figure 7: The upper bound $\psi_{\mathrm{nc}}(\alpha,z)$ in non-conformal
systems, as (a) a function of $\alpha=\mu/T$ for fixed values of $z=m/T$ and
(b) a function of $z$ for fixed values of $\alpha$. For
$\bar{\phi}<\psi_{\mathrm{nc}}$ expansion trajectories are characterised by
$d\mu/d\tau>0$.

However, the introduction of masses for the quarks changes the equation of
state of the system, and accordingly all thermodynamic quantities determining
the upper bound of $\bar{\phi}$ on the r.h.s. of Eq. (40) change, too. As
$\bar{\phi}$ is dimensionless, the bound can depend only on the dimensionless
variables $\alpha=\mu/T$ and $z\equiv m/T$. A plot of
$\psi_{\mathrm{nc}}(\alpha,z)$ (defined by the r.h.s. of the inequality (40))
is shown in Fig. 7, where in panel (a) $\alpha$ is varied while keeping $z$
constant, and vice-versa in panel (b). Clearly, with increasing quark mass the
allowed domain of $\bar{\phi}$ that allows for trajectories with
$d\mu/d\tau>0$ significantly increases. For instance, for typical initial
values $\alpha_{0}\sim 1$ and $m/T_{0}\sim 1$, one needs
$\bar{\phi}_{0}\lesssim-0.36$ which is much less restrictive than the
requirement $\bar{\phi}_{0}\lesssim-0.49$ for the same $\alpha_{0}$ in the
conformal case.

##  V Dissipative expansion trajectories for non-conformal systems with a
conserved charge

###  V.1 Non-conformal Bjorken hydrodynamics

Figure 8: Expansion trajectories from second-order non-conformal hydrodynamics
for a massive quark-gluon gas with non-zero net quark number, in the $T{-}\mu$
plane (a) and $e{-}n$ plane (b), respectively. The green dash-dotted lines
show the corresponding ideal expansion trajectories from Fig. 1b (for
$s/n=14.8$ and 8.2, respectively) for comparison. The other lines are
explained in the text.

We now relax the restriction to conformal symmetry and restore the bulk
viscous pressure $\Pi$ in the second-order viscous hydrodynamic equations
(1-4) and solve them with initial conditions imposed at $\tau_{0}=0.1$ fm/$c$
using the transport coefficients obtained in Appendix A for a massive quark-
gluon gas at $\mu\neq 0$. The resulting expansion trajectories in the
temperature-chemical potential plane (panel a) and the energy-(quark)density
plane (panel b) are shown in Fig. 8. We consider two sets of initial
conditions taken from Fig. 1, with $(T_{0},\mu_{0})\approx(0.8,0.5)$ GeV (set
I) and $\approx(0.7,0.83)$ GeV (set II), respectively.121212These correspond
to initial energy and net-quark densities of
$(e_{0},n_{0})\approx(667\,\mathrm{GeV/fm}^{3},70/\mathrm{fm}^{3})$ (set I)
and $\approx(476\,\mathrm{GeV/fm}^{3},94/\mathrm{fm}^{3})$ (set II). The
corresponding ideal expansion trajectories from Fig. 1b are shown as green
dash-dotted lines for orientation. In both panels the blue and black curves
show dissipative expansion trajectories for two choices of the coupling
strength, parametrized by the relaxation time constant $C$ as given in the
legend, while solid and dashed lines distinguish between different initial
conditions for the dissipative flows: dashed lines use local equilibrium
initial conditions $\bar{\pi}_{0}=\bar{\Pi}_{0}=0$ whereas solid lines
correspond to $\bar{\pi}_{0}\approx-0.437$, $\bar{\Pi}_{0}\approx 0.012$ for
set I and $\bar{\pi}_{0}\approx-0.439$, $\bar{\Pi}_{0}\approx 0.010$ for set
II, such that for both sets
$\bar{\phi}_{0}\equiv\bar{\pi}_{0}-\bar{\Pi}_{0}=-0.45$.131313The specific
choice of these somewhat arbitrary-looking values is motivated by the
comparison with kinetic theory in Sec. V.2 below. As explained in footnote 9,
the three-parameter Romatschke-Strickland distributions (23-25) limit the
range of initial dissipative flows that can be generated; in particular, this
parametrization takes away the freedom of choosing $\bar{\pi}_{0}$ and
$\bar{\Pi}_{0}$ independently once their difference $\bar{\phi}_{0}$ is
specified.

Similar to the conformal case shown in Fig. 2, the dashed and solid lines in
Fig. 8 form two qualitatively different classes of dissipative expansion
trajectories. The differences between them are much more pronounced in the
$T{-}\mu$ planel (panel (a)) than in the $e{-}n$ plane (panel (b)): the dashed
lines, corresponding to equilibrium initial conditions, first graze along the
ideal trajectories before viscous heating drives them left towards larger
equilibrium specific entropies. Since larger $C$ values (longer relaxation
times) describe more weakly coupled systems with larger shear and bulk
viscosities, they lead to stronger viscous entropy production. The solid
lines, on the other hand, which correspond to large negative initial values
for the (normalized) shear and bulk stress combination
$\bar{\phi}\equiv\bar{\pi}-\bar{\Pi}$, exhibit viscous cooling: in panel (a)
they move initially right, i.e. towards larger chemical potentials, as
anticipated in Sec. IV.3, and cool more rapidly than the ideal fluid. Again,
this phenomenon lasts longer for more viscous systems described by larger $C$
values.141414In the $e{-}n$ plane (panel (b)) the differences arising from
different coupling strengths $C$ are almost unrecognizable: the blue and black
dashed trajectories are indistinguishable, and also their solid analogues
(corresponding to off-equilibrium initial conditions) are very close to each
other.  The decrease in equilibrium specific entropy along these curves is due
to their initialisation with large negative $\phi$, as dictated by the non-
conformal version of Eq. (14),

$\displaystyle\frac{d(s_{\mathrm{eq}}\tau)}{d\tau}=\frac{\phi}{T}=\frac{\pi-\Pi}{T}\,.$
(43)

As described in Secs. IV.1 and IV.2, it reflects the transfer of negative
viscous entropy between the non-equilibrium and equilibrium contributions to
the total specific entropy.

Figure 9: Scaled time evolution of (a) shear and (b) bulk inverse Reynolds
numbers obtained from second-order hydrodynamics, corresponding to the
expansion trajectories (of same color and line-style) in set II of Fig. 8. The
dotted curves show the Navier-Stokes limit. (Note that, for fixed
$\tau_{0}=0.1$ fm/$c$, different $C$ values result in different scaled initial
times $\tau_{0}/\tau_{R}.$)

As the system evolves one expects the shear and bulk viscous pressures to
approach their Navier-Stokes limits (which is positive for shear stress and
negative for the bulk viscous pressure). This is illustrated in Fig.
9.151515To avoid clutter we show the evolution of $\bar{\pi}$ and $\bar{\Pi}$
only for set II. The results for set I are qualitatively similar. When this
happens, $\phi$ turns positive and $s_{\mathrm{eq}}/n$ increases again,
forcing the solid curves in Fig. 8 to turn left, towards regions of smaller
chemical potential. As the overall result of dissipation is an increase of
total entropy, the solid lines eventually have to cross the ideal expansion
trajectories and end up to their left where $s_{\mathrm{eq}}/n>(s/n)_{0}$.

As the system approaches local thermal equilibrium $s_{\mathrm{eq}}/n$
saturates, for different initial conditions generally at different values, and
the expansion trajectories settle on isentropic trajectories corresponding to
these asymptotic $s/n$ values. Therefore, once $T/m$ drops below about 0.3,
all trajectories again turn to the right, due to the nonzero mass of the
carriers of the conserved baryon charge as explained in Sec. III. A deeper
understanding of this on a microscopic level will be obtained in the next
subsection where the same phenomena are seen to arise in non-conformal kinetic
theory.

We close this subsection with some additional discussion of Fig. 9. Remember
that the dashed trajectories start with equilibrium initial conditions whereas
the solid lines correspond to large negative initial shear combined with a
much smaller initial bulk viscous pressure. Panel (a) shows that for both
types of initial conditions $\bar{\pi}$ initially grows rapidly to
significantly large positive values, on a rapid time scale
$\tau\lesssim\tau_{R}$. i.e. before microscopic collisions become effective.
This is caused by the rapid longitudinal expansion which red-shifts the
longitudinal momenta and reduces the longitudinal pressure of the system
Jaiswal:2021uvv ; Jaiswal:2019cju ; Kurkela:2019set .161616This argument holds
strictly only as long as $\bar{\Pi}$ is small which applies here. At much
later times $\tau\gg\tau_{R}$ the shear stress joins its first-order Navier-
Stokes limit (here shown for $C=10/4\pi$)
$\bar{\pi}_{\mathrm{NS}}=4/3\times[\beta_{\pi}/(e+P)]/(\tau/\tau_{R})$ where
the dimensionless quantity $\beta_{\pi}/(e{+}P)$ is a function of $\mu/T$ and
$m/T$ that approaches 1/5 in the massless limit.

Panel (b) of Fig. 9 shows that the bulk viscous pressure evolves in a rather
complex fashion which, for the case of equilibrium initial conditions, even
features oscillations driven by the bulk-shear coupling terms in Eqs. (3,4)
(qualitatively consistent with hydrodynamic results obtained in
Jaiswal:2014isa ; Jaiswal:2022udf  for a single-species massive Boltzmann gas
at vanishing chemical potential).171717Due to the smallness of $\bar{\Pi}$ its
back-reaction on the evolution of $\bar{\pi}$ is negligible. Also, it
approaches its first-order Navier-Stokes limit
$\bar{\Pi}_{\mathrm{NS}}=-[\beta_{\Pi}/(e{+}P)]/(\tau/\tau_{R})$ (shown as
dotted lines for both choices of $C$ in Fig. 9b) only at much later times than
the shear stress (even outside the time range shown here). We checked that
this is, too, is caused by bulk-shear coupling Denicol:2014vaa  with a shear
stress that exceeds the magnitude of $\bar{\Pi}$ by more than an order of
magnitude.

###  V.2 Kinetic theory for a non-conformal quark-gluon gas with net baryon
charge

####  V.2.1 Formal solution

The RTA Boltzmann equation (20) for a system of gluons and massive quarks with
Boltzmann statistics undergoing Bjorken flow has been solved in Refs.
Florkowski:2014sfa ; Florkowski:2014txa . For reasons of uniformity of
notation we briefly review this solution, generalizing it along the way to
include quantum statistics in the distribution functions. The solution will be
expressed in terms of macroscopic quantities, i.e. the temperature, chemical
potential, shear, and bulk viscous pressures. The initial conditions for the
shear and bulk viscous pressures are generated from Eqs. (23-25). With these
initial distributions, the exact evolutions of the temperature and chemical
potential are obtained by the Landau matching conditions (26,27), generalized
to non-zero quark mass:

$\displaystyle
e_{\mathrm{eq}}(T(\tau),\mu(\tau))=D(\tau,\tau_{0})\frac{\Lambda_{0}^{4}}{4\pi^{2}}\tilde{H}_{e}\left[\frac{\tau_{0}}{\tau\sqrt{1+\xi_{0}}},\frac{m}{\Lambda_{0}},\frac{\nu_{0}}{\Lambda_{0}}\right]$
$\displaystyle+\frac{1}{4\pi^{2}}\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}}{\tau_{R}(\tau^{\prime})}\,D(\tau,\tau^{\prime})\,T(\tau^{\prime})^{4}\,\tilde{H}_{e}\left[\frac{\tau^{\prime}}{\tau},\frac{m}{T(\tau^{\prime})},\frac{\mu(\tau^{\prime})}{T(\tau^{\prime})}\right],$
(44) $\displaystyle
n_{\mathrm{eq}}(T(\tau),\mu(\tau))=\frac{\tau_{0}}{\tau}\,n_{\mathrm{in}}.$
(45)

Here

$\displaystyle\tilde{H}_{e}(y,z,\alpha)$
$\displaystyle\equiv\int_{0}^{\infty}\\!\\!\\!du\,u^{3}\Bigg{[}g_{g}\,\tilde{f}^{g}_{\mathrm{eq}}(u)\,H_{e}(y,0)$
$\displaystyle+\,g_{q}\left(\tilde{f}^{q}_{\mathrm{eq}}(u,z,\alpha)+\tilde{f}^{\bar{q}}_{\mathrm{eq}}(u,z,\alpha)\right)\,H_{e}(y,z/u)\Bigg{]},$
(46)

with $\tilde{f}_{\mathrm{eq}}$ from Eq. (8-10) above and

$\displaystyle H_{e}(y,z)\equiv
y\left(\sqrt{y^{2}+z^{2}}+\frac{1+z^{2}}{\sqrt{y^{2}-1}}\,\tanh^{-1}\sqrt{\frac{y^{2}-1}{y^{2}+z^{2}}}\right).$
(47)

Eqs. (44-45) are solved by numerical iteration to obtain solutions for
$T(\tau)$ and $\mu(\tau)$. Once these are known the distribution function
$f(\tau;p_{T},w)$ itself can be determined at any $\tau$ using Eq. (IV.2), and
all macroscopic hydrodynamic quantities can be obtained by taking appropriate
moments of $f$. For example, the out-of-equilibrium entropy density is given
by Eq. (35). The shear and bulk viscous stresses are obtained from the
effective longitudinal and transverse pressures, $P_{L}$ and $P_{T}$, using
$\pi\equiv 2(P_{T}-P_{L})/3$ and $\Pi\equiv(P_{L}+2P_{T}-3P)/3$, where

$\displaystyle
P_{L,T}=D(\tau,\tau_{0})\frac{\Lambda_{0}^{4}}{4\pi^{2}}\tilde{H}_{L,T}\left[\frac{\tau_{0}}{\tau\sqrt{1+\xi_{0}}},\frac{m}{\Lambda_{0}},\frac{\nu_{0}}{\Lambda_{0}}\right]$
$\displaystyle+\frac{1}{4\pi^{2}}\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}}{\tau_{R}(\tau^{\prime})}\,D(\tau,\tau^{\prime})\,T^{4}(\tau^{\prime})\,\tilde{H}_{L,T}\left[\frac{\tau^{\prime}}{\tau},\frac{m}{T(\tau^{\prime})},\frac{\mu(\tau^{\prime})}{T(\tau^{\prime})}\right],$
(48)

with

$\displaystyle\tilde{H}_{L,T}(y,z,\alpha)\equiv\int_{0}^{\infty}\\!\\!\\!du\,u^{3}\Bigg{[}g_{g}\tilde{f}^{g}_{\mathrm{eq}}(u)H_{L,T}(y,0)$
$\displaystyle+\,g_{q}\left(\tilde{f}^{q}_{\mathrm{eq}}(u,z,\alpha)+\tilde{f}^{\bar{q}}_{\mathrm{eq}}(u,z,\alpha)\right)H_{L,T}\left(y,\frac{z}{u}\right)\Bigg{]},$
(49)

where

$\displaystyle H_{L}(y,z)$
$\displaystyle=\frac{y^{3}}{(y^{2}-1)^{3/2}}\Bigg{(}\sqrt{(y^{2}-1)(y^{2}+z^{2})}$
$\displaystyle-(z^{2}+1)\,\tanh^{-1}\sqrt{\frac{y^{2}-1}{y^{2}+z^{2}}}\Bigg{)},$
(50) $\displaystyle H_{T}(y,z)$
$\displaystyle=\frac{y}{2(y^{2}-1)^{3/2}}\Bigg{(}-\sqrt{(y^{2}-1)(y^{2}+z^{2})}$
$\displaystyle+(z^{2}+2y^{2}-1)\,\tanh^{-1}\sqrt{\frac{y^{2}-1}{y^{2}+z^{2}}}\Bigg{)}.$
(51)

####  V.2.2 Numerical results

Figure 10a shows the expansion trajectories obtained from non-conformal
kinetic theory, for the same initial conditions as for the hydrodynamic
solutions presented in the preceding subsection, and the associated scaled-
time evolution of the shear and bulk viscous stresses is shown in Fig. 11.
Line styles and colors have the same meaning as for the corresponding
hydrodynamic solutions shown in Figs. 8, 9.

Figure 10: (a) Expansion trajectories from non-conformal kinetic theory with
quark mass $m=1$ GeV. Initial conditions are the same and line colors and
styles have the same meaning as in Fig. 8a. (b) Solid lines: evolution of the
full (non-equilibrium + equilibrium) specific entropy $s_{\mathrm{tot}}/n$
obtained from kinetic theory for the identically colored solid lines in panel
(a). Dotted lines show the corresponding equilibrium contributions
$s_{\mathrm{eq}}/n$ for comparison.

The initial Romatschke-Strickland parameters for the non-equilibrium initial
conditions used to generate the solid lines are listed Table 2.

Set |  $\Lambda_{0}$ (GeV) |  $\nu_{0}$ (GeV) | $\xi_{0}$ | $\bar{\phi}_{0}$ | $\psi_{nc,0}$
---|---|---|---|---|---
I | $0.342$ | $0.597$ | $-0.981$ | $-0.45$ | $-0.301$
II | $0.209$ | $1.001$ | $-0.986$ | $-0.45$ | $-0.269$

Table 2: Initial parameters $(\Lambda_{0},\nu_{0},\xi_{0})$ used to generate
the blue and black solid curves in the sets I and II of Fig. 10. As before
$\bar{\phi}_{0}$ denotes the corresponding initial value of
$\bar{\phi}=\bar{\pi}-\bar{\Pi}$. $\psi_{nc,0}$ denotes the maximum value of
$\bar{\phi}_{0}$ that leads to trajectories for which initially
$d\mu/d\tau>0$.

The qualitative and quantitative similarities between the kinetic theory
results shown in Figs. 10a, 11 and the hydrodynamic results shown earlier in
Figs. 8, 9 are obvious. Some of the quantitative differences are discussed in
Appendix C but they do not impact our qualitative observations. Additional
insights into the mechanisms at work in Fig. 10a can be gleaned from Fig.
10b.181818Lacking an expression for the non-equilibrium entropy density in
non-conformal second-order dissipative hydrodynamics, there was no analogous
plot provided for the hydrodynamic solutions presented in Sec. V.2.1. In
kinetic theory Eq. (35) provides an unambiguous definition of the total
entropy density, and it can be easily split into an equilibrium part and a
non-equilibrium correction. For the expansion trajectories corresponding to
the set of far-from-equilibrium initial conditions listed in Table 2 and shown
as solid lines in Fig. 10a, panel (b) plots the scaled time evolution of the
full specific entropy using solid lines, together with the corresponding
equilibrium contributions $s_{\mathrm{eq}}/n$ using dotted lines. Note that,
up to overall larger specific entropies at smaller $\mu/T$ values, the
qualitative characteristics of the curves shown in Fig. 10b are the same for
initial condition sets I and II, and we will therefore discuss them together.

These expansion trajectories are characterised by large initial dissipative
fluxes, resulting in large negative non-equilibrium corrections to the initial
specific entropy (i.e. the solid lines start at much lower $s/n$ values than
the dotted lines). The large negative deformation parameters $\xi_{0}$, close
to their absolute lower limits $\xi_{\mathrm{min}}=-1$ (see Table 2), imply
that initially the longitudinal local momentum distribution is much wider than
the transverse one and almost flat. Accordingly, the total $(s/n)_{0}$ is much
smaller than its initial equilibrium contribution $(s_{\mathrm{eq}}/n)_{0}$.
Note that, for both sets I and II, the starting values of $s/n$ (both
equilibrium and non-equilibrium parts) are the same for the blue and black
lines; they start, however, at different $\tau/\tau_{R}$ due to different
values for $C$. At first sight, the steep initial drop of all the dotted
$s_{\mathrm{eq}}/n$ lines seems to indicate a rapid approach towards local
thermal equilibrium, on a time scale $\tau\lesssim\tau_{R}$, but once
$s_{\mathrm{eq}}/n$ reaches the solid line showing the total $s/n$, it doesn’t
stay there but “bounces off”, settling on the solid line for good only much
later at $\tau\gg\tau_{R}$.

The clue to understanding this phenomenon is the observation that the initial
“fake” equilibrium state (when $s_{\mathrm{tot}}/n$ and $s_{\mathrm{eq}}/n$
coincide for the first time) happens at $\tau<\tau_{R}$, i.e. before
microscopic collisions have had much of an effect. It is rather a consequence
of approximate free-streaming at very early times in Bjorken flow
Jaiswal:2021uvv ; Jaiswal:2019cju ; Kurkela:2019set ; Chattopadhyay:2021ive
which red-shifts the longitudinal momenta, causing the momentum distribution
to eventually cross from its initially strongly prolate ($P_{L}>P_{T}$) to an
oblate ($P_{L}<P_{T}$) ellipsoidal shape.191919Note that a non-interacting
(i.e. free-streaming) system initialised with Romatschke-Strickland
distributions characterised by two energy scales $(\Lambda_{0},\nu_{0})$ and a
negative $\xi_{0}$ passes through such a “fake thermal equilibrium” state.
This happens at $\tau_{\mathrm{fte}}=\tau_{0}/\sqrt{1+\xi_{0}}$ when the free-
streaming distribution function takes a thermal form:
$f^{i}(\tau_{\mathrm{fte}};p_{T},p_{z})=f_{\mathrm{RS}}^{i}\bigl{(}\tau;p_{T},(\tau_{\mathrm{fte}}/\tau_{0})p_{z}\bigr{)}=f_{\mathrm{eq}}^{i}(\Lambda_{0},\nu_{0};p)$.
Near the crossing point, when the momentum distribution is spherical, the full
specific entropy $s_{\mathrm{tot}}/n$ approximately202020i.e. up to deviations
caused by a non-exponential energy dependence and a non-equilibrium
normalization value. coincides with its equilibrium part $s_{\mathrm{eq}}/n$
but since collisions require more time they cannot keep the system in this
“fake” equilibrium state — instead, the momentum distribution becomes
anisotropic again by turning oblate. Only much later are collisions able to
locally isotropize the momentum distribution and keep it there. At that stage
the dotted lines for $s_{\mathrm{eq}}$ merge with the solid lines for
$s_{\mathrm{tot}}$, and both measures for $s/n$ saturate at their common
asymptotic, global thermal equilibrium value.

The time when $s_{\mathrm{eq}}/n$ approaches its minimum roughly marks the
point where viscous cooling turns into viscous heating and the initially
rightward moving expansion trajectories in Fig. 10a turn again towards the
left of the phase diagram.

Figure 11: Same as Fig. 9 but for the solutions obtained from kinetic theory.

The solid lines show that the total specific entropy grows in two spurts: The
first of these (the viscous cooling stage) occurs early, before the system
reaches its “fake” equilibrium state, and is driven by longitudinal momentum-
stratification via free-streaming. The second spurt is caused by local
thermalization via microscopic interactions which lead to viscous
heating.212121We note that the first spurt is larger for the more strongly
coupled (less viscous) system with $C=3/4\pi$ but the second spurt and overall
entropy production is higher for the less strongly coupled (more viscous)
system with $C=10/4\pi$. From the work of Israel and Stewart Israel:1979wp  it
is known that, in viscous hydrodynamics, the rate of entropy production is
quadratic in the bulk and shear stresses,
$\propto\pi^{2}/2\eta+\Pi^{2}/\zeta$. The evolution of $\pi$ and $\Pi$ for the
initial condition set II is shown in Fig. 11. Comparison of the solid lines
for set II in Fig.10b with those in Fig. 11 confirms that the plateau for
$s/n$ indeed lies close to the time when $\pi$ passes through zero, although
slightly shifted by the contribution from the bulk viscous entropy production
rate.

###  V.3 (Absence of) early-time universality

In Refs. Jaiswal:2021uvv ; Chattopadhyay:2021ive  it was shown that in non-
conformal RTA kinetic theory with Bjorken flow there is an early-time, far-
from-equilibrium attractor for the longitudinal pressure $P_{L}=P+\Pi-\pi$
(driven by the approximate free-streaming dynamics at early times) while no
such attractor exists for the shear and bulk viscous stresses $\pi$ and $\Pi$
separately. That the early-time evolution of $\pi$ and $\Pi$ is not ruled by a
universal attractor at early times is seen in Fig. 11 which shows strong
sensitivity to the initial conditions, with convergence to the first-order
Navier-Stokes attractor only at very late times $\tau/\tau_{R}\gg 1$ (i.e.
small Knudsen number).222222This brief analysis does not rule out the
existence of an early-time attractor in the full phase-space of variables,
$(T,\mu,\pi,\Pi)$, as often considered in the theory of non-autonomous
dynamical systems. Attractors using such generalized definitions have been
pursued for conformal and non-conformal fluids at vanishing chemical potential
Behtash:2017wqg ; Kamata:2022jrc .  In contrast,

Figure 12: Evolution of the scaled effective longitudinal pressure $P_{L}/P$
as a function of the inverse Knudsen number $\tau/\tau_{R}$. The red-dotted
line is the Navier-Stokes solution obtained with initial $(T_{0},\mu_{0})$
corresponding to set II of Fig. 8.

an early-time, far-from-equilibrium attractor for $P_{L}/P$, to which all
initial conditions converge on very short time scales $\tau\lesssim\tau_{R}$,
is seen in Fig. 12. Figure 14 in Appendix C demonstrates that this attractor
is not reproduced by the hydrodynamic approximation studied in Sec. V.1. All
this is consistent with earlier work Jaiswal:2021uvv ; Chattopadhyay:2021ive
and extends it to systems with non-zero net baryon number.

##  VI Conclusions

In this paper we explored the evolution of phase trajectories of a system of
quarks and gluons undergoing Bjorken expansion, comparing second-order
Chapman-Enskog hydrodynamics and kinetic theory. Using this bare bones model
we were able to reproduce the findings in Ref. Dore:2020jye  which, under
certain circumstances, suggest at first sight a loss of entropy, violating the
second law of thermodynamics. This phenomenon is accompanied by a novel
feature of far-off-equilibrium expansion dynamics which we call viscous
cooling: the dissipative fluid cools more rapidly than an ideal one, contrary
to the more common internal friction effect in dissipative fluids known as
viscous heating. We found that viscous cooling arises only for far-off-
equilibrium initial conditions where both the shear and bulk viscous stresses
have the opposite signs of their Navier-Stokes values, to which they relax
over time. Its most dramatic manifestation, where the cooling system initially
moves towards larger $\mu_{B}$ whereas an ideal fluid would have evolved
towards smaller $\mu_{B}$, requires non-zero quark masses, large negative
shear stress and (relatively) large positive initial bulk viscous
pressures.232323We stress that the natural scales for the (smaller) bulk and
(larger) shear viscous pressures differ by about an order of magnitude in our
system, see Figs. 9 and 11.

As the origin of this phenomenon we identified a transfer of negative entropy
from the (large) dissipative correction in the initial state to the
equilibrium sector. For the massless (conformal) case we succeeded in deriving
a macroscopic expression for the dissipative entropy correction, and hence we
were able to check this mechanism in both the macroscopic hydrodynamic and
microsopic kinetic theory approaches. For the massive (non-conformal) system
our detailed quantitative analysis relied entirely on the unambiguous
microscopic definition (35) from kinetic theory. On a qualitative level,
though, the sign (i.e. the direction) of this entropy transfer can in both
cases be obtained from the macroscopic equation (43). We confirmed that the
total entropy (i.e. the sum of the equilibrium contribution and the non-
equilibrium correction) always increases with time, as required by the
H-theorem.

It is worth contemplating how the phenomenon of viscous (or more generally:
dissipative) cooling generalizes to systems undergoing arbitrary 3-dimensional
expansion, without Bjorken symmetry. A particularly interesting problem
arising in this context is to identify a general criterion for far-off-
equilibrium initial conditions (for both the flow velocity profile and the
dissipative fluxes) that lead to dissipative cooling. For systems that
ultimately move towards local thermal equilibrium we expect dissipative
cooling to be a transient phenomenon that turns to dissipative heating once
the dissipative fluxes approach their Navier-Stokes values and become weak.
But there may exist (possibly externally driven) flow profiles where this does
not happen. It will be very interesting to establish conditions for
dissipative cooling that can be realized in laboratory experiments and make
the phenomenon directly observable. These questions are the subject of ongoing
research.

Acknowledgements: The authors thank Lipei Du, Amaresh Jaiswal, and Sunil
Jaiswal for fruitful discussions. This work was supported by the U.S.
Department of Energy, Office of Science, Office for Nuclear Physics under
Awards No. DE-SC0004286 (U.H.) and DE-FG02-03ER41260 (T.S. and C.C.).

##  Appendix A Second-order hydrodynamics for a massive quark-gluon gas

The energy-momentum tensor of a gas of quarks, anti-quarks, and gluons is

$\displaystyle T^{\mu\nu}=\sum_{i=1}^{3}\int
dP_{i}\,p^{\mu}_{i}\,p^{\nu}_{i}\,f^{i}=e\,u^{\mu}u^{\nu}-\left(P+\Pi\right)\Delta^{\mu\nu}+\pi^{\mu\nu},$
(52)

where $i$ differentiates between the species which have masses $m_{1}=m_{2}=m$
(quarks and antiquarks) and $m_{3}=0$ (gluons), and $dP_{i}\equiv
g_{i}d^{3}p/[(2\pi)^{3}E_{p,i}]$ with $E_{p,i}=\sqrt{p^{2}+m_{i}^{2}}$. In the
constitutive relation, $e$ is the energy density in the fluid rest frame, $P$
and $\Pi$ are the equilibrium and bulk viscous pressures, respectively, and
$\pi^{\mu\nu}$ is the shear stress tensor. The projection tensors
$u^{\mu}u^{\nu}$ and $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ select the
components of $T^{\mu\nu}$ along the temporal and spatial directions in the
local fluid rest frame, defined by $u^{\mu}_{\mathrm{LRF}}=(1,0,0,0)$. The
conserved net-quark charge current is given by

$\displaystyle N^{\mu}=n\,u^{\mu}+n^{\mu}=N^{\mu}_{1}-N^{\mu}_{2}$ (53)

where $n$ is the LRF net quark density and $n^{\mu}$ the quark number
diffusion current, with

$\displaystyle N^{\mu}_{i}=\int dP_{i}\,p^{\mu}_{i}\,f^{i}.$ (54)

In the Landau matching scheme the energy and net quark densities are written
in terms of the local equilibrium distributions as

$\displaystyle e$
$\displaystyle=u_{\alpha}\,u_{\beta}\,\sum_{i=1}^{3}\int\,dP_{i}\,p^{\alpha}_{i}\,p^{\beta}_{i}\,f_{\mathrm{eq}}^{i},$
(55) $\displaystyle n$
$\displaystyle=u_{\alpha}\,\sum_{i=1}^{2}\,(-1)^{i-1}\,\int\,dP_{i}\,p^{\alpha}_{i}\,f_{\mathrm{eq}}^{i},$
(56)

which defines the values for the inverse temperature and chemical potential in
the latter:

$\displaystyle f_{\mathrm{eq}}^{i}=\frac{1}{\exp\left(\beta
u{\,\cdot\,}p_{i}-\alpha_{i}\right)-\theta_{i}}.$ (57)

Here $\beta=1/T$ is the inverse temperature and $\alpha_{i}=\mu_{i}/T$ the
normalized chemical potential of species $i$ ($\alpha_{1}=-\alpha_{2}=\alpha$,
$\alpha_{3}=0$), and $\theta_{1}=\theta_{2}=-1$, $\theta_{3}=1$ distinguishes
between Fermi-Dirac and Bose-Einstein statistics. The dynamics of the system
is governed by conservation of energy-momentum and charge current:
$\partial_{\mu}T^{\mu\nu}=0$ and $\partial_{\mu}N^{\mu}=0$. Projection of the
former along the fluid velocity, $u_{\mu}\partial_{\nu}T^{\mu\nu}=0$, and the
latter yield the comoving evolution of energy and number densities, whereas
projection of the former orthogonal to the flow,
$\Delta^{\gamma}_{\nu}\partial_{\mu}T^{\mu\nu}=0$, determines the acceleration
of a fluid element:

$\displaystyle\dot{e}$
$\displaystyle=-\left(e+P\right)\theta-\Pi\,\theta+\pi^{\mu\nu}\,\sigma_{\mu\nu},$
(58) $\displaystyle\dot{n}$ $\displaystyle=-n\theta-\partial_{\mu}n^{\mu}.$
(59) $\displaystyle\left(e{+}P\right)\dot{u}^{\gamma}$
$\displaystyle=\nabla^{\gamma}P_{\mathrm{eq}}-\Pi\,\dot{u}^{\gamma}+\nabla^{\gamma}\Pi-\Delta^{\gamma}_{\nu}\partial_{\mu}\pi^{\mu\nu}.$
(60)

Here we defined the comoving time derivative $\dot{A}=u^{\mu}\partial_{\mu}A$,
the spatial derivative in the fluid rest-frame
$\nabla^{\mu}=\Delta^{\mu\nu}\partial_{\nu}$, the fluid expansion rate
$\theta=\nabla_{\mu}u^{\mu}$, and the velocity shear tensor
$\sigma^{\mu\nu}=\Delta^{\mu\nu}_{\alpha\beta}\,\nabla^{\alpha}u^{\beta}$,
where

$\displaystyle\Delta^{\mu\nu}_{\alpha\beta}=\left(\Delta^{\mu}_{\alpha}\,\Delta^{\nu}_{\beta}+\Delta^{\mu}_{\beta}\,\Delta^{\nu}_{\alpha}\right)/2-\frac{1}{3}\,\Delta^{\mu\nu}\,\Delta_{\alpha\beta}$
(61)

projects any rank-2 tensor on its traceless part and locally spatial
components.

In Eqs. (58-60), the first terms on r.h.s. stem from ideal hydrodynamic
constitutive relations, while the others arise due to dissipation. We shall
eventually express comoving derivatives of energy density and net-quark
density in terms of derivatives of inverse temperature $\beta$ and normalised
chemical potential $\alpha$. For this, we use the definitions given by Eqs.
(55) to write

$\displaystyle\dot{e}$
$\displaystyle=-\dot{\beta}\,I_{3,0}^{(0)}-\dot{\alpha}\sum_{i=1}^{2}\,(-1)^{i}\,I_{2,0}^{(0),i},$
$\displaystyle\dot{n}$
$\displaystyle=\dot{\beta}\sum_{i=1}^{2}\,(-1)^{i}\,I_{2,0}^{(0),i}+\dot{\alpha}\sum_{i=1}^{2}I_{1,0}^{(0),i},$
(62)

where we used the notation $I_{n,q}^{(r)}=\sum_{i=1}^{3}I_{n,q}^{(r),i}$,
with242424The superscript index $r$ is somewhat redundant in the notation as
the moment depends only on the difference $n{-}r$, and one could have simply
defined, $J_{n,q}^{i}=I_{n,q}^{(0),i}$. However, we continue using it for ease
of comparison of our transport coefficients with those obtained previously for
a single component massive Boltzmann gas at vanishing chemical potential
Jaiswal:2014isa .

$\displaystyle I_{n,q}^{(r),i}\equiv\frac{1}{(2q+1)!!}\int
dP_{i}\,\left(u\cdot
p_{i}\right)^{n-r-2q}\,\left(\Delta_{\alpha\beta}\,p^{\alpha}_{i}\,p^{\beta}_{i}\right)^{q}\,{\cal
F}^{i}_{\mathrm{eq}}.$ (63)

The function ${\cal F}^{i}_{\mathrm{eq}}$ stems from taking a derivative of
$f_{\mathrm{eq}}^{i}$ and is defined as ${\cal
F}^{i}_{\mathrm{eq}}=f_{\mathrm{eq}}^{i}\tilde{f}_{\mathrm{eq}}^{i}$ where
$\tilde{f}_{\mathrm{eq}}^{i}=1+a_{i}f_{\mathrm{eq}}^{i}$.252525For later use
we note that the function ${\cal F}^{i}_{\mathrm{eq}}$ has the following
properties: ${\cal F}^{i}_{\mathrm{eq}}=-1/(u\cdot p_{i})\times\partial
f^{i}_{\mathrm{eq}}/\partial\beta$, ${\cal F}^{q}_{\mathrm{eq}}=\partial
f^{q}_{\mathrm{eq}}/\partial\alpha$, and ${\cal
F}^{\bar{q}}_{\mathrm{eq}}=-\partial
f^{\bar{q}}_{\mathrm{eq}}/\partial\alpha$.

Using Eqs. (A) along with the time evolution equations for $e$ and $n$ in Eqs.
(58-59) we obtain

$\displaystyle\dot{\beta}$ $\displaystyle=-{\cal G}_{\beta}\,\theta-{\cal
H}_{\beta}\,\left(\Pi\,\theta-\pi^{\mu\nu}\,\sigma_{\mu\nu}\right),$
$\displaystyle\dot{\alpha}$ $\displaystyle=-{\cal G}_{\alpha}\,\theta-{\cal
H}_{\alpha}\,\left(\Pi\,\theta-\pi^{\mu\nu}\,\sigma_{\mu\nu}\right).$ (64)

The functions ${\cal G}$ and ${\cal H}$ are defined as

$\displaystyle{\cal G}_{\beta}$ $\displaystyle=\frac{1}{{\cal
D}}\Bigl{[}\left(e+P_{\mathrm{eq}}\right)\,\sum_{i=1}^{2}I_{1,0}^{(0),i}+n\,\sum_{i=1}^{2}(-1)^{i}I_{2,0}^{(0),i}\Bigr{]},$
(65) $\displaystyle{\cal G}_{\alpha}$ $\displaystyle=-\frac{1}{{\cal
D}}\Bigl{[}\left(e+P_{\mathrm{eq}}\right)\sum_{i=1}^{2}(-1)^{i}I_{2,0}^{(0),i}+n\,I_{3,0}^{(0)}\Bigr{]},$
(66) $\displaystyle{\cal H}_{\beta}$
$\displaystyle=\frac{\sum_{i=1}^{2}\,I_{1,0}^{(0),i}}{{\cal
D}},\,\,\,\,\,{\cal
H}_{\alpha}=-\frac{\sum_{i=1}^{2}(-1)^{i}\,I_{2,0}^{(0),i}}{{\cal D}},$ (67)
$\displaystyle{\cal D}$
$\displaystyle\equiv\Bigl{(}\sum_{i=1}^{2}(-1)^{i}\,I_{2,0}^{(0),i}\Bigr{)}^{2}-I_{3,0}^{(0)}\,\sum_{i=1}^{2}I_{1,0}^{(0),i}.$
(68)

In the conformal case, the coefficients ${\cal G}_{\beta}$ and ${\cal
G}_{\alpha}$ simplify to ${\cal G}_{\beta}=-\beta/3$, ${\cal G}_{\alpha}=0$,
consistent with Jaiswal:2015mxa .262626We use the conformal relations,
$e(\beta,\alpha)$ and $n(\beta,\alpha)$, given in Eqs. (6-7) to simplify some
of the terms appearing in ${\cal G}_{\beta}$ and ${\cal G}_{\alpha}$:
$\sum_{i=1}^{2}(-1)^{i}I_{2,0}^{(0),i}=-\partial e/\partial\alpha=-3n/\beta$
and $I_{3,0}^{0}=-\partial e/\partial\beta=4e/\beta$. To compute the bulk
viscous pressure and shear stress tensor of the system we start from their
respective definitions

$\displaystyle\Pi$
$\displaystyle=-\frac{\Delta_{\alpha\beta}}{3}\,\sum_{i=1}^{3}\,\int
dP_{i}\,p^{\alpha}_{i}\,p^{\beta}_{i}\,\delta f^{i},$
$\displaystyle\pi^{\mu\nu}$
$\displaystyle=\Delta^{\mu\nu}_{\alpha_{\beta}}\,\sum_{i=1}^{3}\,\int\,dP_{i}\,p^{\alpha}_{i}\,p^{\beta}_{i}\,\delta
f^{i}.$ (69)

These off-equilibrium corrections are obtained by solving the RTA Boltzmann
equation

$\displaystyle p^{\mu}_{i}\,\partial_{\mu}\,f^{i}=-\frac{u\cdot
p_{i}}{\tau_{R}}\,\left(f^{i}-f^{i}_{\mathrm{eq}}\right)$ (70)

iteratively in powers of $\delta f^{i}=f^{i}-f_{\mathrm{eq}}^{i}$. The leading
order correction is

$\displaystyle\delta f^{i}_{(1)}=-\frac{\tau_{R}}{u\cdot
p_{i}}p^{\mu}\partial_{\mu}f_{\mathrm{eq}}^{i}=-\tau_{R}\left(\dot{f}_{\mathrm{eq}}^{i}+\frac{1}{u\cdot
p_{i}}p^{\mu}\nabla_{\mu}f_{\mathrm{eq}}^{i}\right).$ (71)

With Eq. (57) this can be written as

$\displaystyle\delta f^{i}_{(1)}$ $\displaystyle=\tau_{R}\Big{[}(u\cdot
p_{i})\dot{\beta}-\dot{\alpha_{i}}+\frac{\beta\,\theta}{3(u\cdot
p_{i})}\Delta_{\mu\nu}\,p^{\mu}_{i}\,p^{\nu}_{i}\Big{]}{\cal
F}_{\mathrm{eq}}^{i}$ $\displaystyle+\frac{\tau_{R}\,\beta}{(u\cdot
p_{i})}p^{\mu}_{i}\,p^{\nu}_{i}\,\sigma_{\mu\nu}\,{\cal F}_{\mathrm{eq}}^{i}$
$\displaystyle+\tau_{R}\Big{[}\beta\left(\dot{u}\cdot
p_{i}\right)+p^{\mu}_{i}\nabla_{\mu}\beta-\frac{p^{\mu}_{i}}{u\cdot
p_{i}}\nabla_{\mu}\alpha_{i}\,\Big{]}{\cal F}_{\mathrm{eq}}^{i}.$ (72)

The terms in the first square bracket give rise to the bulk viscous pressure,
the second term involving $\sigma^{\mu\nu}$ yields the shear stress tensor,
and the terms in the second line result in the first-order quark diffusion
current. Since it is not needed for Bjorken flow we shall neglect the quark
diffusion current in the following. Also, we convert time derivatives of
$\beta$ and $\alpha$ into velocity gradients using Eqs. (A) to obtain

$\displaystyle\delta f^{i}_{(1)}$ $\displaystyle=\tau_{R}\,\theta\Big{[}-{\cal
G}_{\beta}\,(u\cdot p_{i})+{\cal
G}^{i}_{\alpha}+\frac{\beta\Delta_{\mu\nu}}{3(u\cdot
p_{i})}\,p^{\mu}_{i}\,p^{\nu}_{i}\Big{]}{\cal F}_{\mathrm{eq}}^{i}$
$\displaystyle+\frac{\tau_{R}\,\beta}{(u\cdot
p_{i})}p^{\mu}_{i}\,p^{\nu}_{i}\,\sigma_{\mu\nu}\,{\cal F}_{\mathrm{eq}}^{i}$
$\displaystyle+\tau_{R}\,\left(\Pi\,\theta-\pi^{\mu\nu}\,\sigma_{\mu\nu}\right)\,\left({\cal
H}_{\alpha}^{i}-{\cal H}_{\beta}\,(u\cdot p_{i})\right);$ (73)

here we defined ${\cal G}_{\alpha}^{1}=-{\cal G}_{\alpha}^{2}={\cal
G}_{\alpha},{\cal G}_{\alpha}^{3}=0$ and ${\cal H}_{\alpha}^{1}=-{\cal
H}_{\alpha}^{2}={\cal H}_{\alpha},{\cal H}_{\alpha}^{3}=0$, with ${\cal
G}_{\alpha}$ and ${\cal H}_{\alpha}$ given by Eqs. (66-67). Note that the
terms in the first line of the r.h.s. of Eq. (73) are of first order in
velocity gradients whereas those in the second line are of second order. Using
Eq. (73) in the definitions (69) of the bulk viscous pressure and shear stress
tensor and retaining terms only up to first order in gradients, we obtain the
Navier-Stokes expressions for $\Pi$ and $\pi^{\mu\nu}$:

$\displaystyle\Pi=-\,\tau_{R}\,\beta_{\Pi}\,\theta,\qquad\pi^{\mu\nu}=2\tau_{R}\,\beta_{\pi}\,\sigma^{\mu\nu}.$
(74)

The coefficients $\beta_{\Pi}$ and $\beta_{\pi}$ are given by272727We have
checked that in the conformal limit $\beta_{\pi}\to 4P/5$ as obtained in
Jaiswal:2015mxa  for a conformal quark-gluon gas at finite $\mu_{B}$. Also, in
this limit $e=3P$ and $c_{s}^{2}=1/3$ such that $\beta_{\Pi}\to 0$ as
expected.

$\displaystyle\beta_{\Pi}$
$\displaystyle=\frac{5}{3}\,\beta\,I_{4,2}^{(1)}\,+\sum_{i=1}^{3}\,{\cal
G}_{\alpha}^{i}\,I_{2,1}^{(0),i}-{\cal G}_{\beta}\,I_{3,1}^{(0)}$
$\displaystyle=\frac{5}{3}\beta_{\pi}-\left(e+P\right)c_{s}^{2},$
$\displaystyle\beta_{\pi}$ $\displaystyle=\beta\,I_{4,2}^{(1)}.$ (75)

We simplified the coefficient $\beta_{\Pi}$ using the definition of squared
speed of sound,282828To obtain the following relation, we write
$dP/de=\partial P/\partial\beta\times d\beta/de+\partial
P/\partial\alpha\times d\alpha/de$, and then use the condition of fixed
specific entropy, i.e, $d(s/n)_{\mathrm{eq}}=0\implies n\,de=(e+P)\,dn$, to
express the derivatives $d\beta$ and $d\alpha$ in terms of $de$.

$\displaystyle c_{s}^{2}\equiv\left(\frac{\partial P}{\partial
e}\right)_{s/n}$ $\displaystyle=\frac{1}{(e+P){\cal D}}\Bigg{[}\frac{\partial
P}{\partial\beta}\,\,\left((e+P)\,\frac{\partial
n}{\partial\alpha}+n\,\frac{\partial n}{\partial\beta}\right)$
$\displaystyle+\frac{\partial P}{\partial\alpha}\left((e+P)\,\frac{\partial
e}{\partial\alpha}+n\,\frac{\partial e}{\partial\beta}\right)\Bigg{]},$ (76)

where ${\cal D}$ is given by Eq. (68). Moreover, using the properties of
${\cal F}^{i}_{\mathrm{eq}}$ one finds the following relations:

$\displaystyle I_{2,1}^{(0),1}=-\frac{\partial P_{1}}{\partial\alpha},\qquad
I_{2,1}^{(0),2}=\frac{\partial P_{2}}{\partial\alpha},\qquad
I_{3,1}^{(0)}=\frac{\partial P}{\partial\beta}.$ (77)

Here $P_{1}$ and $P_{2}$ are the partial pressures of quarks and anti-quarks,
respectively. This yields292929To write ${\cal G}_{\beta}$ and ${\cal
G}_{\alpha}$ in terms of the quantities appearing in Eq. (A) we used Eqs. (A)
to identify $I_{3,0}^{(0)}=-\partial e/\partial\beta$,
$\sum_{i=1}^{2}(-1)^{i}\,I_{2,0}^{(0),i}=-\partial e/\partial\alpha=\partial
n/\partial\beta$, and $\sum_{i=1}^{2}I_{1,0}^{(0),i}=\partial
n/\partial\alpha$.

$\displaystyle\sum_{i=1}^{3}{\cal G}_{\alpha}^{i}\,I_{2,1}^{(0),i}-{\cal
G}_{\beta}\,I_{3,1}^{(0)}=-{\cal G}_{\alpha}\,\frac{\partial
P}{\partial\alpha}-{\cal G}_{\beta}\,\frac{\partial
P}{\partial\beta}=-(e{+}P)\,c_{s}^{2}.$ (78)

To derive second-order relaxation-type evolution equations for the bulk and
shear stresses one takes comoving time-derivative of Eqs. (69) and then uses
the RTA Boltzmann equation; one finds

$\displaystyle\dot{\Pi}=-\frac{\Pi}{\tau_{R}}$ (79)
$\displaystyle+\frac{\Delta_{\alpha\beta}}{3}\sum_{i=1}^{3}\int
dP_{i}\,p^{\alpha}_{i}\,p^{\beta}_{i}\bigg{[}-\frac{\delta
f^{i}_{(1)}}{\tau_{R}}+\frac{1}{u\cdot
p_{i}}\,p^{\gamma}_{i}\nabla_{\gamma}\delta f^{i}\bigg{]},$
$\displaystyle\dot{\pi}^{\langle\mu\nu\rangle}=-\frac{\pi^{\langle\mu\nu\rangle}}{\tau_{R}}$
(80) $\displaystyle-\Delta^{\mu\nu}_{\alpha\beta}\sum_{i=1}^{3}\int
dP_{i}\,p^{\alpha}_{i}\,p^{\beta}_{i}\bigg{[}-\frac{\delta
f^{i}_{(1)}}{\tau_{R}}+\frac{1}{u\cdot
p_{i}}\,p^{\gamma}_{i}\nabla_{\gamma}\delta f^{i}\bigg{]}.$

The momentum integration over $\delta f^{i}_{(1)}$ is performed after using
Eq. (73); as we are interested in computing the dissipative evolution to
second order in gradients we also include the second-order terms in Eq. (73):

$\displaystyle\dot{\Pi}$
$\displaystyle=-\frac{\Pi}{\tau_{R}}\,-\beta_{\Pi}\,\theta+{\cal
K}_{0}\left(\Pi\,\theta-\pi^{\mu\nu}\,\sigma_{\mu\nu}\right)+\frac{1}{3}\Delta_{\alpha\beta}\,{\cal
K}^{\alpha\beta},$ (81) $\displaystyle\dot{\pi}^{\langle\mu\nu\rangle}$
$\displaystyle=-\frac{\pi^{\mu\nu}}{\tau_{R}}+2\beta_{\pi}\,\sigma^{\mu\nu}\,-\Delta^{\mu\nu}_{\alpha\beta}\,{\cal
K}^{\alpha\beta},$ (82)

with ${\cal K}_{0}\equiv\Big{(}{\cal
H}_{\beta}\,I_{3,1}^{(0)}-\sum_{i=1}^{3}\,{\cal
H}_{\alpha}^{i}\,I_{2,1}^{(0),i}\Big{)}$.303030Using ${\cal
H}_{\beta}=(1/{\cal D})\times\partial n/\partial\alpha$, ${\cal
H}_{\alpha}=(1/{\cal D})\times\partial e/\partial\alpha$, along with the
definitions (77), the coefficient ${\cal K}_{0}$ simplifies to ${\cal
K}_{0}=(\partial P/\partial e)_{n}$. We introduced a rank two tensor ${\cal
K}^{\alpha\beta}$ with contributions from a rank four tensor,
$\rho^{\alpha\beta\gamma\delta}_{2}$, and from gradients of a rank three
tensor, $\rho^{\alpha\beta\gamma}_{1}$:

$\displaystyle{\cal
K}^{\alpha\beta}\equiv\left(\nabla_{\gamma}\,\rho^{\alpha\beta\gamma}_{1}+\,\sigma_{\gamma\delta}\,\rho^{\alpha\beta\gamma\delta}_{2}+\frac{\theta}{3}\,\Delta_{\gamma\delta}\,\rho^{\alpha\beta\gamma\delta}_{2}\right).$
(83)

These tensors are defined in analogy to $I^{\mu_{1}\cdots\mu_{n}}_{q}$, by
$\rho^{\mu_{1}\dots\mu_{n}}_{q}\equiv\sum_{i=1}^{3}\rho^{\mu_{1}\cdots\mu_{n},i}_{q}$.
The $\rho$-tensors for each species stem from deviations of the corresponding
distribution functions from equilibrium:

$\displaystyle\rho^{\mu_{1}\cdots\mu_{n},i}_{q}=\int\frac{dP_{i}}{(u\cdot
p_{i})^{q}}\,p^{\mu_{1}}_{i}\cdots p^{\mu_{n}}_{i}\,\delta f^{i}.$ (84)

In the rest of this Appendix we compute those
$\rho^{\mu_{1}\cdots\mu_{n}}_{q}$ that are required to obtain ${\cal
K}^{\alpha\beta}$. We start by noting from Eq. (83) that ${\cal
K}^{\alpha\beta}$ is at least of second order in gradients. This implies that
we can simply use the first-order terms of $\delta f^{i}_{(1)}$ in Eq. (73).
However, it is customary to re-write $\delta f^{i}_{(1)}$ such that it
contains only dissipative fluxes, $\Pi$ and $\pi^{\mu\nu}$, and not gradients
of velocity (like $\theta$ and $\sigma^{\mu\nu}$). This is implemented by
substituting the Navier-Stokes equations (74) into the expression (73) for
$\delta f^{i}_{(1)}$ Jaiswal:2013npa . We denote the resulting correction as
$\delta f^{i}_{(1),r}$ where the subscript $r$ denotes re-summation:

$\displaystyle\delta f^{i}_{(1),r}$
$\displaystyle=\frac{\Pi}{\beta_{\Pi}}\Big{[}{\cal G}_{\beta}\,(u\cdot
p_{i})-{\cal G}^{i}_{\alpha}-\frac{\beta}{3(u\cdot
p_{i})}\,\Delta_{\mu\nu}\,p^{\mu}_{i}\,p^{\nu}_{i}\Big{]}{\cal
F}^{i}_{\mathrm{eq}}$ $\displaystyle+\frac{\beta}{2\,\beta_{\pi}\,(u\cdot
p_{i})}p^{\mu}_{i}\,p^{\nu}_{i}\,\pi_{\mu\nu}\,{\cal F}^{i}_{\mathrm{eq}}.$
(85)

With this we express the $\rho$-tensors as momentum integrals over ${\cal
F}^{i}$:

$\displaystyle\rho^{\mu_{1}\cdots\mu_{n}}_{q}$
$\displaystyle=\frac{\Pi}{\beta_{\Pi}}\,{\cal
G}_{\beta}\,I^{\mu_{1}\cdots\mu_{n}}_{q-1}-\frac{\Pi}{\beta_{\Pi}}\sum_{i=1}^{2}\,{\cal
G}^{i}_{\alpha}\,I^{\mu_{1}\cdots\mu_{n},i}_{q}$
$\displaystyle-\frac{\Pi\beta}{3\beta_{\Pi}}\Delta_{\gamma\delta}\,I^{\mu_{1}\cdots\mu_{n}\,\gamma\,\delta}_{q+1}+\frac{\beta}{2\beta_{\pi}}\,\pi_{\gamma\delta}\,I^{\mu_{1}\cdots\mu_{n}\,\gamma\,\delta}_{q+1}.$
(86)

To obtain the evolution of the bulk viscous pressure we evaluate the
projection of the tensor ${\cal K}^{\alpha\beta}$ along
$\Delta_{\alpha\beta}$. For the sake of clarity we show the contributions from
the term involving the gradient of $\rho^{\alpha\beta\gamma}_{1}$ and those
involving $\rho^{\alpha\beta\gamma\delta}_{2}$ separately:

$\displaystyle\frac{1}{3}\,\Delta_{\alpha\beta}\,\nabla_{\gamma}\,\rho^{\alpha\beta\gamma}_{1}={\cal
K}_{1}\,\Pi\,\theta+\frac{2}{3}\,\pi^{\mu\nu}\,\sigma_{\mu\nu},$
$\displaystyle\frac{\Delta_{\alpha\beta}}{3}\,\left(\sigma_{\gamma\delta}\,\rho^{\alpha\beta\gamma\delta}_{2}+\frac{\theta}{3}\,\Delta_{\gamma\delta}\,\rho^{\alpha\beta\gamma\delta}_{2}\right)={\cal
K}_{2}\,\Pi\,\theta+{\cal K}_{3}\,\pi^{\mu\nu}\,\sigma_{\mu\nu}.$ (87)

The coefficients occurring in these expressions are

$\displaystyle{\cal K}_{1}$ $\displaystyle=\frac{5}{3\beta_{\Pi}}\left[{\cal
G}_{\beta}\,I_{3,1}^{(0)}-\sum_{i=1}^{3}\,{\cal
G}^{i}_{\alpha}\,I_{3,1}^{(1),i}-\frac{5}{3}\,\beta_{\pi}\right]=-\frac{5}{3},$
$\displaystyle{\cal K}_{2}$ $\displaystyle=\frac{5}{3\beta_{\Pi}}\,\left[{\cal
G}_{\beta}\,I_{4,2}^{(1)}-\sum_{i=1}^{3}\,{\cal
G}_{\alpha}^{i}\,I_{4,2}^{(2),i}-\frac{7\beta}{3}\,I_{6,3}^{(3)}\right],$
$\displaystyle{\cal K}_{3}$
$\displaystyle=\frac{7\beta}{3\beta_{\pi}}\,I_{6,3}^{(3)}.$ (88)

Similarly, for the shear stress tensor evolution we compute the projection of
${\cal K}^{\alpha\beta}$ along $\Delta^{\mu\nu}_{\alpha\beta}$. This gives
rise to the contributions

$\displaystyle\Delta^{\mu\nu}_{\alpha\beta}\,\nabla_{\gamma}\,\rho^{\alpha\beta\gamma}_{1}={\cal
K}_{4}\,\Pi\,\sigma^{\mu\nu}-2\,\pi^{\langle\mu}_{\gamma}\,\omega^{\nu\rangle\gamma}+2\,\pi^{\langle\mu}_{\gamma}\,\sigma^{\nu\rangle\gamma}$
$\displaystyle+\frac{5}{3}\,\pi^{\mu\nu}\,\theta,$
$\displaystyle\Delta^{\mu\nu}_{\alpha\beta}\,\left(\sigma_{\gamma\delta}\,\rho^{\alpha\beta\gamma\delta}_{2}+\frac{\theta}{3}\,\Delta_{\gamma\delta}\,\rho^{\alpha\beta\gamma\delta}_{2}\right)={\cal
K}_{5}\,\Pi\,\sigma^{\mu\nu}$
$\displaystyle+\frac{4\beta}{\beta_{\pi}}\,I_{6,3}^{(3)}\,\pi^{\langle\mu}_{\gamma}\,\sigma^{\nu\rangle\gamma}+{\cal
K}_{6}\,\pi^{\mu\nu}\,\theta,$ (89)

with

$\displaystyle{\cal K}_{4}$ $\displaystyle=\frac{2}{\beta_{\Pi}}\,\left[{\cal
G}_{\beta}\,I_{3,1}^{(0)}-\sum_{i=1}^{3}{\cal
G}_{\alpha}^{i}\,I_{3,1}^{(1),i}-\frac{5}{3}\,\beta_{\pi}\right]=-2,$
$\displaystyle{\cal K}_{5}$ $\displaystyle=\frac{2}{\beta_{\Pi}}\,\left[{\cal
G}_{\beta}\,I_{4,2}^{(1)}-\sum_{i=1}^{3}{\cal
G}_{\alpha}^{i}\,I_{4,2}^{(2),i}-\frac{7\,\beta}{3}\,I_{6,3}^{(3)}\right]=\frac{6}{5}\,{\cal
K}_{2},$ $\displaystyle{\cal K}_{6}$
$\displaystyle=\frac{7\beta}{3\beta_{\pi}}\,I_{6,3}^{(3)}={\cal K}_{3}.$ (90)

Putting everything together we finally arrive at the following second-order
non-conformal evolution equations for the bulk and shear viscous stresses for
a quark-gluon gas:

$\displaystyle\dot{\Pi}$
$\displaystyle=-\frac{\Pi}{\tau_{R}}-\beta_{\Pi}\,\theta-\delta_{\Pi\Pi}\,\Pi\,\theta+\lambda_{\Pi\pi}\,\pi^{\mu\nu}\,\sigma_{\mu\nu},$
(91) $\displaystyle\dot{\pi}^{\langle\mu\nu\rangle}$
$\displaystyle=-\frac{\pi^{\mu\nu}}{\tau_{R}}+2\,\beta_{\pi}\,\sigma^{\mu\nu}+2\,\pi^{\langle\mu}_{\gamma}\,\omega^{\nu\rangle\gamma}-\tau_{\pi\pi}\,\pi^{\langle\mu}_{\gamma}\,\sigma^{\nu\rangle\gamma}$
$\displaystyle-\delta_{\pi\pi}\,\pi^{\mu\nu}\,\theta+\lambda_{\pi\Pi}\,\Pi\,\sigma^{\mu\nu}.$
(92)

The transport coefficients are given by

$\displaystyle\delta_{\Pi\Pi}$ $\displaystyle=-\left({\cal K}_{0}+{\cal
K}_{1}+{\cal K}_{2}\right)=-{\cal K}_{0}-\frac{5}{9}\chi,$ (93)
$\displaystyle\lambda_{\Pi\pi}$ $\displaystyle=\frac{2}{3}-{\cal K}_{0}+{\cal
K}_{3}=\frac{2}{3}-{\cal
K}_{0}+\frac{7\,\beta}{3\,\beta_{\pi}}\,I_{6,3}^{(3)},$ (94)
$\displaystyle\tau_{\pi\pi}$
$\displaystyle=2+\frac{4\,\beta}{\beta_{\pi}}\,I_{6,3}^{(3)},$ (95)
$\displaystyle\delta_{\pi\pi}$
$\displaystyle=\frac{5}{3}+\frac{7\beta}{3\beta_{\pi}}\,I_{6,3}^{(3)},$ (96)
$\displaystyle\lambda_{\pi\Pi}$ $\displaystyle=-\left({\cal K}_{4}+{\cal
K}_{5}\right)=-\frac{2}{3}\,\chi,\quad\text{with}$ (97) $\displaystyle\chi$
$\displaystyle\equiv\frac{3}{\beta_{\Pi}}\,\Big{[}{\cal
G}_{\beta}\left(I_{3,1}^{(0)}+I_{4,2}^{(1)}\right)-\sum_{i=1}^{3}{\cal
G}_{\alpha}^{i}\,\left(I_{3,1}^{(1),i}+I_{4,2}^{(2),i}\right)$
$\displaystyle\qquad\quad-\frac{\beta}{3}\left(5\,I_{5,2}^{(2)}+7I_{6,3}^{(3)}\right)\Bigr{]}.$
(98)

At vanishing chemical potential and ignoring quantum statistical effects, the
expressions for $\tau_{\pi\pi}$ and $\delta_{\pi\pi}$ are identical to those
found in Jaiswal:2014isa  for a single-species massive Boltzmann gas.
Moreover, in this limit the coefficient ${\cal K}_{0}=(\partial P/\partial
e)_{n}$ becomes equal to the squared speed of sound $c_{s}^{2}$, and the
coefficient $\chi$ reduces to313131For this, we substitute ${\cal
G}_{\beta}=-(e+P)/I_{3,0}^{(0)}=-\beta c_{s}^{2}$ and ${\cal G}_{\alpha}=0$ in
Eq. (98).

$\displaystyle\chi=\frac{\beta}{\beta_{\Pi}}\left[-3\,c_{s}^{2}\left(I_{3,1}^{(0)}+I_{4,2}^{(1)}\right)-5\,I_{5,2}^{(2)}-7\,I_{6,3}^{(3)}\right].$
(99)

Using Eqs. (8,9) in Ref. Jaiswal:2014isa , this expression for $\chi$ is found
to match with their Eq. (40). Accordingly, the transport coefficients obtained
in this section fully reduce to those obtained in Jaiswal:2014isa  for a
single-component Boltzmann gas at $\mu=0$ in the appropriate limit.

Also, in the conformal limit ($m=0$), the coefficient
$I_{6,3}^{(3)}\to-\beta_{\pi}/(7\beta)$. Accordingly, $\tau_{\pi\pi}\to 10/7$,
$\delta_{\pi\pi}\to 4/3$, such that Eq. (92) reduces to the shear stress
evolution equation derived in Jaiswal:2015mxa  for a conformal quark-gluon gas
at finite quark density.

##  Appendix B Second-order entropy current from conformal kinetic theory

We start with the general definitions (33,34) for the entropy current in
kinetic theory, assuming small deviations of the distribution functions from
their equilibrium forms: $f^{i}=f^{i}_{\mathrm{eq}}+\delta f^{i}$. One can
then write

$S^{\mu}=s_{\mathrm{eq}}u^{\mu}+\delta S^{\mu}$ (100)

where $s_{\mathrm{eq}}$ is the equilibrium entropy density of the system and
the deviation $\delta S^{\mu}$ can be Taylor expanded about equilibrium
as323232Different from App. A we here shall not include include any degeneracy
factors in the integration measure but simply write $dP\equiv
d^{3}p/[(2\pi^{3})\,p]$. Also, as all particles are massless, we use a common
four-momentum $p^{\mu}$.

$\displaystyle\delta
S^{\mu}={-}\sum_{i=1}^{3}g_{i}\\!\\!\int\\!\\!dP\,p^{\mu}\Big{(}\phi^{\prime}_{i}[f^{i}_{\mathrm{eq}}]\,\delta
f^{i}+\frac{\phi^{\prime\prime}_{i}[f^{i}_{\mathrm{eq}}]}{2}(\delta
f^{i})^{2}+\cdots\Big{)}.$

Here $\phi^{\prime}_{i}[f^{i}_{\mathrm{eq}}]$ and
$\phi^{\prime\prime}_{i}[f^{i}_{\mathrm{eq}}]$ denote the first and second
derivative of $\phi_{i}[f^{i}]$ w.r.t. $f^{i}$, evaluated at
$f^{i}=f^{i}_{\mathrm{eq}}$.

The leading-order term in $\delta S^{\mu}$ (i.e. the one linear in $\delta f$)
can be simplified using the Landau matching conditions. Taking the first
derivative of Eq. (34),

$\displaystyle\phi^{\prime}_{i}[f^{i}_{\mathrm{eq}}]=\ln\left(\frac{f^{i}_{\mathrm{eq}}}{1+\theta_{i}f^{i}_{\mathrm{eq}}}\right),$
(101)

and inserting the equilibrium distributions for quarks, anti-quarks, and
gluons yields $\phi^{\prime}_{i}[f^{i}_{\mathrm{eq}}]=-\beta\left(u\cdot
p\right)+\alpha_{i}$, where $\alpha_{q}=-\alpha_{\bar{q}}=\alpha$,
$\alpha_{g}=0$. With this the first term in $\delta S^{\mu}$ evaluates to

$\displaystyle\delta S^{\mu}_{1}$ $\displaystyle=\beta u_{\nu}\int
dP\,p^{\mu}\,p^{\nu}\,\left[g_{q}\left(\delta f^{q}+\delta
f^{\bar{q}}\right)+g_{g}\,\delta f^{g}\right]$ $\displaystyle\quad-\alpha\int
dP\,p^{\mu}\,g_{q}\left(\delta f^{q}-\delta f^{\bar{q}}\right)$ (102)
$\displaystyle=\beta\,u_{\nu}\delta
T^{\mu\nu}-\alpha\,n^{\mu}=-\alpha\,n^{\mu};$ (103)

here $n^{\mu}$ is the net-quark diffusion current, and we used the matching
condition $u_{\nu}\delta T^{\mu\nu}=0$. This result could have been
anticipated as to linear-order in dissipation the only available four-vector
is $n^{\mu}$ Muronga:2003ta . Moreover, this shows that the non-equilibrium
correction to the entropy density $u_{\mu}S^{\mu}$ is necessarily of second
order in the dissipative fluxes as $u_{\mu}n^{\mu}=0$. This is also expected,
as for a system with given values of energy and conserved charge densities the
thermal equilibrium state represents a maximum of the entropy density where
its first derivatives with respect to the dissipative flows vanish.

Let us now compute the corrections of ${\cal O}\bigl{(}(\delta f)^{2}\bigr{)}$
in the entropy four-current:

$\displaystyle
S^{\mu}=s_{\mathrm{eq}}\,u^{\mu}-\alpha\,n^{\mu}-\sum_{i=1}^{3}g_{i}\int
dP\,p^{\mu}\,\frac{\phi^{\prime\prime}_{i}[f^{i}_{\mathrm{eq}}]}{2}\,(\delta
f^{i})^{2}.$

With $\tilde{f}^{i}_{\mathrm{eq}}=1+\theta_{i}f^{i}_{\mathrm{eq}}$,
$\phi^{\prime\prime}_{i}[f^{i}_{\mathrm{eq}}]=1/(f^{i}_{\mathrm{eq}}\tilde{f}^{i}_{\mathrm{eq}})\equiv
1/{\cal F}_{\mathrm{eq}}^{i}$, and the ${\cal}(\delta f)^{2}$ correction to
the entropy current is

$\displaystyle\delta S^{\mu}_{2}=-\sum_{i}^{3}\frac{g_{i}}{2}\int
dP\,p^{\mu}\frac{(\delta f^{i})^{2}}{{\cal F}_{\mathrm{eq}}^{i}}.$ (104)

From Ref. Jaiswal:2015mxa  we know that to first order in the Chapman-Enskog
expansion

$\displaystyle\frac{\delta f^{q}}{{\cal F}_{\mathrm{eq}}^{q}}$
$\displaystyle={\cal A}_{p}^{\mu\nu}\pi_{\mu\nu}+{\cal
B}_{p}\,p^{\mu}\,n_{\mu},$ (105) $\displaystyle\frac{\delta f^{\bar{q}}}{{\cal
F}_{\mathrm{eq}}^{\bar{q}}}$ $\displaystyle={\cal
A}_{p}^{\mu\nu}\pi_{\mu\nu}+{\cal\bar{B}}_{p}\,p^{\mu}\,n_{\mu},$ (106)
$\displaystyle\frac{\delta f^{g}}{{\cal F}_{\mathrm{eq}}^{q}}$
$\displaystyle={\cal A}_{p}^{\mu\nu}\pi_{\mu\nu},$ (107)

with

$\displaystyle{\cal A}_{p}^{\mu\nu}$ $\displaystyle\equiv\frac{\beta}{2(u\cdot
p)\beta_{\pi}}p^{\mu}p^{\nu},\quad$ $\displaystyle{\cal B}_{p}$
$\displaystyle\equiv\frac{1}{\beta_{n}}\left(\frac{n}{e{+}P}-\frac{1}{u\cdot
p}\right),$ $\displaystyle{\cal\bar{B}}_{p}$
$\displaystyle\equiv\frac{1}{\beta_{n}}\left(\frac{n}{e{+}P}+\frac{1}{u\cdot
p}\right);$ (108)

here $\beta_{\pi}\equiv 4P/5$ and $\beta_{n}\equiv\frac{1}{3}(\partial
n/\partial\alpha)_{\beta}-n^{2}T/(4P)$. Using these definitions one obtains

$\displaystyle\sum_{i=1}^{3}g_{i}\frac{(\delta f^{i})^{2}}{{\cal
F}_{\mathrm{eq}}^{i}}$
$\displaystyle=\frac{\beta^{2}}{4\beta_{\pi}^{2}}\pi_{\alpha\beta}\,\pi_{\gamma\delta}\,p^{\alpha}p^{\beta}p^{\gamma}p^{\delta}\,\frac{1}{(u\cdot
p)^{2}}\,\left[g_{q}\left({\cal F}_{\mathrm{eq}}^{q}{+}{\cal
F}_{\mathrm{eq}}^{\bar{q}}\right)+g_{g}\,{\cal F}_{\mathrm{eq}}^{g}\right]$
$\displaystyle+\frac{1}{\beta_{n}^{2}}n_{\alpha}\,n_{\beta}\,p^{\alpha}\,p^{\beta}\,\,g_{q}\left[\left(\frac{n^{2}}{\left(e{+}P\right)^{2}}+\frac{1}{(u\cdot
p)^{2}}\right)\left({\cal F}_{\mathrm{eq}}^{q}{+}{\cal
F}_{\mathrm{eq}}^{\bar{q}}\right)-\frac{2n}{(e{+}P)}\frac{1}{u\cdot
p}\left({\cal F}_{\mathrm{eq}}^{q}{-}{\cal
F}_{\mathrm{eq}}^{\bar{q}}\right)\right]$
$\displaystyle+\frac{\beta}{\beta_{\pi}\beta_{n}}\pi_{\alpha\beta}\,n_{\gamma}\,p^{\alpha}p^{\beta}p^{\gamma}\frac{1}{u\cdot
p}\,g_{q}\left[\frac{n}{e{+}P}\left({\cal F}_{\mathrm{eq}}^{q}{+}{\cal
F}_{\mathrm{eq}}^{\bar{q}}\right)-\frac{1}{u\cdot p}\left({\cal
F}_{\mathrm{eq}}^{q}{-}{\cal F}_{\mathrm{eq}}^{\bar{q}}\right)\right].$ (109)

To obtain $\delta S^{\mu}_{2}$ we will substitute the above expression in Eq.
(104). The term quadratic in the shear stress yields

$\displaystyle\delta S^{\mu}_{2,\pi\pi}$
$\displaystyle=-\frac{\beta^{2}}{8\,\beta_{\pi}^{2}}\,\pi_{\alpha\beta}\,\pi_{\gamma\delta}\,\int\frac{dP}{(u\cdot
p)^{2}}\,p^{\mu}p^{\alpha}p^{\beta}p^{\gamma}p^{\delta}\,$
$\displaystyle\qquad\qquad\qquad\qquad\times\left[g_{q}\left({\cal
F}_{\mathrm{eq}}^{q}{+}{\cal F}_{\mathrm{eq}}^{\bar{q}}\right)+g_{g}\,{\cal
F}_{\mathrm{eq}}^{g}\right]$
$\displaystyle=-\frac{\beta^{2}}{4\,\beta_{\pi}^{2}}u^{\mu}\,\pi^{\alpha\beta}\,\pi_{\alpha\beta}\,\left[g_{q}\left(J_{3,2}^{q}{+}J_{3,2}^{\bar{q}}\right)+g_{g}J_{3,2}^{g}\right],$
(110)

where we introduced the thermodynamic integrals

$\displaystyle J_{n,l}^{i}\equiv\frac{1}{\left(2l+1\right)!!}\,\int
dP\,\left(u\cdot p\right)^{n}\,{\cal F}_{\mathrm{eq}}^{i}.$ (111)

Noting that

$\displaystyle{\cal
F}_{\mathrm{eq}}^{i}=f_{\mathrm{eq}}^{i}\,\tilde{f}_{\mathrm{eq}}^{i}=-\frac{1}{u\cdot
p}\,\frac{\partial f_{\mathrm{eq}}^{i}}{\partial\beta}$ (112)

(with the derivative to be taken at fixed $\alpha$) we simplify

$\displaystyle
g_{q}\left(J^{q}_{3,2}+J^{\bar{q}}_{3,2}\right)+g_{g}J^{g}_{3,2}$
$\displaystyle=-\frac{1}{15}\sum_{i}g_{i}\int dP\,\left(u\cdot
p\right)^{2}\frac{\partial f_{\mathrm{eq}}^{i}}{\partial\beta}$
$\displaystyle=-\frac{1}{15}\frac{\partial
e_{\mathrm{eq}}}{\partial\beta}=\frac{\beta_{\pi}}{\beta}$ (113)

such that

$\displaystyle\delta
S^{\mu}_{2,\pi\pi}=-\frac{\beta}{4\beta_{\pi}}\,u^{\mu}\,\pi^{\alpha\beta}\,\pi_{\alpha\beta}.$
(114)

The term in Eq. (109) quadratic in the quark diffusion current is manipulated
similarly:

$\displaystyle\delta
S^{\mu}_{2,nn}=-\frac{g_{q}}{2\beta_{n}^{2}}\,n_{\alpha}n_{\beta}\int
dP\,p^{\mu}p^{\alpha}p^{\beta}\bigg{[}-\frac{2n\,({\cal
F}_{\mathrm{eq}}^{q}{-}{\cal F}_{\mathrm{eq}}^{\bar{q}})}{(e{+}P)(u\cdot
p)}\,$
$\displaystyle\qquad\qquad\qquad+\bigg{(}\frac{n^{2}}{\left(e{+}P\right)^{2}}+\frac{1}{(u\cdot
p)^{2}}\bigg{)}\left({\cal F}_{\mathrm{eq}}^{q}{+}{\cal
F}_{\mathrm{eq}}^{\bar{q}}\right)\bigg{]}$
$\displaystyle=\frac{g_{q}}{2\beta_{n}^{2}}\,u^{\mu}\,n^{\alpha}n_{\alpha}\bigg{[}\frac{n^{2}}{(e{+}P)^{2}}\,\left(J_{3,1}^{q}{+}J_{3,1}^{\bar{q}}\right)+\left(J_{1,1}^{q}{+}J_{1,1}^{\bar{q}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{2n}{e{+}P}\left(J_{2,1}^{q}{-}J_{2,1}^{\bar{q}}\right)\bigg{]}.$
(115)

With Eq. (112) this simplifies further to

$\displaystyle g_{q}\left(J_{3,1}^{q}+J_{3,1}^{\bar{q}}\right)$
$\displaystyle=-\frac{1}{3}\frac{\partial}{\partial\beta}\int dP\,\left(u\cdot
p\right)^{2}g_{q}\left(f^{q}_{\mathrm{eq}}+f^{\bar{q}}_{\mathrm{eq}}\right)$
$\displaystyle=-\frac{1}{3}\frac{\partial\left(e^{q}_{\mathrm{eq}}+e^{\bar{q}}_{\mathrm{eq}}\right)}{\partial\beta}.$
(116)

To compute $g_{q}(J_{1,1}^{q}+J_{1,1}^{\bar{q}})$ we use ${\cal
F}_{\mathrm{eq}}^{q}=(\partial f_{\mathrm{eq}}^{q}/\partial\alpha)$, ${\cal
F}_{\mathrm{eq}}^{\bar{q}}=-(\partial
f_{\mathrm{eq}}^{\bar{q}}/\partial\alpha)$ and obtain

$\displaystyle g_{q}\left(J_{1,1}^{q}{+}J_{1,1}^{\bar{q}}\right)$
$\displaystyle=\frac{1}{3}\,\frac{\partial}{\partial\alpha}\int
dP\,(u{\,\cdot\,}p)\,g_{q}\left(f^{q}_{\mathrm{eq}}{-}f^{\bar{q}}_{\mathrm{eq}}\right)=\frac{1}{3}\frac{\partial
n}{\partial\alpha}.$

Finally we calculate

$\displaystyle g_{q}\left(J_{2,1}^{q}-J_{2,1}^{\bar{q}}\right)$
$\displaystyle=\frac{1}{3}\,\frac{\partial}{\partial\alpha}\int
dP\,\left(u\cdot
p\right)^{2}\,g_{q}\left(f_{\mathrm{eq}}^{q}+f_{\mathrm{eq}}^{\bar{q}}\right)$
$\displaystyle=\frac{1}{3}\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\alpha}.$
(117)

All this allows to express $\delta S^{\mu}_{2,nn}$ in compact form:

$\displaystyle\delta S^{\mu}_{2,nn}$
$\displaystyle=\frac{1}{6\beta_{n}^{2}}\,u^{\mu}\,n^{\alpha}\,n_{\alpha}\bigg{[}-\frac{n^{2}}{(e+P)^{2}}\frac{\partial\left(e^{q}_{\mathrm{eq}}+e^{\bar{q}}_{\mathrm{eq}}\right)}{\partial\beta}+\frac{\partial
n}{\partial\alpha}$
$\displaystyle-\frac{2n}{(e+P)}\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\alpha}\bigg{]}.$
(118)

Finally, we evaluate the cross term between quark diffusion and shear stress
in Eq. (109):

Figure 13: Upper panels compare phase trajectories from kinetic theory (blue
and black) and second-order hydrodynamics (orange and magenta). Panels (a) and
(b) correspond, respectively, to initial conditionsets I and II. The lower
panels compare the evolution of the shear (panel c) and bulk (panel d) inverse
Reynolds numbers for set II of initial conditions.

$\displaystyle\delta
S^{\mu}_{2,n\pi}=-\frac{1}{2}\frac{\beta}{\beta_{\pi}\beta_{n}}\,\pi_{\alpha\beta}\,n_{\gamma}\int\frac{dP}{u\cdot
p}\,p^{\mu}p^{\alpha}p^{\beta}p^{\gamma}$
$\displaystyle\times\Bigl{[}\frac{n}{e{+}P}\left({\cal
F}_{\mathrm{eq}}^{q}{+}{\cal F}_{\mathrm{eq}}^{\bar{q}}\right)-\frac{1}{u\cdot
p}\left({\cal F}^{q}{-}{\cal F}_{\mathrm{eq}}^{\bar{q}}\right)\Bigr{]}$
$\displaystyle=-\frac{\beta}{\beta_{\pi}\beta_{n}}\pi^{\mu\alpha}\,n_{\alpha}\,g_{q}\Bigl{[}\frac{n}{e{+}P}\left(J_{3,2}^{q}{+}J_{3,2}^{\bar{q}}\right)-\left(J_{2,2}^{q}{-}J_{2,2}^{\bar{q}}\right)\Bigr{]}.$

Using

$\displaystyle g_{q}\left(J_{3,2}^{q}+J_{3,2}^{\bar{q}}\right)$
$\displaystyle=-\frac{1}{15}\,\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\beta},$
$\displaystyle g_{q}\left(J_{2,2}^{q}-J_{2,2}^{\bar{q}}\right)$
$\displaystyle=\frac{1}{15}\frac{\partial}{\partial\alpha}\int\,dP\,\left(u\cdot
p\right)^{2}\,g_{q}\left(f_{\mathrm{eq}}^{q}+f_{\mathrm{eq}}^{\bar{q}}\right)$
$\displaystyle=\frac{1}{15}\,\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\alpha}$
(119)

we write

$\displaystyle\delta S^{\mu}_{2,n\pi}$
$\displaystyle=\frac{\beta}{15\,\beta_{\pi}\beta_{n}}\,\pi^{\mu\alpha}\,n_{\alpha}\,\bigg{[}\frac{n}{e+P}\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\beta}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\alpha}\bigg{]}.$
(120)

Collecting all contributions we finally obtain for the second-order entropy
current of a conformal gas of quarks and gluons at finite chemical potential

$\displaystyle S^{\mu}$
$\displaystyle=s_{\mathrm{eq}}\,u^{\mu}-\alpha\,n^{\mu}-\frac{\beta}{4\beta_{\pi}}\,u^{\mu}\,\pi^{\alpha\beta}\,\pi_{\alpha\beta}$
$\displaystyle+c_{nn}\,u^{\mu}\,n^{\alpha}\,n_{\alpha}+c_{n\pi}\,\pi^{\mu\alpha}\,n_{\alpha},$
(121)

with coefficients $c_{nn}$ and $c_{n\pi}$ given by

$\displaystyle c_{nn}$
$\displaystyle=\frac{1}{6\beta_{n}^{2}}\bigg{[}-\left(\frac{n}{e+P}\right)^{2}\frac{\partial\left(e^{q}_{\mathrm{eq}}+e^{\bar{q}}_{\mathrm{eq}}\right)}{\partial\beta}+\frac{\partial
n}{\partial\alpha}$
$\displaystyle\qquad\qquad-\frac{2n}{e+P}\,\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\alpha}\biggr{]},$
(122) $\displaystyle c_{n\pi}$
$\displaystyle=\frac{\beta}{15\beta_{\pi}\beta_{n}}\biggl{[}\frac{n}{e+P}\,\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\beta}+\frac{\partial\left(e_{\mathrm{eq}}^{q}+e_{\mathrm{eq}}^{\bar{q}}\right)}{\partial\alpha}\biggr{]}.$

The last term in Eq. (B) implies that the entropy flux in the fluid rest frame
is not along the direction of net-quark diffusion unless the latter points
along an eigen-direction of the shear stress tensor.333333Note that at
vanishing $\mu_{B}$ the second-order entropy density,
$s=s_{\mathrm{eq}}-\beta/(4\beta_{\pi})\pi^{\alpha\beta}\pi_{\alpha\beta}$,
with $\beta_{\pi}=4P/5$, is identical in form to that obtained for a conformal
Boltzmann gas Chattopadhyay:2014lya .

##  Appendix C Non-conformal kinetic theory versus hydrodynamics

In this Appendix we compare some of the results obtained in Sec. V using non-
conformal kinetic theory (Sec. V.2) and second-order hydrodynamics (Sec. V.1).
The insights gained from this comparison are not new but confirm the findings
reported in Refs. Jaiswal:2021uvv ; Chattopadhyay:2021ive  and extend them to
systems with non-zero net baryon charge.

Figure 13 compares the phase trajectories (upper panels) and the evolution of
the dissipative flows (lower panels) from second-order viscous hydrodynamics
(“hydro”, taken from Fig. 8) and kinetic theory (“KT”, taken from Fig. 10).
While the hydrodynamic and kinetic expansion trajectories agree well for
equilibrium initial conditions (where dissipative effects manifest themselves
by viscous heating), and this agreement gets better with increasing
interaction strength, they differ substantially for the far-off-equilibrium
initial conditions for which the expansion initially leads to viscous cooling.
Clearly, the hydrodynamic approach as a macroscopic approximation of the
underlying microscopic kinetic theory degrades in far-off-equilibrium
situations. This might not have been expected from the evolutions of the
normalized shear viscous stress shown in panel (c) which shows discrepancies
between microscopic and macroscopic approaches that are of similar size for
both sets of initial conditions. On the other hand, the evolution of the
normalized bulk viscous pressure shown in panel (d) also shows much larger
discrepancies between “hydro” and “KT” for the far-off-equilibrium initial
conditions, in particular in the amplitude of its oscillations around zero at
early times. Both features are shared by both weakly and strongly coupled
fluids, at least qualitatively.

Figure 14: Evolution of the scaled longitudinal pressure in kinetic theory
(blue curves) and second-order hydrodynamics (magenta curves). The kinetic
theory trajectories are identical to those in Fig. 12.

Figure 14 illustrates the observation made in Refs. Jaiswal:2021uvv ;
Chattopadhyay:2021ive  that the failure of the (standard) second-order viscous
hydrodynamic approximation is most severe at early times when the expansion
rate is large and the system is far away from equilibrium. During this stage
the expansion is dominated by free-streaming and the exact kinetic evolution
(blue) of the longitudinal pressure is ruled by a “far-from-equilibrium
attractor” that is not reproduced in the hydrodynamic approximation (magenta)
Chattopadhyay:2021ive . The discrepancies are particularly large for far-from-
equilibrium initial condition ($(P/P_{L})_{0}\gg 1$ or $(P/P_{L})_{0}\ll
1$).343434In the latter case the hydrodynamic approximation even produces a
negative longitudinal pressure at early times, violating a fundamental kinetic
theory limit for the system at hand Chattopadhyay:2021ive . Readers interested
in a more detailed discussion are referred to Refs. Jaiswal:2021uvv ;
Chattopadhyay:2021ive .

## References

  * (1) M. A. Stephanov, QCD phase diagram: An Overview, PoS LAT2006 (2006) 024.  arXiv:hep-lat/0701002, doi:10.22323/1.032.0024.
  * (2) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, K. K. Szabo, The Order of the quantum chromodynamics transition predicted by the standard model of particle physics, Nature 443 (2006) 675–678.  arXiv:hep-lat/0611014, doi:10.1038/nature05120.
  * (3) A. Bazavov, et al., Equation of state and QCD transition at finite temperature, Phys. Rev. D 80 (2009) 014504.  arXiv:0903.4379, doi:10.1103/PhysRevD.80.014504.
  * (4) S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, K. K. Szabo, The QCD equation of state with dynamical quarks, JHEP 11 (2010) 077.  arXiv:1007.2580, doi:10.1007/JHEP11(2010)077.
  * (5) A. Bazavov, et al., The chiral and deconfinement aspects of the QCD transition, Phys. Rev. D 85 (2012) 054503.  arXiv:1111.1710, doi:10.1103/PhysRevD.85.054503.
  * (6) M. A. Stephanov, K. Rajagopal, E. V. Shuryak, Signatures of the tricritical point in QCD, Phys. Rev. Lett. 81 (1998) 4816–4819.  arXiv:hep-ph/9806219, doi:10.1103/PhysRevLett.81.4816.
  * (7) M. A. Stephanov, K. Rajagopal, E. V. Shuryak, Event-by-event fluctuations in heavy ion collisions and the QCD critical point, Phys. Rev. D 60 (1999) 114028\.  arXiv:hep-ph/9903292, doi:10.1103/PhysRevD.60.114028.
  * (8) W. Busza, K. Rajagopal, W. van der Schee, Heavy Ion Collisions: The Big Picture, and the Big Questions, Ann. Rev. Nucl. Part. Sci. 68 (2018) 339–376.  arXiv:1802.04801, doi:10.1146/annurev-nucl-101917-020852.
  * (9) A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, N. Xu, Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan, Phys. Rept. 853 (2020) 1–87.  arXiv:1906.00936, doi:10.1016/j.physrep.2020.01.005.
  * (10) P. C. Hohenberg, B. I. Halperin, Theory of dynamic critical phenomena, Reviews of Modern Physics 49 (3) (1977) 435.
  * (11) Y. Hatta, M. A. Stephanov, Proton number fluctuation as a signal of the QCD critical endpoint, Phys. Rev. Lett. 91 (2003) 102003, [Erratum: Phys.Rev.Lett. 91, 129901 (2003)].  arXiv:hep-ph/0302002, doi:10.1103/PhysRevLett.91.102003.
  * (12) M. A. Stephanov, Non-Gaussian fluctuations near the QCD critical point, Phys. Rev. Lett. 102 (2009) 032301.  arXiv:0809.3450, doi:10.1103/PhysRevLett.102.032301.
  * (13) M. A. Stephanov, On the sign of kurtosis near the QCD critical point, Phys. Rev. Lett. 107 (2011) 052301.  arXiv:1104.1627, doi:10.1103/PhysRevLett.107.052301.
  * (14) X. Luo, N. Xu, Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview, Nucl. Sci. Tech. 28 (8) (2017) 112.  arXiv:1701.02105, doi:10.1007/s41365-017-0257-0.
  * (15) L. Adamczyk, et al., Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program, Phys. Rev. C 96 (4) (2017) 044904.  arXiv:1701.07065, doi:10.1103/PhysRevC.96.044904.
  * (16) M. Nahrgang, S. Leupold, C. Herold, M. Bleicher, Nonequilibrium chiral fluid dynamics including dissipation and noise, Phys. Rev. C84 (2011) 024912.  arXiv:1105.0622, doi:10.1103/PhysRevC.84.024912.
  * (17) M. Nahrgang, S. Leupold, M. Bleicher, Equilibration and relaxation times at the chiral phase transition including reheating, Phys. Lett. B711 (2012) 109–116.  arXiv:1105.1396, doi:10.1016/j.physletb.2012.03.059.
  * (18) S. Mukherjee, R. Venugopalan, Y. Yin, Real time evolution of non-Gaussian cumulants in the QCD critical regime, Phys. Rev. C92 (3) (2015) 034912.  arXiv:1506.00645, doi:10.1103/PhysRevC.92.034912.
  * (19) C. Herold, M. Nahrgang, Y. Yan, C. Kobdaj, Dynamical net-proton fluctuations near a QCD critical point, Phys. Rev. C93 (2) (2016) 021902.  arXiv:1601.04839, doi:10.1103/PhysRevC.93.021902.
  * (20) M. Stephanov, Y. Yin, Hydrodynamics with parametric slowing down and fluctuations near the critical point, Phys. Rev. D98 (3) (2018) 036006.  arXiv:1712.10305, doi:10.1103/PhysRevD.98.036006.
  * (21) T. Hirano, R. Kurita, K. Murase, Hydrodynamic fluctuations of entropy in one-dimensionally expanding system, Nucl. Phys. A984 (2019) 44–67.  arXiv:1809.04773, doi:10.1016/j.nuclphysa.2019.01.010.
  * (22) M. Singh, C. Shen, S. McDonald, S. Jeon, C. Gale, Hydrodynamic fluctuations in relativistic heavy-ion collisions, Nucl. Phys. A982 (2019) 319–322.  arXiv:1807.05451, doi:10.1016/j.nuclphysa.2018.10.061.
  * (23) M. Nahrgang, M. Bluhm, T. Schäfer, S. A. Bass, Diffusive dynamics of critical fluctuations near the QCD critical point, Phys. Rev. D99 (11) (2019) 116015.  arXiv:1804.05728, doi:10.1103/PhysRevD.99.116015.
  * (24) Y. Yin, The QCD critical point hunt: emergent new ideas and new dynamicsarXiv:1811.06519.
  * (25) X. An, et al., The BEST framework for the search for the QCD critical point and the chiral magnetic effect, Nucl. Phys. A 1017 (2022) 122343.  arXiv:2108.13867, doi:10.1016/j.nuclphysa.2021.122343.
  * (26) S. Gupta, D. Mallick, D. K. Mishra, B. Mohanty, N. Xu, Limits of thermalization in relativistic heavy ion collisions, Phys. Lett. B 829 (2022) 137021.  doi:10.1016/j.physletb.2022.137021.
  * (27) J. Steinheimer, M. Bleicher, H. Petersen, S. Schramm, H. Stocker, D. Zschiesche, (3+1)-dimensional hydrodynamic expansion with a critical point from realistic initial conditions, Phys. Rev. C 77 (2008) 034901.  arXiv:0710.0332, doi:10.1103/PhysRevC.77.034901.
  * (28) L. Du, X. An, U. Heinz, Baryon transport and the QCD critical point, Phys. Rev. C 104 (6) (2021) 064904.  arXiv:2107.02302, doi:10.1103/PhysRevC.104.064904.
  * (29) C. Shen, Studying QGP with flow: A theory overview, Nucl. Phys. A 1005 (2021) 121788\.  arXiv:2001.11858, doi:10.1016/j.nuclphysa.2020.121788.
  * (30) U. Heinz, K. S. Lee, M. J. Rhoades-Brown, $K^{+}$ and $K^{-}$ slope parameters as a signature for deconfinement at finite baryon density, Phys. Rev. Lett. 58 (1987) 2292–2295.  doi:10.1103/PhysRevLett.58.2292.
  * (31) K.-S. Lee, U. Heinz, The phase structure of strange matter, Phys. Rev. D 47 (1993) 2068–2080.  doi:10.1103/PhysRevD.47.2068.
  * (32) S. J. Cho, K. S. Lee, U. Heinz, Strange matter lumps in the early universe, Phys. Rev. D 50 (1994) 4771–4780.  doi:10.1103/PhysRevD.50.4771.
  * (33) C. M. Hung, E. V. Shuryak, Equation of state, radial flow and freezeout in high-energy heavy ion collisions, Phys. Rev. C 57 (1998) 1891–1906.  arXiv:hep-ph/9709264, doi:10.1103/PhysRevC.57.1891.
  * (34) T. Dore, J. Noronha-Hostler, E. McLaughlin, Far-from-equilibrium search for the QCD critical point, Phys. Rev. D 102 (7) (2020) 074017.  arXiv:2007.15083, doi:10.1103/PhysRevD.102.074017.
  * (35) T. Dore, J. M. Karthein, I. Long, D. Mroczek, J. Noronha-Hostler, P. Parotto, C. Ratti, Y. Yamauchi, Critical lensing and kurtosis near a critical point in the QCD phase diagram in and out-of-equilibriumarXiv:2207.04086.
  * (36) P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, M. Stephanov, QCD equation of state matched to lattice data and exhibiting a critical point singularity, Phys. Rev. C 101 (3) (2020) 034901.  arXiv:1805.05249, doi:10.1103/PhysRevC.101.034901.
  * (37) P. Romatschke, Relativistic (Lattice) Boltzmann Equation with Non-Ideal Equation of State, Phys. Rev. D 85 (2012) 065012.  arXiv:1108.5561, doi:10.1103/PhysRevD.85.065012.
  * (38) G. S. Denicol, H. Niemi, E. Molnar, D. H. Rischke, Derivation of transient relativistic fluid dynamics from the Boltzmann equation, Phys. Rev. D85 (2012) 114047, [Erratum: Phys. Rev. D91, 039902 (2015)].  arXiv:1202.4551, doi:10.1103/PhysRevD.85.114047,10.1103/PhysRevD.91.039902.
  * (39) W. Florkowski, R. Ryblewski, M. Strickland, Testing viscous and anisotropic hydrodynamics in an exactly solvable case, Phys. Rev. C88 (2013) 024903.  arXiv:1305.7234, doi:10.1103/PhysRevC.88.024903.
  * (40) A. Jaiswal, Relativistic dissipative hydrodynamics from kinetic theory in relaxation time approximation, Phys. Rev. C87 (2013) 051901.  arXiv:1302.6311, doi:10.1103/PhysRevC.87.051901.
  * (41) P. Romatschke, Relativistic hydrodynamic attractors with broken symmetries: non-conformal and non-homogeneous, JHEP 12 (2017) 079.  arXiv:1710.03234, doi:10.1007/JHEP12(2017)079.
  * (42) D. Teaney, L. Yan, Second order viscous corrections to the harmonic spectrum in heavy ion collisions, Phys. Rev. C 89 (1) (2014) 014901.  arXiv:1304.3753, doi:10.1103/PhysRevC.89.014901.
  * (43) G. S. Rocha, G. S. Denicol, J. Noronha, Novel Relaxation Time Approximation to the Relativistic Boltzmann Equation, Phys. Rev. Lett. 127 (4) (2021) 042301.  arXiv:2103.07489, doi:10.1103/PhysRevLett.127.042301.
  * (44) G. S. Rocha, G. S. Denicol, Transient fluid dynamics with general matching conditions: A first study from the method of moments, Phys. Rev. D 104 (9) (2021) 096016.  arXiv:2108.02187, doi:10.1103/PhysRevD.104.096016.
  * (45) D. Dash, S. Bhadury, S. Jaiswal, A. Jaiswal, Extended relaxation time approximation and relativistic dissipative hydrodynamics, Phys. Lett. B 831 (2022) 137202.  arXiv:2112.14581, doi:10.1016/j.physletb.2022.137202.
  * (46) S. R. De Groot, W. A. Van Leeuwen, C. G. Van Weert, Relativistic Kinetic Theory. Principles and Applications, North-Holland Publishing Company, 1980\.
  * (47) G. S. Denicol, S. Jeon, C. Gale, Transport Coefficients of Bulk Viscous Pressure in the 14-moment approximation, Phys. Rev. C 90 (2) (2014) 024912.  arXiv:1403.0962, doi:10.1103/PhysRevC.90.024912.
  * (48) J. D. Bjorken, Highly relativistic nucleus-nucleus collisions: the central rapidity region, Phys. Rev. D 27 (1983) 140–151.  doi:10.1103/PhysRevD.27.140.
  * (49) A. Jaiswal, R. Ryblewski, M. Strickland, Transport coefficients for bulk viscous evolution in the relaxation time approximation, Phys. Rev. C90 (2014) 044908.  arXiv:1407.7231, doi:10.1103/PhysRevC.90.044908.
  * (50) M. P. Heller, M. Spalinski, Hydrodynamics beyond the gradient expansion: resurgence and resummation, Phys. Rev. Lett. 115 (2015) 072501.  arXiv:1503.07514, doi:10.1103/PhysRevLett.115.072501.
  * (51) P. Romatschke, Relativistic fluid dynamics far from local equilibrium, Phys. Rev. Lett. 120 (2018) 012301.  arXiv:1704.08699, doi:10.1103/PhysRevLett.120.012301.
  * (52) S. Jaiswal, C. Chattopadhyay, A. Jaiswal, S. Pal, U. Heinz, Exact solutions and attractors of higher-order viscous fluid dynamics for Bjorken flow, Phys. Rev. C 100 (3) (2019) 034901.  arXiv:1907.07965, doi:10.1103/PhysRevC.100.034901.
  * (53) A. Kurkela, W. van der Schee, U. A. Wiedemann, B. Wu, Early- and Late-Time Behavior of Attractors in Heavy-Ion Collisions, Phys. Rev. Lett. 124 (10) (2020) 102301.  arXiv:1907.08101, doi:10.1103/PhysRevLett.124.102301.
  * (54) G. Baym, Thermal equilibration in ultrarelativistic heavy ion collisions, Phys. Lett. B 138 (1984) 18–22.  doi:10.1016/0370-2693(84)91863-X.
  * (55) P. Romatschke, M. Strickland, Collective modes of an anisotropic quark gluon plasma, Phys. Rev. D 68 (2003) 036004.  arXiv:hep-ph/0304092, doi:10.1103/PhysRevD.68.036004.
  * (56) S. Jaiswal, C. Chattopadhyay, L. Du, U. Heinz, S. Pal, Nonconformal kinetic theory and hydrodynamics for Bjorken flow, Phys. Rev. C 105 (2022) 024911.  arXiv:2107.10248, doi:10.1103/PhysRevC.105.024911.
  * (57) C. Chattopadhyay, U. Heinz, S. Pal, G. Vujanovic, Higher order and anisotropic hydrodynamics for Bjorken and Gubser flows, Phys. Rev. C97 (2018) 064909.  arXiv:1801.07755, doi:10.1103/PhysRevC.97.064909.
  * (58) S. Jaiswal, J.-P. Blaizot, R. S. Bhalerao, Z. Chen, A. Jaiswal, L. Yan, From moments of the distribution function to hydrodynamics $-$ the non-conformal casearXiv:2208.02750.