# Exploring freeze out and flow using exact solutions of conformal
hydrodynamics

Owen Bradley  Natural Science Division, Pepperdine University, Malibu, CA
90263, USA  Christopher Plumberg<EMAIL_ADDRESS>Natural
Science Division, Pepperdine University, Malibu, CA 90263, USA

###### Abstract

Exact solutions to the equations of hydrodynamics provide valuable benchmark
tests for numerical hydrodynamic codes and also provide useful insights into
the nature of hydrodynamic flow. In this paper, we introduce two novel,
closely related exact solutions with non-trivial rapidity dependence which are
generalizations of the well-known Gubser flow solution to conformal
hydrodynamics. We then use one of our solutions to explore the consequences of
choosing between two different criteria for implementing the freeze out
process in fluid dynamical simulations of nuclear collisions: freeze out at
constant temperature vs. freeze out at constant Knudsen number. We find that,
employing our exact solution, the differences between these freeze out
criteria are heavily influenced by the presence of strong collective flow. Our
results highlight the importance of accurately describing the freeze out
process in collisions with large flow gradients, particularly in small
systems.

##  I Introduction

Relativistic fluid dynamics has proved to be an extremely powerful tool for
the description of high-energy nuclear collisions Gale  _et al._ (2013); Jeon
and Heinz (2015); Romatschke and Romatschke (2019); Nagle and Zajc (2018);
Adolfsson  _et al._ (2020). To date, hydrodynamic models have been used to
quantitatively describe a host of experimental observables Song  _et al._
(2011); Bozek and Wyskiel-Piekarska (2012); Gardim  _et al._ (2012); Bozek and
Broniowski (2013); Niemi  _et al._ (2016a); Ryu  _et al._ (2015); McDonald
_et al._ (2017); Bernhard  _et al._ (2016); Gardim  _et al._ (2017); Alba  _et
al._ (2018); Giacalone  _et al._ (2018); Eskola  _et al._ (2018); Weller and
Romatschke (2017); Schenke  _et al._ (2019); Everett  _et al._ (2021). With
its ability now to make predictions which are accurate at the percent level
Noronha-Hostler  _et al._ (2016); Niemi  _et al._ (2016b); Adam  _et al._
(2016), the fluid dynamical paradigm has arguably begun to enter the domain of
precision physics.

One of the remaining systematic uncertainties affecting nuclear collision
phenomenology involves the treatment of the freeze out process which
terminates the fluid dynamical evolution. The most common approach is to
freeze out along a hypersurface of fixed temperature Heinz (2004);
Schnedermann  _et al._ (1993) or energy density Petersen  _et al._ (2008),
where the latter allows a straightforward generalization to finite chemical
potentials which are probed by nuclear collisions Denicol  _et al._ (2018);
Plumberg  _et al._ (2023a). An alternative (although less widely used)
criterion for imposing freeze out is to terminate at constant Knudsen number
$\mathrm{Kn}$ Jeon and Heinz (2015). The Knudsen number represents the
competition between mechanisms which drive the system out of equilibrium
(e.g., strong flow gradients) and those which relax the system toward local
equilibrium (e.g., strong interactions or dissipative effects) and thus
provides a measure of how well a system may be described fluid dynamically.
Freeze out at constant Knudsen number can thus be viewed as a choice to
decouple the system in terms of the degree to which the fluid dynamical
description is justifiable at a given stage in the evolution.

The extent to which the use of different freeze-out criteria influences
experimental observables has already received some attention in the context of
hydrodynamic models Schnedermann  _et al._ (1993); Eskola  _et al._ (2008);
Holopainen and Huovinen (2012); Niemi and Denicol (2014); Ahmad  _et al._
(2017); Huovinen and Holopainen (2017), including systems both with and
without event-by-event fluctuations. However, less work has been done on the
differences between these criteria in small systems where large flow gradients
drive fluid dynamics to its breaking point Niemi and Denicol (2014); Heinz and
Moreland (2019); Adolfsson  _et al._ (2020); Summerfield  _et al._ (2021);
Plumberg  _et al._ (2023b), and relatively little progress has been made in
understanding the quantitative implications of these differences for
observables which probe specifically the space-time geometry of nuclear
collisions, such as Hanbury Brown – Twiss interferometry Plumberg (2020a).

In this paper, we begin to address this gap by exploring the extent to which
the shape of the freeze out hypersurface is impacted by the choice of freeze
out criterion and the presence of collective flow in the system. In order to
build intuition and obtain results which can be straightforwardly interpreted,
we have chosen to answer this question in the context of an exact solution to
the equations of relativistic fluid dynamics. The search for such exact
solutions has a long history. Known exact solutions include those obtained by
Bjorken assuming longitudinal boost invariance Bjorken (1983) and a later
generalization by Gubser Gubser (2010); Gubser and Yarom (2011) to include
expansion in the transverse plane as well. Exact solutions have also been
constructed for rotating systems Hatta  _et al._ (2014a, b) and systems which
possess non-trivial dependence on the space-time rapidity Csörgő  _et al._
(2018); Shi  _et al._ (2022), as well as for various generalizations including
viscous effects Marrochio  _et al._ (2015).

In this study, we are interested in exploring the freeze out process using
exact solutions which exhibit strong transverse and longitudinal flow
simultaneously. To do this, we first construct two new, exact, Gubser-like
solutions which contain non-trivial transverse and longitudinal flow, and then
use one of these solutions to compare the geometries of fixed temperature and
fixed Knudsen number hypersurfaces. The first solution (“Solution I")
represents a generalization of Gubser’s ideal hydrodynamic solution and does
not include viscosity, while the second solution (“Solution II") is a
generalization of Gubser’s viscous solution Gubser (2010); Gubser and Yarom
(2011) (and thus allows us to define and study freeze out in terms of the
Knudsen number). Both of our solutions make use of a trick, recently suggested
in Shi  _et al._ (2022) and which we refer to here as a “temporal shift," to
induce a non-trivial rapidity dependence into the collective flow of a known
solution. In this paper we focus our attention on the derivation of the
solutions themselves, as well as their qualitative implications for the
geometry of the freeze out hypersurface. We defer a more thorough discussion
of the quantitative consequences for nuclear collision phenomenology, such as
HBT measurements, to a subsequent study.

We have organized the paper as follows. In Sec. II, we construct Solution I as
a generalization of Gubser’s ideal flow solution and show how to introduce
rapidity dependence by performing a temporal shift. We also discuss a number
of properties of Solution I and consider how its space-time dependence is
influenced by different choices of the free parameters it contains. In Sec.
III, we apply the same temporal shifting procedure to the viscous version of
Gubser’s original solution, in order to obtain our Solution II which possesses
both viscosity and non-trivial $\eta$-dependence. Finally, in Sec. IV, we
apply our Solution II to the question of how the presence of both transverse
and longitudinal flow influence the shape of the freeze out surface. Some
conclusions are presented in Sec. V.

##  II Solution I

Before presenting our derivation of Solution I, we first briefly review the
original ideal solution obtained by Gubser Gubser (2010); Gubser and Yarom
(2011).

###  II.1 The ideal Gubser solution

Relativistic fluid dynamics is predicated on the covariant conservation of
energy and momentum, expressed by

$\displaystyle\nabla_{\mu}T^{\mu\nu}=0$ (1)

where $\nabla_{\mu}$ represents a covariant derivative and $T^{\mu\nu}$ is the
energy-momentum tensor. This yields 4 (3 space, 1 time) equations of motion
which dictate the evolution of the system. For most phenomenological
applications, the initial states of nuclear collisions and their subsequent
evolution are sufficiently complicated that they can only be solved
numerically. However, in certain highly symmetric scenarios, it may be
possible to find exact solutions with at least some phenomenological relevance
to nuclear collisions. The challenge is then to identify a suitable symmetry
group which mimics the dynamics of nuclear collisions sufficiently closely.

The original solution obtained by Gubser was constructed in order to respect
the symmetry group $SO(3)_{q}\times SO(1,1)\times\mathbb{Z}_{2}$, where
$SO(1,1)$ represents boost-invariance along the beam axis, $\mathbb{Z}_{2}$
represents a reflection symmetry under $\eta\to-\eta$, and $SO(3)_{q}$ denotes
a conformal symmetry which includes azimuthal rotations around the beam axis
as a subgroup and is designed in such a way as to generate a non-trivial flow
profile in the radial direction. Since the generators of the $SO(3)_{q}$ and
$SO(1,1)$ commute, one thus obtains a solution to the hydrodynamic equations
(1) which exhibits radial flow while also respecting boost invariance.

The $dS_{3}\otimes\mathbb{R}$ frame is parameterized by the target space
coordinates $\hat{x}^{\mu}=(\rho,\theta,\phi,\eta)$, while the physical Milne
coordinates are $x^{\mu}=(\tau,r,\phi,\eta)$.111Here and below, hatted
quantities (e.g., $\hat{X}$) are defined in the $dS_{3}\otimes\mathbb{R}$
frame, while quantities with out hats (e.g., $X$) are defined in physical
Milne coordinates. In both sets of coordinates, $\phi$ represents the
azimuthal angle around the beam axis and $\eta$ represents the longitudinal
rapidity, while in Milne coordinates, $\tau$ and $r$ represent the Bjorken
proper time and radial coordinate, respectively. Then the
$dS_{3}\otimes\mathbb{R}$ coordinates $\rho$ and $\theta$ are related to Milne
coordinates by

$\displaystyle\sinh\rho$
$\displaystyle=\frac{q^{2}\left(\tau^{2}-r^{2}\right)-1}{2q\tau}$ (2)
$\displaystyle\tan\theta$
$\displaystyle=\frac{2qr}{1+q^{2}\left(\tau^{2}-r^{2}\right)}$ (3)

and thus offer an alternative parametrization of the $(\tau,r)$ directions. In
terms of the $dS_{3}\otimes\mathbb{R}$ coordinates, boost invariance can be
expressed as independence of the rapidity $\eta$, while the $SO(3)_{q}$
symmetry amounts to independence of the coordinates $(\theta,\phi)$. The
$dS_{3}\otimes\mathbb{R}$ frame has the metric

$\displaystyle\hat{g}_{\mu\nu}=\mathrm{diag}\left(1,-\cosh^{2}\rho,-\cosh^{2}\rho\sin^{2}\theta,-1\right)$
(4)

while in Milne coordinates the metric is

$\displaystyle g_{\mu\nu}=\mathrm{diag}\left(1,-1,-r^{2},-\tau^{2}\right)$ (5)

To obtain Gubser’s original solution, one assumes the fluid is stationary in
the $dS_{3}\otimes\mathbb{R}$ frame:

$\displaystyle\hat{u}_{\mu}=(1,0,0,0)$ (6)

By transforming the trivial solution (6) back to Milne coordinates, Gubser
obtained the solution Gubser (2010); Gubser and Yarom (2011)

$\displaystyle u^{\tau}$ $\displaystyle=\tau\frac{\partial\rho}{\partial\tau}$
$\displaystyle=\frac{1+q^{2}\left(\tau^{2}+r^{2}\right)}{\sqrt{1+2q^{2}\left(\tau^{2}+r^{2}\right)+q^{4}\left(\tau^{2}-r^{2}\right)^{2}}}$
(7) $\displaystyle u^{r}$ $\displaystyle=-\tau\frac{\partial\rho}{\partial r}$
$\displaystyle=\frac{2q^{2}\
r\tau}{\sqrt{1+2q^{2}\left(\tau^{2}+r^{2}\right)+q^{4}\left(\tau^{2}-r^{2}\right)^{2}}}$
(8)

with a corresponding energy density $\epsilon$ given by

$\displaystyle\epsilon$
$\displaystyle=\frac{\hat{\epsilon}_{0}}{\tau^{4/3}}\frac{(2q)^{8/3}}{\left[1+2q^{2}(\tau^{2}+r^{2})+q^{4}(\tau^{2}-r^{2})^{2}\right]^{4/3}}$
(9)

Eqs. (7)-(9) define Gubser’s ideal solution to conformal hydrodynamics. It
contains two free parameters: $q$, which represents an inverse lengthscale
characterizing of the flow profile, and $\hat{\epsilon}_{0}$, which sets the
scale for the energy density.

###  II.2 Derivation of Solution I

We now want to generalize Gubser’s solution by relaxing the assumption of
boost invariance. In order to break the boost invariance of the Gubser
solution, we may introduce $\rho$ and $\eta$ dependence without violating the
$SO(3)_{q}$ symmetry. Explicitly, we take the following ansatz for the flow
(cf. Gubser (2010); Gubser and Yarom (2011)):

$\displaystyle\hat{u}_{\mu}=(\cosh\xi(\rho,\eta),0,0,\sinh\xi(\rho,\eta))$
(10)

This ansatz reduces to a stationary fluid (and therefore to Gubser’s solution)
whenever $\xi=0$. We emphasize that boost invariance is broken by the
inclusion of $\eta$-dependence, but the $SO(3)_{q}$ symmetry is unaffected.
Note that (10) is properly normalized using the $dS_{3}\otimes\mathbb{R}$
metric:

$\displaystyle\hat{u}^{\mu}\hat{u}^{\nu}\hat{g}_{\mu\nu}=1$ (11)

The Gubser solution requires a conformal equation of state. For simplicity, we
assume in this study that $\hat{p}=\hat{\epsilon}/3$. In this case, the ideal
energy-momentum tensor
$\hat{T}^{\mu\nu}=\hat{\epsilon}\hat{u}^{\mu}\hat{u}^{\nu}+\hat{p}\left(\hat{u}^{\mu}\hat{u}^{\nu}-\hat{g}^{\mu\nu}\right)$
has the following non-vanishing components:

$\displaystyle\hat{T}^{\rho\rho}$
$\displaystyle=\frac{\hat{\epsilon}}{3}\left(4\left(\hat{u}^{\rho}\right)^{2}-1\right)$
(12) $\displaystyle\hat{T}^{\rho\eta}$
$\displaystyle=\frac{4\hat{\epsilon}}{3}\hat{u}^{\rho}\hat{u}^{\eta}$ (13)
$\displaystyle=\hat{T}^{\eta\rho}$ (14) $\displaystyle\hat{T}^{\theta\theta}$
$\displaystyle=\frac{\hat{\epsilon}}{3}\left(4\left(\hat{u}^{\theta}\right)^{2}+\operatorname{sech}^{2}\rho\right)=\frac{\hat{\epsilon}}{3}\operatorname{sech}^{2}\rho$
(15) $\displaystyle\hat{T}^{\phi\phi}$
$\displaystyle=\frac{\hat{\epsilon}}{3}\left(4\left(\hat{u}^{\phi}\right)^{2}+\operatorname{sech}^{2}\rho\csc^{2}\theta\right)=\frac{\hat{\epsilon}}{3}\operatorname{sech}^{2}\rho\csc^{2}\theta$
(16) $\displaystyle\hat{T}^{\eta\eta}$
$\displaystyle=\frac{\hat{\epsilon}}{3}\left(4\left(\hat{u}^{\eta}\right)^{2}+1\right)$
(17)

This energy-momentum tensor must be covariantly conserved:

$\displaystyle\nabla_{\mu}\hat{T}^{\mu\nu}=\partial_{\mu}\hat{T}^{\mu\nu}+\Gamma^{\mu}_{\alpha\mu}\hat{T}^{\alpha\nu}+\Gamma^{\nu}_{\alpha\mu}\hat{T}^{\alpha\mu}=0$
(18)

Evaluating the Christoffel symbols for the metric (4) yields the following
non-vanishing components:

$\displaystyle\Gamma^{\rho}_{\theta\theta}$
$\displaystyle=\cosh\rho\sinh\rho,\qquad\Gamma^{\rho}_{\phi\phi}=\sin^{2}\theta\cosh\rho\sinh\rho$
(19) $\displaystyle\Gamma^{\theta}_{\rho\theta}$
$\displaystyle=\Gamma^{\phi}_{\rho\phi}=\tanh\rho,\qquad\Gamma^{\theta}_{\phi\phi}=-\sin\theta\cos\theta,\qquad\Gamma^{\phi}_{\theta\phi}=\cot\theta$
(20)

Then (1) implies two independent equations in the $\rho$ and $\eta$
directions:

$\displaystyle 0$
$\displaystyle=\partial_{\mu}\hat{T}^{\mu\rho}+\Gamma^{\mu}_{\alpha\mu}\hat{T}^{\alpha\rho}+\Gamma^{\rho}_{\alpha\mu}\hat{T}^{\alpha\mu}$
$\displaystyle=\partial_{\rho}\hat{T}^{\rho\rho}+\partial_{\eta}\hat{T}^{\rho\eta}+2\tanh\rho\,\hat{T}^{\rho\rho}+\cosh\rho\sinh\rho\left(\hat{T}^{\theta\theta}+\sin^{2}\theta\hat{T}^{\phi\phi}\right)$
(21) $\displaystyle 0$
$\displaystyle=\partial_{\mu}\hat{T}^{\mu\eta}+\Gamma^{\mu}_{\alpha\mu}\hat{T}^{\alpha\eta}+\Gamma^{\eta}_{\alpha\mu}\hat{T}^{\alpha\mu}$
$\displaystyle=\partial_{\rho}\hat{T}^{\rho\eta}+\partial_{\eta}\hat{T}^{\eta\eta}+2\tanh\rho\,\hat{T}^{\rho\eta}$
(22)

Defining $\hat{\epsilon}=\exp(4\mathcal{T})$ and making use of Eq. (10), the
conservation equations (21)-(22) simplify to

$\displaystyle-\tanh\rho$
$\displaystyle=\frac{\partial\mathcal{T}}{\partial\rho}+\cosh(2\xi)\left(\tanh\rho+\frac{\partial\xi}{\partial\eta}+2\frac{\partial\mathcal{T}}{\partial\rho}\right)+\sinh(2\xi)\left(\frac{\partial\xi}{\partial\rho}+2\frac{\partial\mathcal{T}}{\partial\eta}\right)$
(23) $\displaystyle 0$
$\displaystyle=-\frac{\partial\mathcal{T}}{\partial\eta}+\sinh(2\xi)\left(\tanh\rho+\frac{\partial\xi}{\partial\eta}+2\frac{\partial\mathcal{T}}{\partial\rho}\right)+\cosh(2\xi)\left(\frac{\partial\xi}{\partial\rho}+2\frac{\partial\mathcal{T}}{\partial\eta}\right)$
(24)

So far no assumptions have been made. In order to find an exact solution to
these equations, we now consider the special case where both the flow rapidity
$\xi$ and $\mathcal{T}$ are independent of $\eta$.222Note that $\eta$
dependence may still be included below by applying a temporal shift, as we
describe in Sec. II.4. In this case, they can depend only on $\rho$, so that
the equations simplify further to

$\displaystyle 0$
$\displaystyle=2\cosh^{2}\xi\tanh\rho+\left(1+2\cosh\left(2\xi\right)\right)\frac{d\mathcal{T}}{d\rho}+\sinh\left(2\xi\right)\frac{d\xi}{d\rho}$
(25) $\displaystyle 0$
$\displaystyle=\sinh\left(2\xi\right)\left(\tanh\rho+2\frac{d\mathcal{T}}{d\rho}\right)+\cosh\left(2\xi\right)\frac{d\xi}{d\rho}$
(26)

The second equation can be solved to find $\mathcal{T}$ in terms of $\rho$ and
$\xi(\rho)$; the solution is

$\displaystyle\mathcal{T}(\rho)$
$\displaystyle=c_{1}-\frac{1}{2}\ln\cosh\rho-\frac{1}{4}\ln\sinh\left(2\xi(\rho)\right),$
(27)

with $c_{1}$ an arbitrary constant. Substituting this solution for
$\mathcal{T}$ back into the first equation yields

$\displaystyle\sinh(2\xi)\tanh\rho-\left(2+\cosh(2\xi)\right)\frac{d\xi}{d\rho}=0$
(28)

This equation is trivially satisfied when $\xi=0$. In this case, one would
have found the solution for $\mathcal{T}$ to be

$\displaystyle\mathcal{T}(\rho)$
$\displaystyle=\tilde{c}_{1}-\frac{2}{3}\ln\cosh\rho,$ (29)

for some constant $\tilde{c}_{1}$, implying that

$\displaystyle\hat{\epsilon}\propto\left(\cosh\rho\right)^{-8/3},$ (30)

which is simply Gubser’s original solution Gubser and Yarom (2011). The ansatz
(10) thus provides a generalization of Gubser flow which reduces to Gubser’s
ideal solution when $\xi=0$, as anticipated.

To solve Eq. (28) when $\xi\neq 0$, we recast it into a more manageable form
using the changes of variables $\mathcal{x}=\cosh^{2}\rho$ and
$\xi=\frac{1}{2}\ln{\mathcal{y}(\mathcal{x})}$:

$\displaystyle\frac{d\mathcal{y}}{d\mathcal{x}}=\frac{\mathcal{y}\left(\mathcal{y}^{2}-1\right)}{\mathcal{x}\left(1+\mathcal{y}\left(\mathcal{y}+4\right)\right)}$
(31)

It has the solution

$\displaystyle\mathcal{y}(\mathcal{x})\equiv{\tilde{\mathcal{y}}\left(c_{2}\mathcal{x}\right)},$
(32)

where

$\displaystyle{\tilde{\mathcal{y}}\left(w\right)}$
$\displaystyle=1-\frac{w}{3}+\frac{2^{4/3}w\left(w-9\right)}{a(w)}+2^{2/3}a(w),$
(33) $\displaystyle a(w)$
$\displaystyle=\left[-w\left(54+w(2w-27)+3\sqrt{3\left(108-w^{2}\right)}\right)\right]^{1/3},$
(34)

and $c_{2}$ is another arbitrary constant. Our solution in the
$dS_{3}\otimes\mathbb{R}$ frame can therefore be written

$\displaystyle\mathcal{r}$
$\displaystyle=\cosh^{2}\rho=\frac{1+2q^{2}\left(\tau^{2}+r^{2}\right)+q^{4}\left(\tau^{2}-r^{2}\right)^{2}}{4q^{2}\tau^{2}}$
(35) $\displaystyle\hat{\epsilon}(\rho)$
$\displaystyle=\frac{2e^{4c_{1}}}{\mathcal{r}}\left(\frac{\tilde{\mathcal{y}}(c_{2}\mathcal{r})}{\tilde{\mathcal{y}}(c_{2}\mathcal{r})^{2}-1}\right)$
(36) $\displaystyle\xi(\rho)$
$\displaystyle=\pm\frac{1}{2}\ln{\tilde{\mathcal{y}}(c_{2}\mathcal{r})}$ (37)

and $c_{1}$, $c_{2}<0$ are undetermined constants. Note that the energy
density $\hat{\epsilon}$ is conformal; the physical energy density is
$\epsilon=\hat{\epsilon}/\tau^{4}$ in Milne coordinates. Note also that the
rapidity $\xi$ can be chosen positive or negative by the reflection symmetry
which is apparent in the above equations. Here we choose only the positive
‘branch’ of $\xi$ for the sake of definiteness.

The flow $u^{\mu}$ (10) can then be mapped to Milne coordinates by the
transformation rules Gubser and Yarom (2011)

$\displaystyle u^{\tau}\left(\tau,r\right)$
$\displaystyle=\tau\frac{\partial\rho}{\partial\tau}\cosh\xi$ (38)
$\displaystyle u^{r}\left(\tau,r\right)$
$\displaystyle=-\tau\frac{\partial\rho}{\partial r}\cosh\xi$ (39)
$\displaystyle u^{\eta}\left(\tau,r\right)$
$\displaystyle=-\frac{1}{\tau}\sinh\xi$ (40)

We thus arrive at an exact solution for the energy density $\epsilon$ and the
flow profile $u^{\mu}$, given in Eqs. (33)-(40). Boost invariance has clearly
not yet been broken, since the solution is still independent of $\eta$.

###  II.3 Asymptotics

Our exact solution (33)-(40) is singular in the limit $c_{2}\to 0^{-}$ unless
$c_{1}$ is chosen in an appropriate way (cf. (33) and (36)). To avoid this, we
define

$\displaystyle
c_{1}=\frac{1}{12}\ln\left(-2c_{2}\hat{\epsilon}_{0}^{3}\right)$ (41)

where the parameter $\hat{\epsilon}_{0}$ is the same one used in Gubser’s
solution Gubser (2010). This combination ensures that the energy density
$\hat{\epsilon}$ is regular and reduces to Gubser’s original solution as
$c_{2}\to 0^{-}$. One can see this most easily by expanding the solution for
$\hat{\epsilon}(\rho)$ in a power series for small $c_{2}<0$:

$\displaystyle\frac{\hat{\epsilon}(\rho)}{\hat{\epsilon}_{0}}=\frac{1}{\left(\cosh\rho\right)^{8/3}}-\frac{\alpha}{\left(\cosh\rho\right)^{4/3}}+\frac{3\alpha^{2}}{4}+O(\alpha^{3}\left(\cosh\rho\right)^{4/3})$
(42)

where we have also defined $\alpha\equiv(-c_{2})^{2/3}/(3\cdot 2^{1/3})$. The
general term in this expansion is clearly of order
$\alpha^{n}\left(\cosh\rho\right)^{4(n-2)/3}$, so that the series reduces to
(30) when $\alpha^{3/4}\ll\cosh\rho$, i.e., as $c_{2}\to 0^{-}$.

On the other hand, for sufficiently large $\rho$ one finds that the
$\tilde{\mathcal{y}}$ decreases inversely with its argument:
$\tilde{\mathcal{y}}(w)\propto 1/w$. Plugging this behavior into (36), we find
that the conformal energy density scales for large $\rho$ as

$\displaystyle\hat{\epsilon}(\rho)\sim\mathcal{r}^{-2}\sim e^{-4\rho}$ (43)

One therefore expects a regime in $\rho$ where the leading behavior of
$\hat{\epsilon}(\rho)$ transitions from that of Gubser’s solution (30) to a
steeper exponential fall-off in $\rho$ (43). The transition should happen in
the neighborhood of some critical value $\rho=\rho_{*}$ where the sub-leading
corrections in (42) become comparable to the leading term:

$\displaystyle\frac{1}{\left(\cosh\rho_{*}\right)^{8/3}}\approx\frac{\alpha}{\left(\cosh\rho_{*}\right)^{4/3}}\longrightarrow\rho_{*}\approx\operatorname{arccosh}\left(\alpha^{-3/4}\right)=\operatorname{arccosh}\left(\frac{54^{1/4}}{\sqrt{-c_{2}}}\right)$
(44)

Interestingly, we note that the definition (44) implies a constant value for
$\xi\left(\rho_{*}\right)\equiv\xi_{*}\approx 1.188$, since the $c_{2}$
dependence cancels in Eq. (37). Thus, while the precise “Gubser time"
$\rho_{*}$ at which the sub-leading corrections become important depends on
the value of $c_{2}$, it occurs at the same flow rapidity $\xi_{*}$,
regardless of $c_{2}$.

Figure 1: The conformal energy density $\hat{\epsilon}$ (panel (a)) and the
flow rapidity $\xi$ (panel (b)) vs. the time-like coordinate $\rho$, for
different choices of the free parameter $c_{2}$. Gubser’s original solution
(corresponding to $c_{2}=0$ is given by the thick (solid, black) line in both
panels. The horizontal dashed black line in the righthand panel shows that
$\xi_{*}=\xi\left(\rho_{*}\right)\approx 1.188$, independent of $c_{2}$. See
text for discussion.

This behavior is illustrated in Fig. 1, which shows the conformal energy
density $\hat{\epsilon}$ (panel (a)) and the flow rapidity $\xi$ (panel (b))
for several choices of $c_{2}$. For $\hat{\epsilon}$, one indeed observes that
for sufficiently small $\rho$, our solution is close to Gubser’s original
solution (solid black line), while for sufficiently large $\rho$, the solution
instead scales like $e^{-4\rho}$ (black dotted line). Moreover, the point
$\rho_{*}$ where the transition occurs depends on $c_{2}$ roughly as implied
by (44). For example, for $c_{2}=-10^{-3}$, we find using Eq. (44) that the
critical value $\rho_{*}\approx 5.14$, in good agreement with Fig. 1. Several
other values of $\rho_{*}$ are marked for the various $c_{2}$ shown.

Fig. 1(b) shows that $\xi$ is a monotonically increasing function for positive
$\rho$, and by symmetry under $\eta\to-\eta$ is monotonically decreasing for
negative $\rho$ (cf. (35)). As $c_{2}\to 0^{-}$, however, the solution again
converges smoothly back to Gubser’s boost invariant solution. Note that, since
$\xi\neq 0$ whenever $c_{2}\neq 0$, $u^{\eta}$ is always non-zero except when
$c_{2}=0$.

Finally, we emphasize that, aside from the reflection symmetry $\eta\to-\eta$,
our unshifted solution preserves all of Gubser’s original symmetries
(including rotational invariance around the beam axis and boost invariance).
In the next subsection, we will show how to modify this solution to break
boost invariance as well.

###  II.4 The shifted solution

The new solution found above consists of a non-trivial solution for $\epsilon$
and $u^{\mu}$ in Milne coordinates, which introduces a non-vanishing
$u^{\eta}$ while remaining independent of $\eta$. In this sense, the solution
is analogous to an azimuthally invariant system which nevertheless rotates in
the azimuthal direction. However, our original goal was to construct a
solution without boost invariance. As pointed out in Ref. Shi  _et al._
(2022), additional solutions possessing non-trivial $\eta$ dependence can be
generated by simply translating a known solution by a constant amount in the
Minkowski temporal coordinate: $t\to t^{\prime}=t+t_{0}$. We call this
transformation a “temporal shift." This is equivalent to a change of Milne
coordinates given by

$\displaystyle\tau\to\tau^{\prime}\\!\left(\tau,\eta,r\right)$
$\displaystyle=\sqrt{\tau^{2}+2t_{0}\tau\cosh\eta+t_{0}^{2}}$ (45)
$\displaystyle\eta\to\eta^{\prime}\\!\left(\tau,\eta,r\right)$
$\displaystyle=\operatorname{arctanh}\left(\frac{\tau\sinh\eta}{\tau\cosh\eta+t_{0}}\right)$
(46) $\displaystyle r\to r^{\prime}\\!\left(\tau,\eta,r\right)$
$\displaystyle=r$ (47)

Notice that the transformation on the $r$ coordinate is trivial, since the
temporal shift affects only $\tau$ and $\eta$, but we include a shifted
$r^{\prime}$ as well for clarity below.

By definition, the energy density $\epsilon$ transforms as a scalar:333Here
and below, quantities marked with primes (′) correspond to the unshifted
frame; unprimed quantities are defined in the shifted frame.

$\displaystyle\epsilon\left(\tau,r,\eta\right)=\epsilon^{\prime}\\!\left(\tau^{\prime},r\right)$
(48)

where here and below, $\tau^{\prime}$ is defined as in Eq. (45). Since the
flow velocity is a four-vector, it transforms as

$\displaystyle u^{\tau^{\prime}}$
$\displaystyle=\frac{\partial\tau^{\prime}}{\partial\tau}u^{\tau}+\frac{\partial\tau^{\prime}}{\partial\eta}u^{\eta}+\frac{\partial\tau^{\prime}}{\partial
r}u^{r}$ (49) $\displaystyle u^{\eta^{\prime}}$
$\displaystyle=\frac{\partial\eta^{\prime}}{\partial\tau}u^{\tau}+\frac{\partial\eta^{\prime}}{\partial\eta}u^{\eta}+\frac{\partial\eta^{\prime}}{\partial
r}u^{r}$ (50) $\displaystyle u^{r^{\prime}}$ $\displaystyle=\frac{\partial
r^{\prime}}{\partial\tau}u^{\tau}+\frac{\partial
r^{\prime}}{\partial\eta}u^{\eta}+\frac{\partial r^{\prime}}{\partial r}u^{r}$
(51)

Thanks to the triviality of the transformation in the $r$ direction, the
radial component of the flow velocity transforms only by shifting the
coordinates in Eq. (39) appropriately:

$\displaystyle u^{r}(\tau,r,\eta)$
$\displaystyle=u^{r^{\prime}}\\!\left(\tau^{\prime},r\right)$ (52)

The remaining components of the flow velocity are then found by solving (49)
and (50) for $u^{\tau}$ and $u^{\eta}$. The result is (with $c_{1}$ set as in
Eq. (41)):

$\displaystyle\mathcal{r}^{\prime}$
$\displaystyle=\cosh^{2}\rho^{\prime}=\frac{1+2q^{2}\left({\tau^{\prime}}^{2}+r^{2}\right)+q^{4}\left({\tau^{\prime}}^{2}-r^{2}\right)^{2}}{4q^{2}{\tau^{\prime}}^{2}}$
(53) $\displaystyle\frac{\hat{\epsilon}(\rho^{\prime})}{\hat{\epsilon}_{0}}$
$\displaystyle=\frac{\left(-16c_{2}\right)^{1/3}}{\mathcal{r}^{\prime}}\left(\frac{\tilde{\mathcal{y}}(c_{2}\mathcal{r}^{\prime})}{\tilde{\mathcal{y}}(c_{2}\mathcal{r}^{\prime})^{2}-1}\right)$
(54) $\displaystyle\xi(\rho^{\prime})$
$\displaystyle=\frac{1}{2}\ln{\tilde{\mathcal{y}}(c_{2}\mathcal{r}^{\prime})}$
(55) $\displaystyle u^{\tau}(\tau,\eta,r)$
$\displaystyle=\left(\frac{\tau+t_{0}\cosh\eta}{{\tau^{\prime}}}\right)u^{\tau^{\prime}}\\!\left({\tau^{\prime}},r\right)-\left(t_{0}\sinh\eta\right)u^{\eta^{\prime}}\\!\left({\tau^{\prime}},r\right)$
(56) $\displaystyle u^{r}(\tau,\eta,r)$
$\displaystyle=u^{r^{\prime}}({\tau^{\prime}},r)$ (57) $\displaystyle
u^{\eta}(\tau,\eta,r)$
$\displaystyle=\left(1+\frac{t_{0}}{\tau}\cosh\eta\right)u^{\eta^{\prime}}\\!\left({\tau^{\prime}},r\right)-\frac{t_{0}}{\tau}\left(\frac{\sinh\eta}{{\tau^{\prime}}}\right)u^{\tau^{\prime}}\\!\left({\tau^{\prime}},r\right)$
(58)

Eqs. (53)-(58), together with (45), constitute our complete new Solution I
which generalizes Gubser flow to include non-boost-invariant longitudinal
flow, where the unshifted flow components $u^{\tau^{\prime}}$,
$u^{r^{\prime}}$, and $u^{\eta^{\prime}}$ have been defined in (38)-(40).

###  II.5 Properties of Solution I

We now consider how our Solution I obtained above depends on the parameters
$c_{2}$ and $t_{0}$. For definiteness, all other parameters are fixed to their
values as defined in Refs. Gubser (2010); Gubser and Yarom (2011).
Specifically, we take

$q=\left(4.3\text{ fm}\right)^{-1}\text{ and }\hat{\epsilon}_{0}=880$

which are chosen in order to have a reasonable description of central Au+Au
collisions with $\sqrt{s_{NN}}=200A$ GeV.

Figure 2: The physical energy density $\epsilon$ as a function of $\tau$, for
various radius $r$ and rapidity $\eta$. Each panel shows a different
combination of free parameters $t_{0}$ and $c_{2}$: $t_{0}=0.1$ fm/$c$ and
$c_{2}=0$ (panel (a)); $t_{0}=0.1$ fm/$c$ and $c_{2}=-0.1$ (panel (b));
$t_{0}=0.5$ fm/$c$ and $c_{2}=0$ (panel (c)); $t_{0}=0.5$ fm/$c$ and
$c_{2}=-0.1$ (panel (d)).

We begin by considering how the physical energy density depends on $c_{2}$ and
$t_{0}$. This is shown in Figure 2. We observe several qualitative behaviors.
First, by comparing panels (a) and (b) with the panels (c) and (d), we see
that increasing $t_{0}$ has the effect of steepening the rapidity dependence
of the shifted solution. This is intuitively plausible, since setting
$t_{0}=0$ must lead back to an $\eta$-independent solution. A similar effect
is observed by making $c_{2}$ non-zero: the rapidity dependence is weakened as
the value of $c_{2}$ is decreased. Notice further that the effects of
modifying $t_{0}$ or $c_{2}$ affect $\epsilon$ in different ways: increasing
$t_{0}$ causes $\epsilon$ to fall faster at late $\tau$, while decreasing
$c_{2}$ leads to a faster decrease of $\epsilon$ in the radial direction and
only weakly affects the $\tau$ dependence. We only show curves with $\eta\geq
0$ since the energy density is always symmetric under $\eta\to-\eta$ (cf. (54)
and (45)).

Figure 3: The transverse velocity $u^{r}/u^{\tau}$ as a function of $\eta$,
for various $\tau$ and $r$. Each panel shows a different combination of free
parameters $t_{0}$ and $c_{2}$: $t_{0}=0.1$ fm/$c$ and $c_{2}=0$ (panel (a));
$t_{0}=0.1$ fm/$c$ and $c_{2}=-0.1$ (panel (b)); $t_{0}=0.5$ fm/$c$ and
$c_{2}=0$ (panel (c)); $t_{0}=0.5$ fm/$c$ and $c_{2}=-0.1$ (panel (d)). See
text for discussion.

Our solution’s flow profile is shown in Figs. 3 and 4. In Fig. 3, we show the
ratio $u^{r}/u^{\tau}$ plotted against rapidity $\eta$, for different
combinations of $r$, $\tau$, $t_{0}$, and $c_{2}$. We observe that increasing
$t_{0}$ (cf. panels (a) and (b) with panels (c) and (d)) again reduces the
solution’s width in $\eta$. Decreasing $c_{2}$, however, manifestly breaks the
reflection symmetry $\eta\to-\eta$ as anticipated above (cf. panels (a) and
(c) with panels (b) and (d)): having chosen the ‘positive’ branch of $\xi$ for
our solution in Eq. (37), we see that the resulting transverse flow is
stronger for $\eta<0$. This is a consequence of the fact that a positive flow
rapidity $\xi$ yields a negative rapidity component $u^{\eta}$, by Eq. (40).
Choosing $c_{2}<0$ therefore leads to an asymmetric flow in $\eta$, by (56)
and (58).

Figure 4: The longitudinal velocity $\tau u^{\eta}/u^{\tau}$ as a function of
$\eta$, for various $\tau$ and $r$. Each panel shows a different combination
of free parameters $t_{0}$ and $c_{2}$: $t_{0}=0.1$ fm/$c$ and $c_{2}=0$
(panel (a)); $t_{0}=0.1$ fm/$c$ and $c_{2}=-0.1$ (panel (b)); $t_{0}=0.5$
fm/$c$ and $c_{2}=0$ (panel (c)); $t_{0}=0.5$ fm/$c$ and $c_{2}=-0.1$ (panel
(d)). See text for discussion.

In Fig. 4, we show the corresponding results for the longitudinal flow
velocity, $\tau u^{\eta}/u^{\tau}$. In panels (a) and (c), we see that the
only $\eta$ dependence arises from the shift to $t_{0}\neq 0$; since
$c_{2}=0$, there is no additional rapidity dependence which arises from
breaking the reflection symmetry. In this case $\tau u^{\eta}/u^{\tau}$
depends only on $\tau$ and $\eta$, but not on $r$ (cf. Eqs. (56) and (58)).
One also observes that increasing $t_{0}$ steepens the rapidity dependence, as
was noted above; setting $t_{0}=0$ would result in $u^{\eta}=0$ for all $\tau$
and $r$. A more nuanced picture of the flow emerges once we let $c_{2}<0$
(shown in Fig. 4(b) and (d)). In this case, the location with the longitudinal
flow vanishes shifts to $\eta<0$, where the size of the shift depending on
both $\tau$ and $r$. One therefore finds that the regions where the rapidity
flow ($u^{\eta}$) is weakest are also those where the radial flow ($u^{r}$) is
strongest. This is intuitively what one would expect for an incompressible
fluid. Note, however, that this rough correspondence is not exact: even in the
symmetric case ($c_{2}=0$), $u^{r}/u^{\tau}$ does not reach a maximum at
$\eta=0$ where the longitudinal flow vanishes. This is due to the complicated
nature of Gubser flow and the effects of translating in the temporal direction
(cf. panels (a) and (b) in Figs. 3 and 4).

Our solution therefore exhibits some curious properties: the energy density is
always symmetric under $\eta\to-\eta$, while the flow profile is asymmetric
whenever $c_{2}<0$. This suggests that scenarios with $c_{2}<0$ may not bear
much phenomenological relevance to the modeling of nuclear collisions, where
one would expect to find both flow profiles and energy densities to be
asymmetric in $\eta$ Shi  _et al._ (2022). Nonetheless, this solution may also
be useful in generating extreme flow profiles which can be used to subject
hydrodynamic codes to valuable performance tests. In the specific case that
$c_{2}=0$, however, the temporal shift induces non-trivial $\eta$ dependence
without violating the $\mathbb{Z}_{2}$ symmetry in $\eta$, and thus amounts to
a shifted, non-boost-invariant generalization of Gubser’s original, boost-
invariant solution for ideal hydrodynamics. In the next section, we show how
to apply the same ‘temporal shift’ trick to the viscous version of Gubser’s
solution, for which the shear viscosity is non-zero.

##  III Solution II

The solution derived above in Sec. II was obtained assuming only ideal
hydrodynamics. However, Gubser Gubser (2010) originally derived a version of
his boost-invariant solution to the Navier-Stokes equations with a finite
specific shear viscosity $\tilde{\eta}/s\neq 0$.444Here and below, we denote
the shear viscosity by $\tilde{\eta}$, in order to distinguish it from the
space-time rapidity $\eta$. Gubser’s viscous solution is thus a generalization
of our Solution I (with $c_{2}=0$) to a finite shear viscosity. In this
section, we apply a temporal shift to introduce rapidity dependence also into
this viscous solution. This will allow us to define and examine the Knudsen
number in Sec. IV.

Gubser’s viscous solution is most conveniently written in terms of the
temperature $T_{\text{G}}$, where the subscript ‘G’ denotes ‘Gubser.’
$T_{\text{G}}$ is related to the energy density by $\epsilon\equiv
f_{*}T_{\text{G}}^{4}$ (where $f_{*}\equiv 11$) and is defined by Gubser
(2010); Gubser and Yarom (2011)

$\displaystyle T_{\text{G}}(\tau,r)$ $\displaystyle=\frac{\hbar
c}{f_{*}^{1/4}\tau}\hat{T}\\!\left(\rho(\tau,r)\right)$ (59)
$\displaystyle\hat{T}(\rho)$
$\displaystyle=\frac{\hat{T}_{0}}{\left(\cosh\rho\right)^{2/3}}\left(1+\frac{H_{0}}{9\hat{T}_{0}}\left(\sinh\rho\right)^{3}{{}_{2}F_{1}}\left(\frac{3}{2},\frac{7}{6},\frac{5}{2};-\sinh^{2}\rho\right)\right),$
(60)

where $\rho$ is defined by Eq. (2) and $H_{0}\propto\eta/s$. We set
$\hat{T}_{0}=5.55$ and $H_{0}=0.33$, corresponding to a choice of
$\eta/s=0.134$ Gubser (2010).

The flow $u_{\text{G}}^{\mu}$ in the viscous case is the same as that obtained
in the ideal case (7)-(9). One easily checks that
$\left(u_{\text{G}}^{\tau}\right)^{2}-\left(u_{\text{G}}^{r}\right)^{2}=1$.
Gubser’s complete viscous solution for $T_{\text{G}}$ and $u_{\text{G}}^{\mu}$
is thus given by Eqs. (7), (8), and (59).

To obtain the shifted version of Gubser’s viscous solution (denoted by the
subscript ‘sG’), we apply the same transformation discussed in Sec. II.4. Like
the energy density $\epsilon$, the temperature $T_{\text{G}}$ must transform
as a scalar. This implies that Gubser’s viscous solution with rapidity
dependence is given by

$\displaystyle T_{\text{sG}}(\tau,r,\eta)$
$\displaystyle=T_{\text{G}}(\tau^{\prime},r)$ (61)

In analogy with (57), the radial component of the shifted flow
$u_{\text{sG}}^{\mu}$ is

$\displaystyle u_{\text{sG}}^{r}(\tau,r,\eta)$
$\displaystyle=u_{\text{G}}^{r^{\prime}}\left(\tau^{\prime},r\right)$ (62)

The remaining components of the flow velocity are then found by solving
(49)-(50) for the shifted components as before, and recalling that
$u_{\text{G}}^{\eta^{\prime}}=0$ as a result of the boost invariance and
reflection symmetry of the Gubser solution. The result is

$\displaystyle u_{\text{sG}}^{\tau}(\tau,r,\eta)$
$\displaystyle=u_{\text{G}}^{\tau^{\prime}}\left(\tau^{\prime},r\right)\left(\frac{\tau+t_{0}\cosh\eta}{\sqrt{\tau^{2}+2t_{0}\tau\cosh\eta+t_{0}^{2}}}\right)$
(63) $\displaystyle u_{\text{sG}}^{\eta}(\tau,r,\eta)$
$\displaystyle=-u_{\text{G}}^{\tau^{\prime}}\left(\tau^{\prime},r\right)\left(\frac{t_{0}\sinh\eta}{\tau\sqrt{\tau^{2}+2t_{0}\tau\cosh\eta+t_{0}^{2}}}\right)$
(64)

Again, one easily confirms that the shifted flow still possesses the correct
normalization:

$\displaystyle\left(u_{\text{sG}}^{\tau}\right)^{2}-\left(u_{\text{sG}}^{r}\right)^{2}-\left(\tau
u_{\text{sG}}^{\eta}\right)^{2}=1$ (65)

Eqs. (61)-(64) thus provide a generalization of Gubser’s viscous solution
which relaxes the assumption of boost invariance by introducing a non-trivial
dependence on the space-time rapidity $\eta$. It coincides with the first
solution found above in Secs. II only if one takes $c_{2}=0$ and
$\tilde{\eta}/s=0$. In Sec. IV, we will consider several properties of this
generalized solution and their implications for the relationship between the
freeze out hypersurface and flow in relativistic nuclear collisions.

##  IV Flow and Freeze out

Having introduced the two new exact solutions derived above, we turn finally
to discuss how we can use them to better understand the role of flow in
influencing the freeze out process. Since we wish to compare freeze out at
constant temperature with that at constant Knudsen number, we will focus our
attention on Solution II.

In general, the criteria defining freeze out are normally formulated Jeon and
Heinz (2015) in terms of quantities which depend locally on the hydrodynamic
state of the system and which reflect (directly or indirectly) the extent to
which the system in a given region may be represented as a strongly coupled
plasma of deconfined color charges. The two criteria we consider here are:
(i), freeze out at constant temperature; and (ii), freeze out at constant
Knudsen number. Freeze out at constant temperature has historically been the
more widely used criterion, and may be interpreted either as a
phenomenological choice of where to ‘switch’ from a partonic to a hadronic
description of the system Ryu  _et al._ (2015); Weil  _et al._ (2016) or as a
way of characterizing the system’s passage through the
confinement/deconfinement transition Everett  _et al._ (2021). On the other
hand, as noted earlier, freeze out at constant Knudsen number characterizes
the partonic-hadronic transition in terms of the degree to which a fluid
dynamical description is applicable. Insofar as constant $T$ hypersurfaces
closely approximate those of constant $\mathrm{Kn}$ in realistic hydrodynamic
simulations Niemi and Denicol (2014), the distinction between them is largely
academic. Our goal here, however, is to explore in the context of an
analytical solution whether this identification always holds, and if not,
under what kinds of conditions it may fail.

The Knudsen number Kn is defined by Niemi and Denicol (2014); Jeon and Heinz
(2015)

$\displaystyle\mathrm{Kn}\equiv\frac{\ell_{\mathrm{mfp}}}{\ell_{\mathrm{macro}}}$
(66)

where $\ell_{\mathrm{mfp}}$ is the mean free path of the fluid and
$\ell_{\mathrm{macro}}$ is an appropriate macroscopic length scale compared to
which $\ell_{\mathrm{mfp}}$ should be sufficiently small in order to justify
the applicability of fluid dynamics.

In order to compute the Knudsen number, we must estimate both lengthscales
entering Eq. (66). In terms of the quantities employed here,
$\ell_{\mathrm{mfp}}$ can be estimated by Niemi and Denicol (2014); Jeon and
Heinz (2015)

$\ell_{\mathrm{mfp}}\sim\frac{5\tilde{\eta}}{sT}$

where $T$ is the temperature, and $\tilde{\eta}/s$ is the ratio of the shear
viscosity to entropy density. This choice of $\ell_{\mathrm{mfp}}$ corresponds
to the value of the shear relaxation time $\tau_{\pi}$ using the 14-moment
approximation of the relativistic Boltzmann equation in the massless limit
Denicol  _et al._ (2012); Molnár  _et al._ (2014). As noted above, since a
finite $\ell_{\mathrm{mfp}}\neq 0$ implies a nonvanishing shear viscosity
ratio $\tilde{\eta}/s\neq 0$, we limit our discussion in this section to the
shifted viscous Gubser solution obtained in Sec. III.

Many choices are possible for the macroscopic lengthscale
$\ell_{\mathrm{macro}}$ Niemi and Denicol (2014); here, we define it in terms
of the scalar expansion rate $\theta$,555Another possible choice is the
spatial gradient of the energy density, $\epsilon$, which is expected to be of
the same order of magnitude as the definition we use here Niemi and Denicol
(2014). which defines a corresponding lengthscale
$\ell_{\mathrm{macro}}=L_{\theta}$ given by

$\displaystyle\frac{1}{L_{\theta}}$
$\displaystyle\equiv\theta\left(\tau,r,\eta\right)=\partial_{\mu}u^{\mu}+\Gamma^{\mu}_{\mu\alpha}u^{\alpha}=\frac{1+q^{2}\left(r^{2}+5{\tau^{\prime}}^{2}\right)}{{\tau^{\prime}}\sqrt{1+2q^{2}\left({\tau^{\prime}}^{2}+r^{2}\right)+q^{4}\left({\tau^{\prime}}^{2}-r^{2}\right)^{2}}}$
(67)

where $\tau^{\prime}$ is defined in (45), and the non-vanishing Christoffel
symbols in Milne coordinates are

$\displaystyle\Gamma^{\tau}_{\eta\eta}=\tau,\qquad\Gamma^{\eta}_{\eta\tau}=\Gamma^{\eta}_{\tau\eta}=\frac{1}{\tau},\qquad\Gamma^{r}_{\phi\phi}=r,\qquad\Gamma^{\phi}_{\phi
r}=\Gamma^{\phi}_{r\phi}=\frac{1}{r}$ (68)

$L_{\theta}$ thus characterizes the typical lengthscale over which the
strength of the system flow changes appreciably. As pointed out above and in
Niemi and Denicol (2014), these scales should be much larger than the relevant
microscopic scale which determines the strength or frequency of interaction,
i.e., the mean free path $\ell_{\mathrm{mfp}}$, in order for fluid dynamics to
be applicable. The Knudsen number may therefore be defined as

$\displaystyle\mathrm{Kn}$
$\displaystyle\equiv\mathrm{Kn}_{\theta}\equiv\frac{\ell_{\mathrm{mfp}}}{L_{\theta}}=\frac{5\tilde{\eta}\theta}{sT}$
(69)

From (69), it follows that that the Knudsen number is small enough
($\mathrm{Kn}\ll 1$) to warrant the use of fluid dynamics only when either
$\tilde{\eta}/s$ (i.e., out of equilibrium viscous corrections) or $\theta$
(i.e., the flow gradients) is sufficiently small. In the former case, one may
have a violently expanding system which nevertheless is describable using
fluid dynamics, as long as the mean free path is sufficiently short and the
frequency of interactions sufficiently high. In the latter case, one may have
a relatively large mean free path but nevertheless maintain local equilibrium
if the system evolves sufficiently slowly. Conversely, when either
$\tilde{\eta}/s$ or $\theta$ becomes too large relative to the other, then the
interactions become too infrequent to maintain local equilibrium and thus
fluid dynamics no longer provides a justifiable description of the system. It
is at this point when freeze out in terms of $\mathrm{Kn}$ should be
implemented.

In the remainder of this section, we consider two questions, using our
Solution II as a concrete model:

  1. 1.

How much do the two freeze out criteria (constant temperature, constant
Knudsen number) differ in the freeze out contours which they produce?

  2. 2.

To what extent are these differences influenced by the strength of collective
flow and the size of its gradients throughout the system?

Figure 5: Pairs of freeze out contours using criteria of constant temperature
(solid black lines) and constant Knudsen number $\mathrm{Kn}$ (dashed black
lines), plotted on top of the transverse flow $v_{r}=u^{r}/u^{\tau}$ in the
$\tau$-$r$ plane. Both panels take $t_{0}=0.5$ fm$/c$ and either $\eta=0$
(panel (a), where $u^{\eta}=0$) or $\eta=2$ (panel (b)). The Knudsen number
contours correspond to $\mathrm{Kn}=$ 0.50, 0.75, 1.00, and 1.25 (ordered by
increasing $\tau$ at $r=0$), and the temperature contours have been adjusted
to match the Knudsen number contours at $r=0$. A solid red contour at $T=130$
MeV is included for reference.

To address the first question, in Fig. 5 we plot contours of fixed Knudsen
number (69) for Solution II (with $t_{0}=0.5\text{ fm}/c$) as a function of
$\tau$ and $r$ (displayed as solid, black lines), for $\eta=0$ (panel (a)) and
$\eta=2$ (panel (b)). We also show for reference an isothermal contour of
$T=130$ MeV (solid, red line) corresponding roughly to $\mathrm{Kn}\approx
0.75$, according to the definition $\eqref{Knudsen_number_max}$. The contours
are displayed on top of a density plot of the radial velocity
$v_{r}=u^{r}/u^{\tau}$ which indicates where the transverse flow is the
strongest. Four pairs of $T$/$\mathrm{Kn}$ (solid/dashed) contours are shown,
where each pair is adjusted to coincide when $r=0$. The Knudsen contours
assume values of $\mathrm{Kn}=$ 0.50, 0.75, 1.00, and 1.25 (corresponding at
$r=0$ respectively to $\tau\approx 2.2$, 4.7, 7.3, and 10.3 fm$/c$ and
temperatures of $T\approx 209$ MeV, 115 MeV, 69 MeV, and 44 MeV).

One finds that, for $\eta=0$, $\mathrm{Kn}\lesssim 0.5$ at sufficiently small
$r\lesssim 3$ fm and times which are neither too early nor too late
($\tau\lesssim 2$ fm$/c$), and thus that both sets of contours (constant $T$
and constant $\mathrm{Kn}$) demarcate similar regions in the $(\tau,r)$-plane
where the system may be treated fluid dynamically. Indeed, the contour of
$\mathrm{Kn}=0.5$ clearly agrees quantitatively with the corresponding contour
of constant $T$ at sufficiently small radii $r$ and late times $\tau$. Similar
statements apply to the $\eta=2$ case, except that all freeze out contours are
reached at earlier times and thus the regime where fluid dynamics is valid
shrinks accordingly.

For larger $\mathrm{Kn}$, Fig. 5 shows explicitly that the shapes of constant
temperature contours can differ quite significantly, both qualitatively and
quantitatively, from contours of constant Knudsen number. We find that the
contours of fixed $T$ acquire a positive concavity at $r=0$ at sufficiently
late times and low temperatures.666In the ideal case, which is a reasonable
approximation to the viscous solution at $r=0$, one easily shows that the
concavity of the temperature contour changes sign at a critical proper time of
$\tau_{\mathrm{c}}\equiv 1/(cq)=4.3$ fm$/c$, corresponding to an ideal
temperature of $T(\tau_{\mathrm{c}})=\hbar cq\hat{T}_{0}/f_{*}^{1/4}\approx
140$ MeV. This seems to agree quite well with where the temperature contours
in Fig. 5 become flat at $r=0$. This leads to the formation of a ‘wing-like’
structure at large radii which has been noted previously Heinz and Moreland
(2019) and which has been shown to reflect the violent collective expansion of
the system Plumberg (2020a, b). This is because the interior of the fireball
freezes out before the edges do, thanks to the rapid cooling produced by the
strong collective flow.

Moreover, one finds that, while the contours corresponding to
$\mathrm{Kn}=0.50$ and 0.75 agree qualitatively with one another, significant
discrepancies develop in the contour pairs with $\mathrm{Kn}=$ 1.00 and 1.25
in the region where transverse flow becomes the strongest. This has a
straightforward and intuitive interpretation: since $\mathrm{Kn}$ essentially
differs from $T$ only by a constant factor times $\theta$, the differences
between constant temperature and constant Knudsen contours must arise from the
expansion rate which pushes the constant $\mathrm{Kn}_{\theta}$ contours to
earlier times and smaller radii, compared with constant $T$ contours. Note
that since $u^{\eta}=0$ when $\eta=0$ (by Eq. (58)), in this case the
discrepancies between the contours are completely determined by
$u^{r}/u^{\tau}$, as well as the gradients of both transverse and longitudinal
flow.

In response to our first question above, then, we conclude that different
freeze out criteria such as decoupling at a fixed $T$ vs. a fixed
$\mathrm{Kn}$ may indeed lead to dramatically different contours and
consequently produce different geometries for the freeze out process. In the
present example, the largest differences occur at larger Knudsen numbers
$\mathrm{Kn}\gtrsim 1.00$ and hence smaller temperatures $T\lesssim 100$ MeV.
However, as we have already suggested and now demonstrate, the real origin of
these discrepancies is the effect of powerful collective flow with large
gradients which drives a rapid cooling of the system and forces an earlier
freeze out when it is sufficiently strong.

Figure 6: Panel (a) displays thick contours in the $(\tau,r)$-plane for
$\mathrm{Kn}_{\theta}=0.75$ (solid blue), 1.00 (dashed purple), and 1.25 (dot-
dashed red). For each of these respective contours, we also plot the
temperature $T$ (panel (b)) and radial velocity $u^{r}/u^{\tau}$ (panel (c))
as a function of radius $r$, with the same styling of lines in each panel. For
each thick contour representing $\mathrm{Kn}_{\theta}$, we also plot a thin
contour of the same color and linestyle to represent
$\widetilde{\mathrm{Kn}_{\theta}}$.

In order to facilitate a direct comparison between the contours of fixed $T$
and fixed $\mathrm{Kn}$ and to show how their differences are influenced by
collective flow, in Fig. 6 we plot the temperature $T$ and radial velocity
$v_{r}$ along the contours of constant $\mathrm{Kn}$ themselves, assuming
$t_{0}=0.5\text{ fm}/c$. In Fig. 6(a), we plot portions of the
$\mathrm{Kn}_{\theta}$ as thick $(\tau,r)$-contours at values of 0.75 (solid
blue), 1.00 (dashed purple), and 1.25 (dot-dashed red), extending at late
times from $r=0$ to approximately where $d\tau/dr$ along the contour diverges.
Unlike the constant $T$ contours in Fig. 5, the constant $\mathrm{Kn}$
contours exhibit negative concavity at $r=0$: freeze out occurs last at the
center of the fireball. In Fig. 6(b) and (c), we plot $T$ and $v_{r}$ along
these same contours as functions of $r$ and find that $T$ is decidedly not
constant along the $\mathrm{Kn}$ contours. Instead, thanks to the rapid growth
of transverse flow with $r$ and its competition with the temperature gradients
in the fireball’s interior Schnedermann  _et al._ (1993), one finds that $T$
first increases from $r=0$, before reaching a peak and then decreasing toward
larger radii. By inspecting the curves in Fig. 6(b), we observe that the
relative height of the peak above the central temperature is approximately
$13\%$ for $\mathrm{Kn}_{\theta}=0.75$, $37\%$ for
$\mathrm{Kn}_{\theta}=1.00$, and $72\%$ for $\mathrm{Kn}_{\theta}=1.25$.
Comparison with the curves in Fig. 6(c) then confirms that this relative
difference increases with the change in the maximum transverse velocity of the
fluid, indicating the strong influence of the latter on the former. We can
therefore answer our second question by confirming that strong collective flow
gradients can generate qualitatively significant discrepancies between
different freeze out criteria.777Note that only very slight discrepancies
between freeze out criteria were observed in Ref. Niemi and Denicol (2014) in
the case of p+Pb collisions where the same ‘wing-like’ structure was observed
in the constant $T$ contours, in apparent contrast to the results found here.
In this regard, however, we emphasize several key differences between our
exact solution and the calculations performed in Niemi and Denicol (2014):
(i), Gubser’s solution (and thus our Solution II) assumes Navier-Stokes
hydrodynamics, for which the shear correction $\pi^{\mu\nu}$ is directly
proportional to the shear tensor $\sigma^{\mu\nu}$, whereas Niemi and Denicol
(2014) employed Israel-Stewart hydrodynamics Israel and Stewart (1979), in
which $\pi^{\mu\nu}$ evolves according to a relaxation equation; (ii), Niemi
and Denicol (2014) did not use a constant $\tilde{\eta}/s$ as assumed here;
and (iii), it is unclear whether the Gubser solution’s flow gradients are
really comparable to those present in the p+Pb simulation shown in Niemi and
Denicol (2014). Some combination of these differences may account for the much
larger sensitivity to freeze out criterion seen here.

The tight connection between flow and the freeze out process may be further
clarified by considering how the Knudsen contours would be affected by the
omission of a portion of the flow. To see this, we construct an unphysical
flow which omits the $\eta$ component from $u^{\mu}$:

$\displaystyle\tilde{u}^{\mu}$
$\displaystyle=\left(\tilde{u}^{\tau},\tilde{u}^{r},0,0\right)\equiv\left(\sqrt{1+\left(u^{r}\right)^{2}},u^{r},0,0\right)$
(70)

with $u^{r}$ given by the shifted Gubser solution used in Eq. (67). We
emphasize that Eq. (70) is not a solution of the original equations of motion
(1), although it remains properly normalized by construction:

$\displaystyle\tilde{u}_{\mu}\tilde{u}^{\mu}=1$ (71)

Using $\tilde{u}^{\mu}$, we compute a corresponding unphysical Knudsen number
$\widetilde{\mathrm{Kn}}_{\theta}$ which may in turn be used to define
modified freeze-out contours. Specifically, we define

$\displaystyle\widetilde{\mathrm{Kn}}_{\theta}$
$\displaystyle=\frac{\tilde{\eta}}{sT}\left(\partial_{\mu}\tilde{u}^{\mu}+\Gamma^{\mu}_{\mu\alpha}\tilde{u}^{\alpha}\right)$
(72)

We use the same definition of $T$ for both $\mathrm{Kn}_{\theta}$ and
$\widetilde{\mathrm{Kn}}_{\theta}$ as in Eq. (67). The corresponding contours
are shown as thin lines in Fig. 6. We use three of the same values
($\mathrm{Kn}=0.75$, 1.00, 1.25) of the Knudsen numbers as before. By
comparing the respective thick and thin line pairs we may then estimate the
importance of the component $u^{\eta}$ for the shape of the freeze out
surface.

We first observe the effects of omitting $u^{\eta}$ on the shape of the
$\mathrm{Kn}_{\theta}$ contours themselves. By Eq. (58), we see that choosing
$t_{0}=0\text{ fm}/c$ would have resulted in $u^{\eta}=0$ and therefore no
differences between the $\mathrm{Kn}_{\theta}$ and
$\widetilde{\mathrm{Kn}}_{\theta}$ contours. With $t_{0}=0.5\text{ fm}/c$,
however, we see that the unphysical $\widetilde{\mathrm{Kn}}_{\theta}$
contours shift to smaller $\tau$ and $r$ values relative to the
$\mathrm{Kn}_{\theta}$ contours. Omitting $u^{\eta}$ from the shifted viscous
Gubser solution of Sec. III thus effectively enhances the flow gradients in
the system and shrinks the region where fluid dynamics is applicable. As a
result of the artificially shifted contours, the unphysical flow also steepens
the temperature and radial flow dependence along the contours. Looking at the
center panel of Fig. 6, we find that the discrepancies between the central
($r=0$) and peak temperatures along the contours are slightly enhanced along
the shifted $\widetilde{\mathrm{Kn}}_{\theta}$ contours. The effects on
$v_{r}/c$ are more pronounced, producing up to a 6% reduction in the peak
radial velocity for the $\mathrm{Kn}_{\theta}=0.75$ contour. Although we do
not show it here, the effects of omitting $u^{\eta}$ are enhanced further as
$t_{0}$ is increased. We conclude that one effect of setting $t_{0}>0$ is to
‘redirect’ some amount of the unshifted solution’s flow into the rapidity
direction, such that artificially suppressing the resulting $\eta$ component
of the flow generates spurious gradients which exacerbate the differences
between freeze out criteria. Both transverse and longitudinal flow therefore
play essential roles in determining the quantitatively correct shape of the
freeze out surface, as well as the variations in the temperature and radial
flow along it.

Figure 7: Panel (a) displays contours in the $(\tau,r)$-plane for
$\mathrm{Kn}_{\theta}=0.75$ (solid blue), 1.00 (dashed purple), and 1.25 (dot-
dashed red). For each of these respective contours, we also plot the
temperature $T$ (panel (b)) and radial velocity $u^{r}/u^{\tau}$ (panel (c))
as a function of radius $r$, with the same styling of lines in each panel.
Thick contours represent $\mathrm{Kn}_{\theta}$ at $\eta=0$, while thin
contours represent $\eta=2$.

The rapidity dependence itself of the shifted solution can be seen in Fig. 7,
where we compare contours of constant $\mathrm{Kn}_{\theta}$ for $\eta=0$
(thick lines) with $\eta=2$ (thin lines). We emphasize that the solution is
independent of $\eta$ (i.e., boost invariant) when $t_{0}=0$. For $t_{0}>0$,
however, we see in Fig. 7(a) that increasing $\eta$ away from zero shifts
freeze out to earlier $\tau$, without dramatically affecting the radial
dependence. The radial velocity (Fig. 7(c)) likewise shifts to smaller values
without modifying the radial dependence significantly, producing a roughly 8%
in the maximum $v_{r}/c$ at $\eta=2$ relative to mid-rapidity when
$\mathrm{Kn}_{\theta}=0.75$. In Fig. 7(b), we observe that, since both $T$ and
$\mathrm{Kn}$ are Lorentz scalars, their relationship is unaffected by the
shift in $t$, and thus the thick and thin lines coincide for all
$\mathrm{Kn}_{\theta}$ contours, regardless of the values of $t_{0}$ and
$\eta$. Furthermore, increasing the magnitude of $\eta$ further shifts freeze
out to still earlier times. Taken together, these results imply once again
that an accurate description of the system’s shape and flow are both essential
for a quantitatively precise characterization of the freeze out process in
nuclear collisions.

Our shifted version of Gubser’s original solution thus illustrates the
importance of accounting for both transverse and longitudinal flow (and, more
precisely, their gradients) on the shape of the freeze out surface. In
particular, the constant temperature contours and constant Knudsen contours
imply qualitatively similar geometries for the freeze out hypersurface when
the flow is relatively weak. However, in the presence of strong collective
flow, such as one might expect to find in small collision systems Kalaydzhyan
and Shuryak (2015); Nagle and Zajc (2018); Heinz and Moreland (2019), which
freeze out criterion one adopts is likely to have a larger systematic effect
on experimental observables. The effects could be especially pronounced for
observables which are directly sensitive to the geometry of the freeze out
hypersurface, such as the Hanbury Brown – Twiss interferometric radii Plumberg
(2020a, b). We plan to explore the nature of this sensitivity in a subsequent
study.

##  V Conclusions

In this work, we have developed two novel exact solutions conformal
hydrodynamics. Both solutions apply a constant shift of the Minkowski time
coordinate (a “temporal shift") to induce a non-trivial rapidity dependence
into an already known solution. In the first case (our “Solution I"), we have
generalized Gubser’s solution of ideal relativistic hydrodynamics to allow for
flow in the rapidity direction and combined this with a temporal shift to
generate a solution with rapidity dependence as well. In the second case
(“Solution II"), we have applied a temporal shift to Gubser’s original
solution of viscous relativistic hydrodynamics, yielding a solution for which
we can define out-of-equilibrium quantities like the Knudsen number. Our new
solutions admit a number of interesting applications, including providing
useful benchmark tests for 3+1D hydrodynamic codes, as well as allowing
analytically tractable insights into the nature of far-from-equilibrium
dynamics in nuclear collisions.

As a concrete instance of this latter application, we have used the Solution
II constructed here to clarify two related questions: first, the amount by
which alternative freeze out criteria may differ in the freeze out
hypersurfaces which they produce, and second, the extent to which these
differences are influenced by the strength of collective flow and the size of
its gradients throughout the system.

We have answered the first question by considering the qualitative differences
which emerge in our solution between contours of fixed temperature $T$ and
fixed Knudsen number $\mathrm{Kn}$. We have found that large discrepancies
tend to arise as the system falls progressively away from local equilibrium:
in our solution, this tends to occur at temperatures $T\lesssim 100$ MeV. We
anticipate that these differences should lead to systematic effects on the
implementation of the freeze out process and the evaluation of, e.g., Cooper-
Frye integrals Cooper and Frye (1974); Cooper  _et al._ (1975) in more
realistic numerical simulations.

To address the second question, we directly plotted the local temperature $T$
and radial velocity $v_{r}$ along contours of constant $\mathrm{Kn}$ using our
Solution II. This allowed us to quantify the typical magnitude of these
effects and to explore the role played by collective flow on the
discrepancies. We have also considered how these results are affected by the
artificial suppression of $u^{\eta}$ in the full solution, as a way of
assessing their sensitivity to the flow itself. We found that our solution
predicts a strong temperature variation (of up to 72%) along constant
$\mathrm{Kn}$ contours in regions where the flow is sufficiently strong. Even
in regions where the $T$ and $\mathrm{Kn}$ contours were in qualitative
agreement, we still observed up to 13% variation in $T$ (relative to
$\mathrm{Kn}$) in those regions where the bulk of particle production occurs
Schnedermann  _et al._ (1993). This underscores the fact that systematic
effects can and do arise on the basis of which freeze out criterion one adopts
in hydrodynamic simulations.

Since the discrepancies we observe between different freeze out criteria are
somewhat more pronounced than those seen in previous studies (e.g., Niemi and
Denicol (2014); Ahmad  _et al._ (2017)), it is important to determine whether
the results obtained here can be reproduced in the context of more realistic
numerical simulations and to assess what, if any, are the systematic effects
on predictions for experimental observables. Doing so could shed valuable
light on the physical mechanisms underlying the robustness of fluid dynamical
behavior Heller and Spalinski (2015); Heinz and Moreland (2019) and help to
identify specific observables which are especially well suited to probing the
breakdown of hydrodynamics. In this respect, HBT interferometry is an example
of an observable which should be ideal for probing the importance of these
effects for the geometry of the freeze out hypersurface, and we plan to
explore its sensitivity to the freeze out criterion in a subsequent study.

###### Acknowledgements.

The authors sincerely thank Travis Dore, Ulrich Heinz, and Jorge Noronha for
several useful comments and suggestions for improvements to the manuscript.

## References

  * Gale  _et al._ (2013) C. Gale, S. Jeon, and B. Schenke, Int. J. Mod. Phys. A 28, 1340011 (2013), arXiv:1301.5893 [nucl-th] .
  * Jeon and Heinz (2015) S. Jeon and U. Heinz, Int. J. Mod. Phys. E 24, 1530010 (2015), arXiv:1503.03931 [hep-ph] .
  * Romatschke and Romatschke (2019) P. Romatschke and U. Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Nagle and Zajc (2018) J. L. Nagle and W. A. Zajc, Ann. Rev. Nucl. Part. Sci. 68, 211 (2018), arXiv:1801.03477 [nucl-ex] .
  * Adolfsson  _et al._ (2020) J. Adolfsson _et al._ , Eur. Phys. J. A 56, 288 (2020), arXiv:2003.10997 [hep-ph] .
  * Song  _et al._ (2011) H. Song, S. A. Bass, U. Heinz, T. Hirano, and C. Shen, Phys. Rev. Lett. 106, 192301 (2011), [Erratum: Phys.Rev.Lett. 109, 139904 (2012)], arXiv:1011.2783 [nucl-th] .
  * Bozek and Wyskiel-Piekarska (2012) P. Bozek and I. Wyskiel-Piekarska, Phys. Rev. C 85, 064915 (2012), arXiv:1203.6513 [nucl-th] .
  * Gardim  _et al._ (2012) F. G. Gardim, F. Grassi, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. Lett. 109, 202302 (2012), arXiv:1203.2882 [nucl-th] .
  * Bozek and Broniowski (2013) P. Bozek and W. Broniowski, Phys. Rev. C 88, 014903 (2013), arXiv:1304.3044 [nucl-th] .
  * Niemi  _et al._ (2016a) H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016a), arXiv:1505.02677 [hep-ph] .
  * Ryu  _et al._ (2015) S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * McDonald  _et al._ (2017) S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, and C. Gale, Phys. Rev. C 95, 064913 (2017), arXiv:1609.02958 [hep-ph] .
  * Bernhard  _et al._ (2016) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Gardim  _et al._ (2017) F. G. Gardim, F. Grassi, M. Luzum, and J. Noronha-Hostler, Phys. Rev. C 95, 034901 (2017), arXiv:1608.02982 [nucl-th] .
  * Alba  _et al._ (2018) P. Alba, V. Mantovani Sarti, J. Noronha, J. Noronha-Hostler, P. Parotto, I. Portillo Vazquez, and C. Ratti, Phys. Rev. C 98, 034909 (2018), arXiv:1711.05207 [nucl-th] .
  * Giacalone  _et al._ (2018) G. Giacalone, J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 97, 034904 (2018), arXiv:1711.08499 [nucl-th] .
  * Eskola  _et al._ (2018) K. J. Eskola, H. Niemi, R. Paatelainen, and K. Tuominen, Phys. Rev. C 97, 034911 (2018), arXiv:1711.09803 [hep-ph] .
  * Weller and Romatschke (2017) R. D. Weller and P. Romatschke, Phys. Lett. B 774, 351 (2017), arXiv:1701.07145 [nucl-th] .
  * Schenke  _et al._ (2019) B. Schenke, C. Shen, and P. Tribedy, Phys. Rev. C 99, 044908 (2019), arXiv:1901.04378 [nucl-th] .
  * Everett  _et al._ (2021) D. Everett _et al._ (JETSCAPE), Phys. Rev. Lett. 126, 242301 (2021), arXiv:2010.03928 [hep-ph] .
  * Noronha-Hostler  _et al._ (2016) J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 93, 034912 (2016), arXiv:1511.06289 [nucl-th] .
  * Niemi  _et al._ (2016b) H. Niemi, K. J. Eskola, R. Paatelainen, and K. Tuominen, Phys. Rev. C 93, 014912 (2016b), arXiv:1511.04296 [hep-ph] .
  * Adam  _et al._ (2016) J. Adam _et al._ (ALICE), Phys. Rev. Lett. 116, 222302 (2016), arXiv:1512.06104 [nucl-ex] .
  * Heinz (2004) U. W. Heinz, in  _2nd CERN-CLAF School of High Energy Physics_ (2004) pp. 165–238, arXiv:hep-ph/0407360 .
  * Schnedermann  _et al._ (1993) E. Schnedermann, J. Sollfrank, and U. W. Heinz, NATO Sci. Ser. B 303, 175 (1993).
  * Petersen  _et al._ (2008) H. Petersen, J. Steinheimer, G. Burau, M. Bleicher, and H. Stöcker, Phys. Rev. C 78, 044901 (2008), arXiv:0806.1695 [nucl-th] .
  * Denicol  _et al._ (2018) G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 98, 034916 (2018), arXiv:1804.10557 [nucl-th] .
  * Plumberg  _et al._ (2023a) C. Plumberg _et al._ , in  _30th International Conference on Ultrarelativstic Nucleus-Nucleus Collisions_ (2023) arXiv:2312.07415 [hep-ph] .
  * Eskola  _et al._ (2008) K. J. Eskola, H. Niemi, and P. V. Ruuskanen, Phys. Rev. C 77, 044907 (2008), arXiv:0710.4476 [hep-ph] .
  * Holopainen and Huovinen (2012) H. Holopainen and P. Huovinen, J. Phys. Conf. Ser. 389, 012018 (2012), arXiv:1207.7331 [hep-ph] .
  * Niemi and Denicol (2014) H. Niemi and G. S. Denicol, (2014), arXiv:1404.7327 [nucl-th] .
  * Ahmad  _et al._ (2017) S. Ahmad, H. Holopainen, and P. Huovinen, Phys. Rev. C 95, 054911 (2017), arXiv:1608.03444 [nucl-th] .
  * Huovinen and Holopainen (2017) P. Huovinen and H. Holopainen, Acta Phys. Polon. Supp. 10, 913 (2017).
  * Heinz and Moreland (2019) U. W. Heinz and J. S. Moreland, J. Phys. Conf. Ser. 1271, 012018 (2019), arXiv:1904.06592 [nucl-th] .
  * Summerfield  _et al._ (2021) N. Summerfield, B.-N. Lu, C. Plumberg, D. Lee, J. Noronha-Hostler, and A. Timmins, Phys. Rev. C 104, L041901 (2021), arXiv:2103.03345 [nucl-th] .
  * Plumberg  _et al._ (2023b) C. Plumberg, D. Almaalol, T. Dore, J. Noronha, and J. Noronha-Hostler, EPJ Web Conf. 276, 06002 (2023b).
  * Plumberg (2020a) C. Plumberg, Phys. Rev. C 102, 054908 (2020a), arXiv:2008.01709 [nucl-th] .
  * Bjorken (1983) J. D. Bjorken, Phys. Rev. D 27, 140 (1983).
  * Gubser (2010) S. S. Gubser, Phys. Rev. D 82, 085027 (2010), arXiv:1006.0006 [hep-th] .
  * Gubser and Yarom (2011) S. S. Gubser and A. Yarom, Nucl. Phys. B 846, 469 (2011), arXiv:1012.1314 [hep-th] .
  * Hatta  _et al._ (2014a) Y. Hatta, J. Noronha, and B.-W. Xiao, Phys. Rev. D 89, 051702 (2014a), arXiv:1401.6248 [hep-th] .
  * Hatta  _et al._ (2014b) Y. Hatta, J. Noronha, and B.-W. Xiao, Phys. Rev. D 89, 114011 (2014b), arXiv:1403.7693 [hep-th] .
  * Csörgő  _et al._ (2018) T. Csörgő, G. Kasza, M. Csanád, and Z. Jiang, Universe 4, 69 (2018), arXiv:1805.01427 [nucl-th] .
  * Shi  _et al._ (2022) S. Shi, S. Jeon, and C. Gale, Phys. Rev. C 105, L021902 (2022), arXiv:2201.06670 [hep-ph] .
  * Marrochio  _et al._ (2015) H. Marrochio, J. Noronha, G. S. Denicol, M. Luzum, S. Jeon, and C. Gale, Phys. Rev. C 91, 014903 (2015), arXiv:1307.6130 [nucl-th] .
  * Weil  _et al._ (2016) J. Weil _et al._ (SMASH), Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .
  * Denicol  _et al._ (2012) G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Phys. Rev. D 85, 114047 (2012), [Erratum: Phys.Rev.D 91, 039902 (2015)], arXiv:1202.4551 [nucl-th] .
  * Molnár  _et al._ (2014) E. Molnár, H. Niemi, G. S. Denicol, and D. H. Rischke, Phys. Rev. D 89, 074010 (2014), arXiv:1308.0785 [nucl-th] .
  * Plumberg (2020b) C. Plumberg, (2020b), arXiv:2010.11957 [nucl-th] .
  * Israel and Stewart (1979) W. Israel and J. M. Stewart, Annals Phys. 118, 341 (1979).
  * Kalaydzhyan and Shuryak (2015) T. Kalaydzhyan and E. Shuryak, Phys. Rev. C 91, 054913 (2015), arXiv:1503.05213 [hep-ph] .
  * Cooper and Frye (1974) F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
  * Cooper  _et al._ (1975) F. Cooper, G. Frye, and E. Schonberg, Phys. Rev. D 11, 192 (1975).
  * Heller and Spalinski (2015) M. P. Heller and M. Spalinski, Phys. Rev. Lett. 115, 072501 (2015), arXiv:1503.07514 [hep-th] .