11institutetext: Center for Nuclear Theory, Department of Physics and
Astronomy,
Stony Brook University, Stony Brook, New York 11794, USA

# Hydrodynamic attractors in heavy ion collisions: a review

Alexander Soloviev

(Received: date / Revised version: date)

###### Abstract

A review of the recent progress of relativistic hydrodynamic attractors is
presented, with a focus on applications in heavy ion collisions and the quark
gluon plasma. Pedagogical introductions to the effective descriptions relevant
for attractors in high energy physics, namely hydrodynamics, holography and
kinetic theory, are followed by highlights of some recent advances.

##  1 Introduction

Relativistic hydrodynamics Kovtun:2012rj ; Romatschke:2017ejr ; Shen:2020mgh
is an incredibly useful tool in effective theories, capturing the underlying
symmetries of a system in the long wavelength limit. It has been successfully
applied at a vast swathe of scales: from the smallest in understanding quantum
chromodynamics (QCD) at high energies Heinz:2013th , to condensed matter
systems chaikin2000principles  and to the largest in cosmological applications
shore ; Weinberg:2008zzc .

The application of hydrodynamics is not always readily apparent. In fact, in
the high energy context one of the most unexpected surprises to come from
evaluating the data at the Relativistic Heavy Ion Collider (RHIC) and the
Large Hadron Collider (LHC) was that the quark-gluon plasma (QGP), the extreme
state of matter composed from liberated constituents of neutrons and protons
following collisions of heavy ions, behaves like the “most perfect fluid” and
as such is amenable to a description via relativistic hydrodynamics. Abundant
evidence points to the QGP exhibiting collective flow, such as the observation
of elliptic flow Teaney:2009qa ; Snellings:2011sz  and the particularly low
value of the ratio of the shear viscosity to entropy density Romatschke:2007mq
, close to the value predicted from holography Policastro:2001yc ;
Kovtun:2003wp ; Kovtun:2004de . For further mechanisms and features of the
thermalization of the QGP, see Berges:2020fwq  and references therein.

As such, the hydrodynamic description of the QGP had great phenomenological
success. Active research deals with the hydrodynamization, i.e. the approach
of the post-collision state to a fluid description, see Florkowski:2017olj
for a review. It is worth reiterating how surprising this is since
hydrodynamics is an effective theory based on the standard assumption that the
system is near thermal equilibrium, although what exactly sets the regime of
applicability of hydrodynamics was not clear. In typical heavy ion collisions,
the system is certainly not in or even near equilibrium, which naturally
beckons the question: why is hydrodynamics so unreasonably effective?

Recently, there has been progress in understanding how systems with little
hope of being in equilibrium become nonetheless well-approximated by
hydrodynamics. One proposed solution is due to the emergence of attractor
phenomena. In a word, the study of attractors stems from the observation that
solutions of evolution equations with different initial conditions evolve to
the same late time behavior, which is called the (hydrodynamic) attractor.
Quite generally, nonequilibrium solutions of different microscopic theories
decay to attractors, losing knowledge of their initial conditions well before
reaching thermal equilibrium. Specifically, heavy ion collisions have a mix of
hydrodynamic and non-hydrodynamic modes, the latter of which decay rapidly in
the expanding medium, leaving behind just the long wavelength modes, well
described by hydrodynamics. For the heavy ion community, attractors came to
the forefront with the remarkable discovery of the hydrodynamic attractor in
an expanding boost-invariant background by Heller and Spalinski Heller:2015dha
, building on the earlier results of Borel resummed hydrodynamic gradient
expansion in holographic Bjorken flow Heller:2011ju ; Heller:2013fn .

Spurred on by the characterization of the hydrodynamic attractors, attention
was drawn to attractors in different descriptions relevant for heavy ion
collisions. For instance, similar attractors were found in kinetic theory and
in holography, even in non-conformal and non-homogeneous settings
Florkowski:2017jnz ; Romatschke:2017acs , which indicates that attractor
solutions should exist for a wide range of theories independent of the
symmetries of those theories. In hindsight, this should not be so surprising,
since the hydrodynamic limit is found in so many diverse physical settings.
For instance, hydrodynamic attractors have also been explored in holographic
models of QGP Buchel:2016cbj  and in Gauss-Bonnet gravity Casalderrey-
Solana:2017zyh , which could have well been anticipated due to the fluid-
gravity correspondence Rangamani:2009xk .

A key feature of the attractor is that when it is present, the system
“forgets” about its initial conditions. Therefore, attractors can be used
phenomenologically to distinguish between observables which are dependent on
initial conditions and those which are not. This can lead to the
identification of universal behavior in heavy ion collisions. Recently, the
potential of hydrodynamic attractors as a tool for phenomenological studies
has begun to be appreciated. For instance, they have been used to relate
initial state energies to the produced particle multiplicity Giacalone:2019ldn
; Jankowski:2020itt . Moreover, recent studies examined the hydrodynamic
attractor at finite chemical potential Dore:2020fiq ; Dore:2020jye , tracing
the hydrodynamic evolution as trajectories through the QCD phase diagram. One
of the more exciting ideas relevant to heavy ion physics to come to fruition
at least in part due to attractors is the Kurkela-Mazeliauskas-Paquet-
Schlichtling-Teaney (KøMPøST) framework Kurkela:2018wud ; Kurkela:2018vqr ,
essentially linking the early time dynamics, captured by a kinetic
description, to the late time hydrodynamic behavior. Another practical
application of hydrodynamic attractors is as a viable arena to study
resurgence Basar:2015ava ; Aniceto:2018bis . A comprehensive discussion of
resurgence is beyond the scope of the present review.

Of course, it would be remiss to not mention that attractors have a rich,
storied history in dynamical systems. There is considerable value in trying to
understand the attractor studied in heavy ion collisions in the more
mathematical language of dynamical systems. Progress in this direction has
already been made, see Behtash:2017wqg ; Behtash:2019txb ; Behtash:2019qtk .

A related story to the hydrodynamic attractor is the study of non-thermal
fixed points Berges:2013eia ; Berges:2013fga ; Berges:2013lsa  in heavy ion
physics (as well as in a variety of other contexts Berges:2014bba ), which are
essentially an observation that at early enough times, the distribution
function exhibits a special scaling behavior. This early time regime is
currently the subject of numerical study, such as in a simplified model of the
glasma in $2+1$ dimensions Boguslavski:2019fsb . This should be considered as
a distinct feature of the weakly coupled degrees of freedom at early times in
the QGP, different from the hydrodynamic attractors that we focus on here.

Since hydrodynamic attractors have been found in a wide range of (effective)
theories pertinent to heavy ion collisions, many works have focussed on
comparing features of attractor behavior in relativistic hydrodynamics,
kinetic theory and holography Romatschke:2017acs ; Kurkela:2019set . These
theories vary not only in their frameworks, but also in their description of
weak and strong coupling regimes. The benefits of this approach are clear,
giving ample evidence of the importance of the attractors. This review will
present in pedagogical detail the attractor story for the trio of theories
mentioned previously, first by describing an illustrative example, followed by
reviewing current progress in the field.

As such, the organization of the paper is as follows. In Sec. 2, we will
discuss the prototypical example of the hydrodynamic attractor in the simplest
setting of $(0+1)$-dimensional Bjorken expansion, before reviewing further
progress. Next, in Sec. 3, we will present the attractor as it pertains to
holographic systems. Then, we will turn our attention to attractors in kinetic
theory in Sec. 4, especially focussing on the relaxed time approximation.
Finally, in Sec. 5, we turn to the recent applications of attractors to heavy
ion collisions. A table of useful abbreviations is included in Appendix A.

##  2 Hydrodynamic attractors

In this section, we will review the advances in studying the attractor in
expanding relativistic hydrodynamics. As is usually the case, it is best to
start simple and then discuss more complicated examples. We will refresh the
reader with a discussion of relativistic hydrodynamics and various models of
curing the acausality, before reaching the simplest and quintessential example
of an attractor in $0+1d$ Müller-Israel-Stewart (MIS) Bjorken flow. We will
finish by going into some detail of some advances in the field.

###  2.1 Primer on relativistic hydrodynamics

Hydrodynamics represents the course-grained description of microscopic
dynamics. The macroscopic objects of interest are the energy momentum tensor,
$T^{\mu\nu}$, and the conserved currents, $J^{\mu}$. The dynamical equations
are generated by considering the conservation of the energy momentum tensor,

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

as well as the conservation of the currents, schematically
$\nabla_{\mu}J^{\mu}=0$.

Generally, the fluid energy momentum tensor can be decomposed into ideal and
dissipative contributions

$\displaystyle T^{\mu\nu}=T^{\mu\nu}_{\text{ideal}}+\Pi^{\mu\nu}.$ (2)

The ideal energy-momentum tensor is

$\displaystyle T^{\mu\nu}_{\text{ideal}}=\varepsilon
u^{\mu}u^{\nu}+P\Delta^{\mu\nu},$ (3)

where $\varepsilon$ is the energy density, $P$ is the pressure and
$\Delta^{\mu\nu}=u^{\mu}u^{\nu}+g^{\mu\nu}$ is the spatial projection
operator. Note that the four-velocity, $u^{\mu}$, is timelike normalized,
$u^{\mu}u_{\mu}=-1$. The dissipative tensor, $\Pi^{\mu\nu}$, generically
contains higher derivative corrections of $u^{\mu}$ and the temperature, $T$.
This organization is known as the gradient expansion.

We will ultimately need to specify a reference frame for the fluid. Here, we
choose to work in the Landau frame, where the energy density is the eigenvalue
associated with $u^{\mu}$ via $T^{\mu\nu}u_{\nu}=-\varepsilon u^{\mu}$. This
leads to the requirement that the dissipation tensor is orthogonal to the
flow, i.e.

$\displaystyle\Pi^{\mu\nu}u_{\mu}=0.$ (4)

At this point, the standard approach would be to specify the dissipation
tensor order by order in a gradient expansion with arbitrary transport
coefficients. These transport coefficients are part of the required input for
a complete hydrodynamic description, which can be determined from the
underlying microscopic theory. To first order in derivatives, we have two
distinct tensor structures, which we write as follows

$\displaystyle\Pi^{\mu\nu}=-\eta\sigma^{\mu\nu}-\zeta\Delta^{\mu\nu}\nabla_{\alpha}u^{\alpha},$
(5)

where the two transport coefficients are the shear viscosity, $\eta$, and the
bulk viscosity, $\zeta$. Higher orders lead to more possible tensor structures
with additional transport coefficients. The shear tensor is transverse,

$\displaystyle\sigma^{\mu\nu}u_{\mu}=0,$ (6)

traceless,

$\displaystyle g_{\mu\nu}\sigma^{\mu\nu}=0,$ (7)

and is explicitly given by

$\displaystyle\sigma^{\mu\nu}=\Delta^{\mu\alpha}\Delta^{\nu\beta}(\nabla_{\alpha}u_{\beta}+\nabla_{\beta}u_{\alpha})-\frac{2}{3}\Delta^{\mu\nu}\nabla_{\alpha}u^{\alpha}.$
(8)

The causal hydrodynamic evolution of a system can be described via (1). There
is one complication: the four equations in (1) contain five variables: the
energy density, $\varepsilon$, the pressure, $P$, and the four velocity,
$u^{\mu}$. As mentioned before, the last needs to satisfy the timelike
normalization condition $u^{\mu}u_{\mu}=-1,$ meaning that $u^{\mu}$ has three
independent components. Thus, we are shy one equation.

To fill this gap, we will need to specify an equation of state, relating the
pressure to the energy density, $P=P(\varepsilon)$. In general, the equation
of state depends on the system that one is interested in modelling. One of the
simplest choices (which we will largely adopt here) is the conformal equation
of state $\varepsilon=3P$, which means that the energy momentum tensor is
traceless,

$\displaystyle g_{\mu\nu}T^{\mu\nu}=0.$ (9)

In the conformal case, the energy density can only depend on temperature,
which means that $\varepsilon\sim T^{4}.$

####  2.1.1 Truncation schemes

The issue with the choice of (5) for the dissipation tensor is that resulting
equations become acausal Hiscock:1985zz ; Kostadt:2000ty . A straightforward
way to see this, e.g. Romatschke:2017ejr , is by considering metric
perturbations of the form, $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$, and seeing
how the stress tensor (2) changes due to the metric perturbation in linear
response. In essence, we are computing the retarded correlator,
$G^{{\mu\nu},{\alpha\beta}}_{R}$, via

$\displaystyle
T^{\mu\nu}=T_{0}^{\mu\nu}-\frac{1}{2}G^{{\mu\nu},{\alpha\beta}}_{R}h_{\alpha\beta}+\ldots$
(10)

where the dots obscure terms higher order in perturbations. Focussing simply
on the shear channel, we see that the retarded correlator is Kovtun:2012rj

$\displaystyle G^{0x,0x}_{R}=\varepsilon+\frac{\eta
k^{2}}{i\omega-\frac{\eta}{\varepsilon+P}k^{2}}.$ (11)

From the denominator, we identify the diffusive pole,
$\omega=-i\frac{\eta}{\varepsilon+P}k^{2}$. The group velocity,
$\left|{d\omega}/{dk}\right|$, is given by

$\displaystyle\left|\frac{d\omega}{dk}\right|=2\frac{\eta}{\varepsilon+P}k,$
(12)

which grows without bound as $k$ gets large.

There are a few ways to cure such acausality. A simple and popular approach is
to implement the MIS equations Muller:1967zza ; Israel:1976tn ; Israel:1979wp
(for more recent developments, see Baier-Romatschke-Son-Starinets-Stephanov
(BRSSS) Baier:2007ix ), where the dissipation tensor satisfies a relaxation
equation. Here, as we will be interested in the conformal case of non-zero
shear viscosity, but vanishing bulk viscosity (as well as higher order
transport coefficients), the MIS equation in this case is then

$\displaystyle\left(\tau_{\pi}u^{\alpha}\nabla_{\alpha}+1\right)\Pi^{\mu\nu}=-\eta\sigma^{\mu\nu},$
(13)

where $\tau_{\pi}$ is the relaxation time. The previous case (5) corresponds
to vanishing relaxation time $\tau_{\pi}$ (as well as $\zeta=0$).

How does this cure the acausality? For small only time-dependent
perturbations, we see that

$\displaystyle\delta\Pi^{\mu\nu}=\frac{\eta}{i\tau_{\pi}\omega-1}\delta\sigma^{\mu\nu},$
(14)

which can be compared directly to (5) in the limit that
$\omega,\zeta\rightarrow 0$. This simply amounts to replacing
$\eta\rightarrow\frac{\eta}{i\tau_{\pi}\omega-1}$. Revisiting the computation
for the retarded correlator, we see that the diffusive mode is found when

$\displaystyle
i\omega(1-i\tau_{\pi}\omega)-\frac{\eta}{\varepsilon+P}k^{2}=0.$ (15)

Solving for $\omega$ and computing the group velocity, we find one of the
modes has dispersion

$\displaystyle\left|\frac{d\omega}{dk}\right|=\frac{2\frac{\eta}{\varepsilon+P}k}{\left|\sqrt{\frac{4\eta\tau_{\pi}}{\varepsilon+P}k^{2}-1}\right|}.$
(16)

The group velocity now saturates to a maximum value, namely

$\displaystyle\lim_{k\rightarrow\infty}\left|\frac{d\omega}{dk}\right|=\sqrt{\frac{\eta}{\tau_{\pi}(\varepsilon+P)}}$
(17)

and as such, there is no issue with causality, so long as the right hand side
of (17) is smaller than 1.

A more recent improvement to the MIS equations is known as Denicol-Niemi-
Molnar-Rischke (DNMR) effective theory Denicol:2012cn , which is a second
order hydrodynamic theory using Grad’s 14-moment approximation for the single
particle distribution function. To first order in gradients, DNMR is similar
to MIS (13), in that (taking $u^{\mu}=(1,0,0,0)$) we see that the general form
the equations take to this order is Strickland:2017kux

$\displaystyle\partial_{\tau}{\Pi}=\frac{4\eta}{3\tau\tau_{\pi}}-\beta_{\pi\pi}\frac{\Pi}{\tau}-\frac{\Pi}{\tau_{\pi}},$
(18)

where for DNMR we have $\beta_{\pi\pi}=38/21$ and $\tau_{\pi}=\tau_{eq}$
Denicol:2012cn , while for MIS this is $\beta_{\pi\pi}=4/3$ and
$\tau_{\pi}=6\tau_{\rm eq}/5$ Muronga:2003ta .

Another truncation scheme goes under the name of anisotropic hydrodynamics
(aHydro) Strickland:2014pga ; McNelis:2018jho ; Almaalol:2018ynx ;
Alalawi:2020zbx . Essentially, the idea is to formulate hydrodynamics in
momentum space with a view towards taking into account the large anisotropies
found in the QGP. One takes a distribution function with an explicit
anisotropy parameter, $\xi$ (see Sec. 4 for a primer on kinetic theory and
(163) for the form of the anisotropic distribution function). As outlined in
Sec. 4, the Landau matching condition determines the temperature. The
evolution of $\xi$ is determined by considering a higher moment of the
Boltzmann equation. In this way, aHydro resums an infinite number of terms in
an expansion in inverse Reynolds number Strickland:2017kux

$\displaystyle
R^{-1}=\frac{\sqrt{\Pi^{\mu\nu}\Pi_{\mu\nu}}}{P_{0}}\propto\bar{\Pi}(\xi),$
(19)

where $\bar{\Pi}=\Pi/\varepsilon.$ A comparison to the other truncation
schemes can be made by rewriting the evolution equation for $\xi$ as a
function of $\Pi$ and expanding for small anisotropy. This results in
agreement with (18) to linear order in $\Pi$. However, note that at higher
order $\Pi^{2}$, the term $\lambda_{1}\Pi^{2}/(2\tau_{\pi}\eta^{2})$ has a
numerical disagreement of the transport coefficient, $\lambda_{1}^{\rm
aHydro}=\eta\tau_{\pi}/7$, with $\lambda_{1}^{\rm DNMR}=5\lambda_{1}^{\rm
aHydro}$ Romatschke:2011qp ; Strickland:2017kux . The differences between the
two approaches arises from the fact that aHydro includes higher order terms
compared to DNMR at the same level of truncation.

Which scheme to use? In terms of conceptual simplicity, the MIS framework
tends to get more mileage since it is rather straightforward to include
additional higher order transport coefficients. MIS and DNMR differ mainly in
their far-from-equilibrium transient effects, as they have the same late time
Navier-Stokes limit Dore:2020jye . It was found in Strickland:2017kux  that
DNMR is more accurate than MIS when compared to the exact solution in the
relaxed time approximation (RTA), as it captures systematically second order
hydrodynamic contributions, while aHydro approximates the exact attractor the
best. Likewise, in the context of Gubser flow where the exact attractor is
known, aHydro did a better job in capturing the attractor compared to MIS and
DNMR Heinz:2015gka ; Behtash:2017wqg .

As an aside, it is worthwhile to underline that the MIS equations represent a
particularly easy model to work with. However, if one were inclined to include
non-hydrodynamic modes to their model building, such as the oscillatory
behavior of $N=4$ super Yang-Mills theory via a hydrodynamic model, the
Heller-Janik-Spalinski-Witaszczyk (HJSW) model could be considered a natural
generalization to the MIS equations Heller:2014wfa ; Aniceto:2015mto ;
Heller:2020anv . In this case, the dissipation tensor, $\Pi^{\mu\nu}_{\rm
total}=\Pi^{\mu\nu}_{\rm MIS}+\Pi^{\mu\nu}_{\rm HSJW}$, can be decomposed into
a piece evolving with the MIS equations (13) and another piece
satisfying111Note that for the purposes of modelling, one can also study the
case where $\Pi^{\mu\nu}_{\rm MIS}=0.$

$\displaystyle\left((\frac{1}{T}u^{\alpha}\nabla_{\alpha})^{2}+2\Omega_{I}u^{\alpha}\nabla_{\alpha}+|\Omega|^{2}\right)\Pi^{\mu\nu}_{\rm
HSJW}=-\eta|\Omega|^{2}\sigma^{\mu\nu}-\frac{c_{\sigma}}{T}u^{\mu}\nabla_{\mu}(\eta\sigma^{\mu\nu})+\ldots$
(20)

where $\Omega=\Omega_{R}+i\Omega_{I}$ is the complex frequency, meant to mimic
the quasinormal mode frequency. The values of the real and imaginary part can
be fixed by comparing to $N=4$ Super Yang-Mills (SYM) theories Nunez:2003eq ;
Aniceto:2015mto . The presence of these non-hydrodynamic modes is guaranteed
by the appearance of the second derivative in the above. In principle, one can
further generalize by considering higher derivatives on the left hand side of
(20), although how to specify the extra initial conditions is not clear.

###  2.2 $0+1d$ Bjorken flow

Now let’s get our hands dirty with an example computation. We will review the
hydrodynamic attractor in Bjorken flow Heller:2015dha , a $(0+1)$-dimensional
expanding transversally homogeneous system with longitudinal boost-invariance.

Bjorken flow Bjorken:1982qr  describes a simple model for heavy ion
collisions, where the collision axis is along the z-axis and the system is
otherwise homogeneous. Essentially, the nuclei are modelled as infinite flat
sheets in the transverse ($x,y$) directions. Furthermore, the matter produced
in the forward light cone is boost invariant. The background metric is given
in Milne coordinates

$\displaystyle g_{\mu\nu}=\text{diag}(-1,1,1,\tau^{2}).$ (21)

We can relate Milne coordinates $(\tau,x,y,\xi)$ to Minkowski coordinates
$(t,x,y,z)$ by observing that the proper time and the spacetime rapidity are

$\displaystyle\tau$ $\displaystyle=\sqrt{t^{2}-z^{2}},$ (22)
$\displaystyle\xi$ $\displaystyle=\tanh^{-1}(z/t),$ (23)

respectively. Note that the metric (21) describes a space which is Ricci flat.
For future convenience, we also note that the nonzero Christoffel symbols are

$\displaystyle\Gamma^{\tau}_{\phantom{\tau}\xi\xi}=\tau,\quad\Gamma^{\xi}_{\phantom{\xi}\xi\tau}=\tau^{-1}.$
(24)

We will assume all variables depend only on the proper time, $\tau$, which
means that (1) and (13) take a particularly simple form:

$\displaystyle 0$
$\displaystyle=\tau\dot{\varepsilon}+\frac{4}{3}\varepsilon-\phi,$ (25)
$\displaystyle 0$
$\displaystyle=\tau_{\pi}\dot{\phi}+\frac{4\tau_{\pi}\phi}{3\tau}+\phi-\frac{4\eta}{3\tau},$
(26)

where $\phi\equiv-\Pi^{\xi}_{\xi}$. Thus, the conservation equation, (25), and
the MIS equation, (26), represent the set of equations describing the
relativistic evolution of a viscous fluid in a Bjorken expanding background.
We proceed by expressing the shear viscosity and the relaxation time in terms
of the energy density:

$\displaystyle\eta=C_{\eta}\frac{\varepsilon+P}{\varepsilon^{1/4}},\quad\quad\tau_{\Pi}=C_{\tau_{\Pi}}\varepsilon^{-1/4},$
(27)

where $C_{\eta}$ and $C_{\tau_{\Pi}}$ are dimensionless constants. As we will
see shortly, the relaxation time is the hydrodynamic transport coefficient
which indicates the decay time of the non-hydrodynamic mode. Taking insight
from holography Baier:2007ix , we will take these constants to be

$\displaystyle C_{\eta}=\frac{1}{4\pi},\quad C_{\tau_{\Pi}}=\frac{2-\ln
2}{2\pi}.$ (28)

We can then eliminate $\phi$ to arrive to arrive at an expression solely for
the energy density

$\displaystyle
9C_{\tau_{\Pi}}\tau^{2}\ddot{\varepsilon}+(33C_{\tau_{\Pi}}\tau+9\tau^{2}\varepsilon^{1/4})\dot{\varepsilon}+16(C_{\tau_{\Pi}}-C_{\eta})\varepsilon+12\tau\varepsilon^{5/4}=0.$
(29)

We can recover the expressions found in Heller:2015dha  by remembering that in
a conformal theory the energy density is related to the temperature via
$\varepsilon\sim T^{4}$.

Figure 1: Hydrodynamic attractor for $0+1d$ Bjorken flow. Solutions of (31)
with various initial conditions are shown in blue, while the attractor is
shown in magenta. The red and green dashed lines correspond to first and
second order hydrodynamics. Figure from Heller:2015dha .

Following Heller:2011ju ; Heller:2015dha , we can now introduce dimensionless
quantities

$\displaystyle\omega=\tau\varepsilon^{1/4},\quad
f=\tau\partial_{\tau}\log{\omega},$ (30)

which leads to the following equation for $f=f(\omega)$

$\displaystyle C_{\tau_{\Pi}}\omega
ff^{\prime}(\omega)+4C_{\tau_{\Pi}}f^{2}+\left(\omega-\frac{16C_{\tau_{\Pi}}}{3}\right)f-\frac{4C_{\eta}}{9}+\frac{16C_{\tau_{\Pi}}}{9}-\frac{2\omega}{3}=0.$
(31)

We can now numerically solve the above equation for a variety of initial
conditions, which we plot in Fig. 1. Irrespective of the choice of initial
conditions, the large $\omega$ behavior converges to a single curve, which is
known as the hydrodynamic attractor. In this concrete sense, the attractor is
said to “forget” about its initial conditions. The approach to the attractor
can be found if we solve the linearized form of (31) for $\delta f$. We find
that

$\displaystyle\delta
f=e^{-\frac{3\omega}{2C_{\tau_{\Pi}}}}\omega^{C_{\eta}/C_{\tau_{\Pi}}}+\mathcal{O}(\omega^{-1}).$
(32)

From here it clear to see that the perturbation decays exponentially to the
attractor.

Note that another often used dimensionless quantity (see e.g.
Florkowski:2017olj ), known as the pressure anisotropy, is defined as

$\displaystyle\mathcal{A}\equiv\frac{\mathcal{P}_{T}-\mathcal{P}_{L}}{\mathcal{E}/3},$
(33)

which is normalized with respect to the equilibrium pressure,
$P=\mathcal{E}/3$. The two dimensionless quantities are related simply via

$\displaystyle\mathcal{A}=18\left(f(\omega)-\frac{2}{3}\right).$ (34)

We can now examine the hydrodynamic gradient expansion as an asymptotic
expansion in $\omega$ of (31), given by

$\displaystyle f$ $\displaystyle=\sum_{n=0}^{\infty}f_{n}\omega^{-n},$ (35)
$\displaystyle=\frac{2}{3}+\frac{4C_{\eta}}{9\omega}+\frac{8C_{\eta}C_{\tau_{\Pi}}}{27\omega^{2}}+\mathcal{O}\left(\omega^{-3}\right).$
(36)

It is straightforward to see that the resulting series has terms growing
factorially Heller:2013fn , indicating that the series is indeed divergent,
see Fig. 2.

Figure 2:  Coefficients in the expansion of $f$ (35) grow like $n!$ Figure
from Heller:2015dha .

It is worthwhile to point out that the choice of variables
$T(\tau)\leftrightarrow f(\omega)$ has nothing to do with the appearance of
the attractor, which is inherent to the system. For the case considered above,
the attractor can alternatively be captured in a study of the phase space
Heller:2020anv , by studying for instance the evolution of
$(\tau_{0}T(\tau),\tau_{0}^{2}\dot{T}(\tau))$ as a function of $\tau$. As food
for thought, the evolution of a set of initial conditions in phase space can
be viewed as dynamical dimensional reduction of the phase space regions. This
might be a helpful way to think about the attractor in more realistic
scenarios, e.g. in the case of relaxed symmetries.

Another method to characterize the attractor is known as the slow-roll
approximation, originally studied in the cosmological setting Liddle:1994dx .
Formally, it amounts to expanding (31) in derivatives, by adding a small
parameter, $\delta$, to the $ff^{\prime}$ term Strickland:2017kux . One then
expands (31) in powers of $\delta$

$\displaystyle f=\sum_{n=0}^{\infty}f_{n}\delta^{n}.$ (37)

This is a useful computational technique, as at a given order, it contains
derivative terms to all orders, although it should be noted that the slow-roll
approximation diverges Denicol:2017lxn .

The Bjorken case is one of the more simple examples due to the assumed
symmetries and low dimensionality. The discussion presented lends itself
naturally to generalizations, which we will review here after a brief aside.

####  2.2.1 Quick introduction to trans-series

Trans-series are an interesting mathematical development with real physics
application. In some sense, they represent a generalization to Taylor series,
taking into account not just the polynomial powers, but also the non-
perturbative structure, such as exponentials and logarithms. A related topic,
resurgence, is the curing of a perturbatively divergent series with non-
perturbative terms Dorigoni:2014hea . Resurgence can be used to deduce
nonperturbative effects via analytic continuation from a perturbative series.
For an introduction to trans-series with a more mathematical flavor, see
transbegin , and for more on resurgence, see for example Aniceto:2018bis .

Here, we will follow the discussion of trans-series relevant to hydrodynamics
Basar:2015ava , noting that the transseries of $N=4$ SYM in Bjorken flow have
been explored in Aniceto:2018uik .

As we see in (32), there is the non-hydrodynamic mode, exponentially supressed
at late times, analogous to quasinormal modes in gravitational theory
Aniceto:2015mto  and to instantons in quantum field theories Schafer:1996wv .
It becomes clear that instead of the naive expansion (35), one should seek a
trans-series ansatz of the form

$\displaystyle
f(\omega)\sim\sum_{n=0}^{\infty}f^{(n)}(\omega)\sigma^{n}\zeta(\omega)^{n},$
(38)

where $\sigma$ is the trans-series expansion parameter, which counts the “non-
perturbative” order of the trans-series and is generally complex and

$\displaystyle\zeta(\omega)\equiv\omega^{\beta}e^{-S\omega},$ (39)

is the so-called “instanton fugacity”.

Plugging this ansatz into (31) leads to equations at each order in $\sigma.$
For $n=0,$ we have the previous case discussed, namely the naive Taylor
series. There is only freedom to fix one parameter, $f^{(1)}_{0}=1$, and all
coefficients are then determined by their recursion relation, which is a
generic feature of trans-series transcostin . The resulting non-hydrodynamic
terms in the trans-series correspond to exponentially damped modes times a
gradient expansion. It is also worthwhile to point out that the overall
constant is known as the Stokes constant, which when fixed resolves the
ambiguity in the choice of contour when doing the inverse Borel transform
Casalderrey-Solana:2017zyh  (note that in the example we consider in Sec. 2.5,
we would have two Stokes coefficients, while they are infinite for
holography).

Do trans-series tell us anything about early times, i.e. when $\omega$ is
small? The exponentially suppressed terms at asymptotically late times can
become large compared to the powers of $\omega$ for early times, with such
terms jockeying for position in size. Trans-asymptotic rearrangement rewrites
the trans-series to better study the early time behavior transmatch . For the
example considered in Basar:2015ava , we should rewrite the expansion (38) as

$\displaystyle f(\omega)\sim\sum_{k=0}^{\infty}F_{k}(\sigma\zeta)\omega^{-k},$
(40)

where
$F_{k}(\sigma\zeta)\equiv\sum_{k=0}^{\infty}f_{k}^{(n)}(\omega)\sigma^{n}\zeta^{n}.$
Inserting the reordered ansatz into (31), one can solve order by order,
finding that for $k>1$ each equation contains just a linear term in $F_{k}.$
The imaginary part of the parameter $\sigma$ is fixed in terms of the Stokes
parameter by the requirement that $f(\omega)$ is real. The real part of
$\sigma$ is arbitrary Aniceto:2013fka , thus parameterizing the different
possible initial conditions at early times.

###  2.3 Gubser flow

A straightforward generalization of Bjorken flow was introduced by Gubser in
Gubser:2010ze , which has been studied at length in various settings. Like
Bjorken flow, Gubser flow is boost-invariant with the additional aspect that
there is azimuthally symmetric radial expansion. In fact, a radially expanding
Bjorken flow can be mapped to Gubser flow via an appropriate Weyl rescaling.
It can be thought of as a comoving flow in conformal theories. The Gubser flow
attractor has been widely studied as it is closer to the set up relevant in
heavy ion collisions Behtash:2017wqg ; Bemfica:2017wps ; Denicol:2018pak ;
Chattopadhyay:2018apf ; Behtash:2019qtk ; Dash:2020zqx .

In terms of the de Sitter coordinates $(\rho,\theta,\phi,\eta)$, the Gubser
flow metric is

$\displaystyle
g_{\mu\nu}=\text{diag}(-1,\cosh^{2}\rho,\cosh^{2}\rho\sin^{2}\theta,1),$ (41)

which relates to the Milne coordinates via

$\displaystyle\sinh\rho$
$\displaystyle=-\frac{1-q^{2}\tau^{2}+q^{2}r^{2}}{2q\tau},$ (42)
$\displaystyle\tan\theta$
$\displaystyle=\frac{2qr}{1+q^{2}\tau^{2}-q^{2}r^{2}},$ (43)

where $q^{-1}$ is the characteristic length scale of the transverse direction
of the system. The time variable is $-\infty<\rho<\infty$. That the system is
expanding is straightforward to see by computing

$\displaystyle\nabla_{\mu}u^{\mu}=2\tanh\rho.$ (44)

One then finds that the conservation of the energy momentum tensor,
$\nabla_{\mu}T^{\mu\nu}=0,$ satisfies a single equation, namely
Marrochio:2013wla ; Behtash:2017wqg

$\displaystyle
0=\frac{1}{T}\partial_{\rho}T+\frac{2}{3}\tanh{\rho}-\frac{1}{3}\pi\tanh\rho,$
(45)

where the normalized shear stress is
$\pi\equiv-\pi^{\eta\eta}/(\varepsilon+P)$ (note that this parameterization is
different from $\phi$ introduced below (26)). As mentioned above, the
evolution of $\pi$ can be specifed in a number of ways, the simplest of which
is arguably MIS (13), which in Gubser coordinates reads Behtash:2017wqg

$\displaystyle
0=\partial_{\rho}\pi+\frac{\pi}{\tau_{\pi}}+\left(\frac{4}{3}\pi^{2}-\frac{4}{3}\frac{\eta}{sT\tau_{\pi}}\right)\tanh{\rho}.$
(46)

Using conformal symmetry, we can write $\tau_{\pi}=c/T$ with $c=5\eta/s$. As
in the Bjorken case, the two equations can be combined into one equation by
solving for $\pi$ in (45). Similar to the discussion in the Sec. 2.2, a
compact way to write the above equations is by introducing the dimensionless
time $\omega=\tanh\rho/T$ and the anisotropy parameter Behtash:2017wqg

$\displaystyle\mathcal{A}=\frac{1}{\tanh{\rho}}\partial_{\rho}\log T,$ (47)

which is related via the conservation equations to the normalized shear stress
as ${\pi}=3\mathcal{A}+2$. We thus have

$\displaystyle
3\omega(\coth^{2}\rho-1-\mathcal{A}(\omega))\mathcal{A}^{\prime}(\omega)+\frac{4}{3}(3\mathcal{A}(\omega)+2)^{2}+\frac{3\mathcal{A}(\omega)+2}{c\omega}-\frac{4}{15}=0$
(48)

As a useful aside, note that DNMR in Gubser flow yield similar evolution
equation for $A(\omega)$ Behtash:2017wqg

$\displaystyle
3\omega(\coth^{2}\rho-1-\mathcal{A}(\omega))\mathcal{A}^{\prime}(\omega)+\frac{4}{3}(3\mathcal{A}(\omega)+2)^{2}+\left(3\mathcal{A}(\omega)+2\right)\left[\frac{1}{c\omega}-\frac{10}{7}\right]-\frac{4}{15}=0.$
(49)

The qualitative differences between the previously discussed Bjorken flow and
Gubser flow were established in a number of studies, although it is worth
pointing out that the gradient expansion of Gubser flow was found to diverge
factorially similarly to the Bjorken case Denicol:2018pak . One of the major
distinctions between Gubser flow and the previously discussed case of Bjorken
flow is the rich range of behaviors Gubser flow experiences. During its
evolution, a system with Gubser dynamics undergoes early time free-streaming,
progressing to intermediate time hydrodynamization and finally at late times
free streaming Dash:2020zqx . The final stage is absent from $1+1$d Bjorken
flow, as a system with Gubser flow expands more rapidly than in a Bjorken flow
and interactions cannot keep up. In examining the local entropy current and
entropy density for the two flows Chattopadhyay:2018apf ;
Chattopadhyay:2018fzy , it was found that Bjorken entropy eventually
saturates, while for Gubser flow, the entropy continues to grow for late
times, marking the ever-growing departure from thermal equilibrium of the
system.

Seeing how Gubser flow has more features closer to real heavy ion collisions
compared to Bjorken flow, such as transverse expansion, one can speculate that
with further relaxed symmetries and the inclusion of transverse dynamics, the
attractor story will more resemble Gubser flow.

###  2.4 Phase transitions

Since hydrodynamic attractors seem to play a role in simpler models of heavy
ion collisions, it would be prudent to study the interplay between the
attractor and other known features of the QGP, such as spontaneous symmetry
breaking. A step in this direction was recently explored in Mitra:2020hbj  by
studying superfluid Bjorken flow, where a conformal fluid in a MIS framework
was coupled to a $U(1)$ scalar field undergoing a phase transition via a
temperature dependent mass. In this case, the action reads

$\displaystyle S=\int d^{4}x\sqrt{-g}\left[p(T)-\frac{1}{2}D_{\mu}\Sigma
D^{\mu}\Sigma^{\star}-V(\Sigma,T)\right],$ (50)

where $\Sigma=\rho e^{i\psi}$ is the complex scalar field, $\rho$ is the
modulus, $\psi$ is the phase, the gauge covariant derivative is
$D_{\mu}\Sigma=(\nabla_{\mu}+iA_{\mu})\Sigma$, the potential is

$\displaystyle V=\frac{1}{2}m^{2}(T)\rho^{2}+\frac{\lambda}{4}\rho^{4},$ (51)

and the metric is the usual Bjorken metric (21). The phase transition is
characterized by parameterizing the temperature-dependent mass via
${m^{2}\sim(T-T_{c})}$. The complete set of equations arise from the energy
momentum tensor conservation, the conservation of the $U(1)$ current, the
scalar equation of motion and the MIS equations. Note that the energy momentum
tensor and the current can be decomposed into a sum of normal and superfluid
components, e.g.

$\displaystyle j^{\mu}=\frac{1}{\sqrt{-g}}\frac{\delta S}{\delta
A_{\mu}}=j_{n}^{\mu}+j_{\psi}^{\mu},$ (52)

where $j^{\mu}_{n}=nu^{\mu}$, $n$ is the number density and
$j^{\mu}_{\psi}=\rho^{2}D^{\mu}\psi.$

As a momentary aside, it is worthwhile to show that the effective fluid action
Brown:1992kc

$\displaystyle S[g_{\mu\nu}]=\int d^{4}x\sqrt{-g}\;p(T),$ (53)

indeed generates the hydrodynamic equations via diffeomorphism invariance. We
define the temperature as

$\displaystyle T\equiv(-\beta^{\mu}g_{\mu\nu}\beta^{\nu})^{-1/2},$ (54)

$\beta^{\mu}=\frac{1}{T}u^{\mu}$ and the four-velocity, $u^{\mu}$, is timelike
normalized, $u^{\mu}u_{\mu}=-1$. If we vary (53), using the thermodynamic
identity involving the entropy, $s=\frac{\partial p}{\partial T},$ and
$\varepsilon+p=sT$, we find

$\displaystyle\delta S$
$\displaystyle=\int\text{d}^{4}x\big{[}\sqrt{-g}\frac{\partial p}{\partial
T}\delta T+\delta\sqrt{-g}p(T)\big{]},$ (55)
$\displaystyle=\int\text{d}^{4}x\sqrt{-g}\big{[}\frac{1}{2}sT^{3}\beta^{\mu}\beta^{\nu}+\frac{1}{2}p(T)g^{\mu\nu}\big{]}\delta
g_{\mu\nu},$ (56)
$\displaystyle=\frac{1}{2}\int\text{d}^{4}x\sqrt{-g}\big{[}sTu^{\mu}u^{\nu}+p(T)g^{\mu\nu}\big{]}\delta
g_{\mu\nu},$ (57)

where we recognize the ideal fluid energy momentum tensor in the last line

$\displaystyle T^{\mu\nu}=(\varepsilon+p)u^{\mu}u^{\nu}+pg^{\mu\nu}.$ (58)

Under diffeomorphisms $x^{\mu}\rightarrow x^{\mu}+\xi^{\mu}$, we have that

$\displaystyle\delta g_{\mu\nu}=\nabla_{\mu}\xi_{\nu}+\nabla_{\nu}\xi_{\mu},$
(59)

which directly leads to the conservation of the stress tensor after partial
integration, $\nabla_{\mu}T^{\mu\nu}=0.$ Since we recover the usual ideal
hydrodynamic equations from the variation of (53), we will refer to (53) as
the fluid action.

Dissipation is added to the energy momentum tensor, $T^{\mu\nu}=T_{\rm
ideal}^{\mu\nu}+\pi^{\mu\nu}$, and the current $j^{\mu}=j^{\mu}_{\rm
ideal}+q^{\mu},$ with relaxation akin to the MIS formulation,

$\displaystyle\left[\tau_{\pi}u^{\alpha}\nabla_{\alpha}+1\right]\pi^{\mu\nu}$
$\displaystyle=-2\eta\sigma^{\mu\nu}-\zeta\Delta^{\mu\nu}\nabla_{\alpha}u^{\alpha},$
(60)
$\displaystyle\left[\tau_{q}u^{\alpha}\nabla_{\alpha}+1\right](q^{\mu}+j^{\mu}_{\psi})$
$\displaystyle=-\kappa\nabla^{\mu}\frac{\mu}{T},$ (61)

where $\tau_{\pi}$ and $\tau_{q}$ are the relaxation times. Furthermore, the
equation for the condensate, $\rho$, and the phase, $\psi$, also have added
noise

$\displaystyle-\nabla_{\mu}\nabla^{\mu}\rho+\frac{\partial
V}{\partial\rho}+\rho D_{\mu}\psi D^{\mu}\psi$ $\displaystyle=\theta_{1},$
(62) $\displaystyle\nabla_{\mu}j^{\mu}_{\psi}$ $\displaystyle=\theta_{2}.$
(63)

Figure 3: Evolution of the system for a variety of initial conditions, which
are shown in green. The inset shows the approach of the system to the usual
hydrodynamic attractor for certain initial conditions, before veering away and
approaching one of two possible fixed points at the end of the system’s
evolution, shown in red and cyan. Figure from Mitra:2020hbj .

We now determine the form of the dissipation by requiring positive entropy
production. The entropy is given by $s=(\varepsilon+P-\mu N)/T$, where $\mu$
is the chemical potential. The divergence of the entropy current reads

$\displaystyle\partial_{\mu}(su^{\mu}-\frac{\mu}{T}(q^{\mu}+j^{\mu}_{\psi}))=$
$\displaystyle-\frac{1}{T}\Pi^{\mu\nu}\nabla_{\mu}u_{\nu}-((q+j_{\psi})\cdot\nabla)\frac{\mu}{T}$
$\displaystyle-\frac{1}{T}\theta_{1}(u\cdot\nabla)\rho-\frac{1}{T}\theta_{2}(u\cdot
D)\psi.$ (64)

Enforcing positive entropy production leads to the following identification of
dissipative structures

$\displaystyle\Pi^{\mu\nu}=-2\eta\sigma^{\mu\nu}-\zeta\Delta^{\mu\nu}\nabla_{\alpha}u^{\alpha},\quad
q^{\mu}+j^{\mu}_{\psi}=-\kappa\nabla^{\mu}\frac{\mu}{T},$ (65)
$\displaystyle\theta_{1}=-\kappa_{1}u^{\alpha}\nabla_{\alpha}\rho,\quad\theta_{2}=-\kappa_{2}u^{\alpha}D_{\alpha}\psi,$
(66)

where $\eta,$ $\zeta,$ $\kappa,$ $\kappa_{1}$ and $\kappa_{2}$ are positive
functions of $T$ and $\mu$.

For a conformal system, $p\sim T^{4}$, the bulk viscosity vanishes, $\zeta=0$,
and the transport coefficients can be parameterized via

$\displaystyle\eta=\frac{4}{3}C_{\eta}T^{3},\quad\tau_{\pi}=\frac{C_{\tau_{\pi}}}{T},\quad\kappa_{1}=C_{\kappa_{1}}T,\quad\kappa_{2}=C_{\kappa_{2}}T^{3}.$
(67)

Defining the pressure anisotropy $\chi:=-3\pi^{\eta}_{\eta}/(4T^{4})$ and
setting $\psi^{\prime}=0$ as well as $\mu=0,$ for simplicity, the equations of
motion are

$\displaystyle\ddot{\rho}+\frac{\dot{\rho}}{\tau}+\lambda\rho^{3}+\sigma(T-T_{c})$
$\displaystyle=-C_{\kappa_{1}}T\dot{\rho},$ (68)
$\displaystyle\frac{\tau\dot{T}}{T}+\frac{1}{3}(1-\chi)-\sigma\frac{\rho^{2}+2\tau\rho\dot{\rho}}{8T^{3}}$
$\displaystyle=\frac{\tau}{4T^{3}}C_{\kappa_{1}}\dot{\rho}^{2},$ (69)
$\displaystyle\tau\dot{\chi}+\frac{4}{3}\left(\chi-\frac{C_{\eta}}{C_{\tau_{\pi}}}\right)+4\chi\tau\frac{\dot{T}}{T}+\frac{\tau
T}{C_{\tau_{\pi}}}\chi$ $\displaystyle=0.$ (70)

The relatively simple system described by (50) and the above equations turns
out to have a rich set of behaviors. Starting in the unbroken phase with
$T>T_{c}$, the system approaches the usual hydrodynamic attractor
($T\propto\tau^{-1/3}$), as can be seen in Fig. 3. During this time, the
condensate, $\rho\equiv|\Sigma|$, value is (close to) zero. Then, some time
later, the system reheats and exponentially quickly develops a non-zero value
for the condensate. This corresponds to a symmetry-breaking fixed point of the
system. The fixed points at $(\pm\rho_{\star},t_{\star}\equiv T/T_{c})$ can be
computed directly from the equations of motion by recognizing that there
should be no dissipation at the fixed point, i.e. $\chi=0$, and
$\dot{\rho}=0,$ which means the equation of motion for $\rho$ leads to

$\displaystyle\frac{\partial
V}{\partial\rho}=0\rightarrow\rho_{\star}=\sqrt{\frac{\sigma}{\lambda}T_{c}(1-t_{\star})}.$
(71)

Plugging this into the equation for the temperature results in

$\displaystyle 1-t_{\star}=\frac{8\lambda
T_{c}^{2}}{3\sigma^{2}}t_{\star}^{3},$ (72)

which in the range $1\geq t_{\star}\geq 0$ always has a real solution,
irrespective of the choice of constants.

Of course, to make contact to phase transitions found in the evolution of the
QGP, like the chiral phase transition, one would need to include the full
$SU(2)_{L}\times SU(2)_{R}\rightarrow SU(2)_{V}$ symmetry breaking pattern,
following the set-up studied in Grossi:2021gqi . As such, this study serves as
a playground for exploring further interplay between hydrodynamic attractors
and other features relevant in the QGP.

###  2.5 Attractor of hybrid models

The QGP has a rich structure of strongly and weakly coupled degrees of
freedom, interacting in a range of energies. To study the interplay between
these two sectors, as a first step it would be worthwhile to study a hybrid
fluid model in a Bjorken background Mitra:2020mei , where the two fluids
represent the hydrodynamic modes of the hard and soft sectors. The central
idea is to see the effect of coupling two viscous fluids which naturally
exhibit attractor behavior to one another. A priori, it is not clear how the
full system will behave or whether a new attractor will emerge. The discussion
here will deviate from Mitra:2020mei , but arrive at the same coupling
equations studied in the Minkowski background of Kurkela:2018dku .

We can study the coupling of two fluids via the fluid action discussed
previously. The action of the first fluid is given by (53), while the second
fluid is

$\displaystyle\tilde{S}[\tilde{g}_{\mu\nu}]=\int\text{d}^{4}x\sqrt{-\tilde{g}}\;\tilde{p}(T_{2}).$
(73)

What is the combined pressure of these two subsystems, (53) and (73), if they
are brought in contact? Taking inspiration from Dalton’s law dalton_law ,
which states that the total pressure is given by the sum of partial pressures
for fixed volume and temperature, i.e.

$\displaystyle P=\sum_{i}P_{i},$ (74)

as well as Agamat’s law bejan_2006 , which states that the total volume is
given by the sum of partial volumes with fixed pressure and temperature, i.e.

$\displaystyle V=\sum_{i}V_{i},$ (75)

we make the ansatz that the total action is

$\displaystyle\mathcal{S}[G_{\mu\nu},T]$
$\displaystyle=S[g_{\mu\nu},T_{1}]+\tilde{S}[\tilde{g}_{\mu\nu},T_{2}],$ (76)
$\displaystyle=\int\text{d}^{4}x\Big{[}\sqrt{-g}P(T_{1})+\sqrt{-\tilde{g}}\tilde{P}(T_{2})\Big{]},$
(77)

where the combined action is given by the total pressure

$\displaystyle\mathcal{S}=\int\text{d}^{4}x\sqrt{-G}\mathcal{P}(T).$ (78)

Note that if the effective metrics, $g_{\mu\nu}$ and $\tilde{g}_{\mu\nu}$, are
the same as the physical background metric, $G_{\mu\nu}$, we then have that
the total pressure as just the sum of the individual pressures, which is just
a statement of Dalton’s law (74). If the pressures are constant and equal,
then the integral of the volume elements can be identified as the partial
volumes (i.e. $\int d^{d}x\sqrt{-\tilde{g}}=\tilde{V}$), whose sum is the
total volume of the system, which is Agamat’s law (75).

As a perhaps useful aside, there seem to be too many metrics. As will be shown
later, the effective metrics, $g_{\mu\nu}$ and $\tilde{g}_{\mu\nu}$, encode
the interaction between the two subsystems. The physical background metric,
$G_{\mu\nu}$, is the metric that could be made dynamical à la PhysRev.94.1468
via

$\displaystyle S=\int d^{d}x\sqrt{-G}(R[G_{\mu\nu}]+\mathcal{P}(T)).$ (79)

Here, we will consider the background metric $G_{\mu\nu}$ to not have any
dynamics, ultimately choosing the Minkowski metric as the background metric,
$G_{\mu\nu}=\eta_{\mu\nu}$.

We will assume that each subsystem is in local thermal equilibrium with a
temperature (Killing) vector
$\tilde{\beta}^{\mu}=\frac{1}{T_{2}}\tilde{u}^{\mu}$ Israel:1979wp , which
means the temperature is given by
$T_{2}=(-\tilde{\beta}^{\mu}\tilde{g}_{\mu\nu}\tilde{\beta}^{\nu})^{-1/2}$.
This leads to the following identification when we are in hydrostatic
equilibrium tolman

$\displaystyle
T\sqrt{-G_{00}}=T_{1}\sqrt{-g_{00}}=T_{2}\sqrt{-\tilde{g}_{00}}=const.$ (80)

We will refer to $T$ as the physical temperature, while $T_{1}$ and $T_{2}$
will be referred to as the effective temperatures of each subsystem.

We can use the thermodynamic identity
$\text{d}\mathcal{P}=\mathcal{S}\text{d}T$ (and also
$\text{d}P=s_{1}\text{d}T_{1}$ and $\text{d}\tilde{P}=s_{2}\text{d}T_{2}$) to
compute the total entropy density. We thus have

$\displaystyle\frac{d\mathcal{P}}{dT}=\frac{\sqrt{-g}}{\sqrt{-G}}\frac{dP}{dT_{1}}\frac{dT_{1}}{dT}+\frac{\sqrt{-\tilde{g}}}{\sqrt{-G}}\frac{d\tilde{P}}{dT_{2}}\frac{dT_{2}}{dT},$
(81)
$\displaystyle\Rightarrow\mathcal{S}\sqrt{-G}=\frac{\sqrt{-g}}{\sqrt{-g_{00}}}s_{1}+\frac{\sqrt{-\tilde{g}}}{\sqrt{-\tilde{g}_{00}}}s_{2},$
(82)

which agrees with the result of the entropy in Kurkela:2018dku . Thus, the
entropy is essentially additive (upto interactions hidden in the effective
metric determinants).

Finally, we point out that the two fluids are coupled only through the
effective metrics, i.e. varying with respect to the effective metrics,
$g_{\mu\nu}$ and $\tilde{g}_{\mu\nu},$ the microscopic detail is given by

$\displaystyle\nabla_{\mu}t^{\mu}_{\phantom{\nu}\nu}=0,\quad\tilde{\nabla}_{\mu}\tilde{t}^{\mu}_{\phantom{\mu}\nu}=0,$
(83)

where the covariant derivative is defined with respect to each subsystem, i.e.

$\displaystyle
0=\nabla_{\mu}t^{\mu\nu}=\partial_{\mu}t^{\mu\nu}+\Gamma^{\mu}_{\phantom{\mu}\mu\alpha}[g_{\rho\sigma}]t^{\alpha\nu}+\Gamma^{\nu}_{\phantom{\mu}\mu\alpha}[g_{\rho\sigma}]t^{\alpha\mu},$
(84)

and similarly for the tilded equations. Thus, one fluid only feels the other
through the effective metric. In other words, the interactions are encoded by
the effective metrics.

We can now compute the full energy-momentum tensor of (76), where we find

$\displaystyle\mathcal{T}^{\mu\nu}=\frac{\sqrt{-g}}{\sqrt{-G}}t^{\alpha\beta}\frac{\partial
g_{\alpha\beta}}{\partial
G_{\mu\nu}}+\frac{\sqrt{-\tilde{g}}}{\sqrt{-G}}\tilde{t}^{\alpha\beta}\frac{\partial\tilde{g}_{\alpha\beta}}{\partial
G_{\mu\nu}}.$ (85)

We see that the change in the effective metrics, $g_{\mu\nu}$ and
$\tilde{g}_{\mu\nu},$ encompasses the interactions between the two fluids.

We introduce interactions by deforming the effective metrics in such a way
that the total energy momentum tensor, $T^{\mu\nu}$, is conserved in the
physical background, $G_{\mu\nu}$, i.e. $\nabla_{\mu}^{(B)}T^{\mu\nu}=0.$
First, we mention that it is straightforward to verify that (84) can be
rewritten as

$\displaystyle\partial_{\mu}(t^{\mu}_{\phantom{\mu}\nu}\sqrt{-g})$
$\displaystyle=\frac{1}{2}\sqrt{-g}t^{\alpha\beta}\partial_{\nu}g_{\alpha\beta},$
(86)

and similarly for the other subsystem. Then, using the previous identity, we
see that the two subsystems are coupled via the lowest order coupling rules

$\displaystyle g_{\mu\nu}$
$\displaystyle=G_{\mu\nu}+\frac{\sqrt{-\tilde{g}}}{\sqrt{-G}}\Bigl{[}\gamma
G_{\alpha\mu}\tilde{t}^{\alpha\beta}G_{\beta\nu}+\gamma^{\prime}\tilde{t}G_{\mu\nu}\Bigr{]},$
$\displaystyle\tilde{g}_{\mu\nu}$
$\displaystyle=G_{\mu\nu}+\frac{\sqrt{-g}}{\sqrt{-G}}\Bigl{[}\gamma
G_{\alpha\mu}t^{\alpha\beta}G_{\beta\nu}+\gamma^{\prime}tG_{\mu\nu}\Bigr{]},$
(87)

where we use the shorthand $t=t^{\mu\nu}G_{\mu\nu}$ and
$\tilde{t}=\tilde{t}_{\mu\nu}G^{\mu\nu}$. Note that the tensor coupling,
$\gamma$, and trace coupling, $\gamma^{\prime}$, have mass dimension $-4$.
Higher order couplings, like
$\mathcal{O}(\gamma^{2},\gamma\gamma^{\prime},(\gamma^{\prime})^{2}$), can be
included in a straightforward manner, see Kurkela:2018dku  for a general
discussion. Using (86) and (2.5), we can read off the full energy momentum
tensor which is conserved in the background metric:

$\displaystyle
T^{\mu}_{\phantom{\mu}\alpha}\sqrt{-G}=\frac{1}{2}\left[(t^{\mu}_{\phantom{\mu}\alpha}+t_{\alpha}^{\phantom{\mu}\mu})\sqrt{-g}+(\tilde{t}^{\mu}_{\phantom{\mu}\alpha}+\tilde{t}_{\alpha}^{\phantom{\mu}\mu})\sqrt{-\tilde{g}}-\delta^{\mu}_{\alpha}\frac{\sqrt{-g}\sqrt{-\tilde{g}}}{\sqrt{-G}}\left(\gamma
t\cdot\tilde{t}+\gamma^{\prime}t\tilde{t}\right)\right],$ (88)

where
$t\cdot\tilde{t}=t^{\mu\nu}G_{\nu\alpha}\tilde{t}^{\alpha\beta}G_{\beta\mu}$.
Note that stress tensors have their indices raised/lowered by their effective
metric, e.g. $t^{\mu}_{\phantom{\mu}\nu}=t^{\mu\alpha}g_{\alpha\nu}.$ The
coupling rules (2.5) are chosen in this way to ensure that the full system has
a conserved energy-momentum tensor in the physical background,
$G_{\mu\nu}=\eta_{\mu\nu}$, namely

$\displaystyle\partial_{\mu}T^{\mu}_{\phantom{\mu}\nu}=0.$ (89)

Figure 4: Anisotropy as a function of proper time. Note that
$\chi=\mathcal{A}/6$, where $\mathcal{A}$ is the anisotropy parameter defined
in (33). The thick lines denote the attractor, while the thin lines denote
solutions with different initial conditions. The more viscous system is in
red, while the blue lines are the less viscous system. Finally the orange line
denotes the emergent system attractor. Figure from Mitra:2020mei .

To solve the coupling equations (2.5), we make the ansatz that the effective
metrics will take the diagonal form, motivated by the diagonal form of the
background metric, $G_{\mu\nu}=\eta_{\mu\nu}$:

$\displaystyle g_{\mu\nu}$
$\displaystyle=\text{diag}(-a^{2},b^{2},b^{2},c^{2}),$ (90)
$\displaystyle\tilde{g}_{\mu\nu}$
$\displaystyle=\text{diag}(-\and^{2},\tilde{b}^{2},\tilde{b}^{2},\tilde{c}^{2}).$
(91)

Likewise, the two stress tensors are assumed to take the form

$\displaystyle t^{\mu\nu}$
$\displaystyle=(\varepsilon+P)u^{\mu}u^{\nu}+Pg^{\mu\nu}+\Pi^{\mu\nu}$ (92)

and similarly for the other, tilded subsystem. $\Pi^{\mu\nu}$ is taken to
evolve using the MIS equation (13) in each sector with each system having
their own shear viscosities, $\eta,\tilde{\eta}$. To distinguish the fluids
and in mind to describe the weak/strong coupling sectors found in the QGP, we
set the tilded subsystem to be more strongly coupled to have
$C_{\tilde{\eta}}=1/4\pi$ to have a smaller shear viscosity than the more
weakly coupled untilded sector

$\displaystyle C_{\eta}$ $\displaystyle=10C_{\tilde{\eta}},\quad
C_{\tau}=5C_{\eta},$ (93) $\displaystyle C_{\tilde{\tau}}$
$\displaystyle=\frac{2-\ln 2}{2\pi}.$ (94)

Thus the untilded weakly coupled sector can be thought of as governed by
kinetic theory, while the tilded strongly coupled sector is governed by a
holography.

It is worth noting that in the decoupling limit,
$\gamma,\gamma^{\prime}\rightarrow 0$, each subsystem contains an attractor.
As can be seen in Fig. 4, with interactions turned on, the combined system
surprisingly has an emergent attractor as well. This emergent attractor is
parameterized by the energy densities of each subsystem, $\gamma\varepsilon$
and $\gamma\tilde{\varepsilon}$. Interestingly, the scenario modelled here is
reminiscent of bottom-up thermalization Baier:2000sb , in that at early times
the strongly coupled subsystem is dominant, while at later times it is the
weakly coupled subsystem which dominates. Intriguingly, though the individual
subsystems were taken to be conformal
($t^{\mu\nu}g_{\mu\nu}=0=\tilde{t}^{\mu\nu}\tilde{g}_{\mu\nu}$), the full
system is in general not, $T^{\mu\nu}G_{\mu\nu}\neq 0,$ becoming conformal at
late times.

The hydrodynamization time for each subsystem can be defined (see
Heller:2016rtz ; Attems:2017zam  for similar definitions) when the pressure is
sufficiently close to first order hydrodynamics, namely

$\displaystyle\frac{|\phi-\phi_{\rm 1st}|}{P}<0.1,\quad\text{for
}\tau>\tau_{\rm hyd},$ (95)

and similarly for the other sector. It was found in a variety of cases that
the hydrodynamization time for the harder sector is longer than for the
strongly interacting sector. One can compare the ratio of the hard to soft
hydrodynamization times via the ratio $R_{\rm hd}=\tau_{\rm
hd}/\tilde{\tau}_{\rm hd}$, where it is clear that depending on the initial
conditions, $R_{\rm hd}$ can vary by an order of magnitude, perhaps indicating
that the hydrodynamization scenario in small-system collisions Chesler:2015bba
; Loizides:2016tew ; Schlichting:2016sqo  may be quite different to their
heavier ion counterpart.

The hybrid coupling described in this section lends itself naturally as a
framework to study more nuanced examples, closer to the physics of the QGP.
For instance, instead of fluids, one can consider the coupling between an
ultraviolet sector described by kinetic theory, while the low energy degrees
of freedom are described holographically. Such a framework is known as
semiholography, a term coined in Faulkner:2010tq  and introduced to heavy ion
collisions in Iancu:2014ava . It would be interesting to see the behavior of
the attractor with the democratic metric coupling in addition to the scalar
coupling between sectors Banerjee:2017ozx ; Ecker:2018ucc ; Mondkar:2021qsf ,
which drives the energy from the hard to the soft sector.

##  3 Holographic attractor

The connection between hydrodynamic and holographic theories is firmly
established, e.g. see Bhattacharyya:2008jc  for a discussion bridging BRSSS
hydrodynamics to a holographic theory. Moreover, boost-invariant holography
had been studied by numerous groups over the years Heller:2011ju ;
Heller:2012je ; Jankowski:2014lna , going back to the linearized analysis of
Janik:2006gp . The holographic setting was also where the Borel resummed
hydrodynamic gradient expansion was first considered Heller:2011ju ;
Heller:2013fn . As such, after finding the hydrodynamic attractor
Heller:2015dha , the stage was set to explore the attractor story in the
strong coupling context. An important early work in this direction was due to
Romatschke Romatschke:2017vte , whose example we will reproduce and expand on
in Sec. 3.1. We will then turn our attention to some of the novel directions
explored in holographic attactors.

###  3.1 An illustrative example

Here we will follow Romatschke:2017vte . In the holographic context, one aims
to solve a classical gravity problem which is dual to a strongly coupled
quantum field theory. We will want this dual theory to be comparable to the
hydrodynamic theory. The dual theory is encoded in the boundary conditions of
the gravity equations, which we will outline shortly. We begin with the vacuum
Einstein equations

$\displaystyle R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R+\Lambda g_{\mu\nu}=0,$ (96)

where $R_{\mu\nu}$ is the Ricci tensor and $R\equiv R_{\mu\nu}g^{\mu\nu}$ is
the Ricci scalar. The cosmological constant in AdS is given by
$\Lambda=-{(d-1)(d-2)}/{2L^{2}}=-6$, since we work in $d=5$ dimensions and
with the AdS radius, $L$, set to unity.

The metric in Eddington-Finkelstein coordinates is given by

$\displaystyle ds^{2}=2drd\tau-
Ad\tau^{2}+\Sigma^{2}\left(e^{B}dx_{\perp}^{2}+e^{-2B}d\xi^{2}\right),$ (97)

where $\xi$ is the spacetime rapidity. We take the functions $A$, $\Sigma$ and
$B$ to be functions of the proper time, $\tau$, and the holographic radius,
$r$. With this metric ansatz, the Einstein equations (96) can be helpfully
written as a set of nested ordinary differential equations Chesler:2008hg

$\displaystyle 0$
$\displaystyle=\Sigma\dot{\Sigma}^{\prime}+2\Sigma^{\prime}\dot{\Sigma}-2\Sigma^{2},$
(98) $\displaystyle 0$
$\displaystyle=2\Sigma\dot{B}^{\prime}+3(\Sigma^{\prime}\dot{B}+B^{\prime}\dot{\Sigma}),$
(99) $\displaystyle 0$
$\displaystyle=A^{\prime\prime}+3B^{\prime}\dot{B}-12\Sigma^{\prime}\dot{\Sigma}\Sigma^{-2}+4,$
(100) $\displaystyle 0$
$\displaystyle=2\ddot{\Sigma}+\dot{B}^{2}\Sigma-A^{\prime}\dot{\Sigma}^{\prime},$
(101) $\displaystyle 0$ $\displaystyle=2\Sigma^{\prime\prime}+B^{\prime
2}\Sigma,$ (102)

where the prime denotes the radial derivative,
$h^{\prime}\equiv\partial_{r}h$, and the dot denotes the radial derivative
along the null geodesic

$\displaystyle\dot{h}\equiv\partial_{\tau}h+\frac{1}{2}A\partial_{r}h.$ (103)

Note that (101) and (102) are constraint equations, while the other equations
are to be solved numerically.

Near the boundary at $r=\infty,$ the asymptotic expansion of the metric
functions is given by

$\displaystyle A(\tau,r)$
$\displaystyle=r^{2}\sum_{n=0}^{\infty}a_{n}(\tau)r^{-n},$ (104)
$\displaystyle B(\tau,r)$
$\displaystyle=\sum_{n=0}^{\infty}a_{n}(\tau)r^{-n},$ (105)
$\displaystyle\Sigma(\tau,r)$
$\displaystyle=r\sum_{n=0}^{\infty}a_{n}(\tau)r^{-n}.$ (106)

We want the metric at the boundary to take the form of the Bjorken metric,
(21), in order for the dual theory to exhibit boost invariance. This implies
that the boundary conditions at $r\rightarrow\infty$ are

$\displaystyle A=r^{2},\quad
B=-\frac{2}{3}\ln{\frac{1+r\tau}{r}}\quad\text{and}\quad\Sigma^{3}=r^{2}(1+r\tau).$
(107)

This determines the leading coefficients in the expansion of the metric
functions, $a_{0}=1$ and $s_{0}=1.$ Furthermore, it is easy to see that the
metric (97) remains invariant under the diffeomorphism, $r\rightarrow
r+f(\tau).$ This residual gauge freedom can be fixed by setting $a_{1}=0.$

The boundary stress tensor

$\displaystyle T_{\mu\nu}=\text{diag}(\varepsilon,p_{T},p_{T},p_{L})$ (108)

is related to the metric coefficient $a_{4}$ via holographic renormalization
Balasubramanian:1999re ; deHaro:2000vlm  by

$\displaystyle\varepsilon$ $\displaystyle=-\frac{3\kappa}{4}a_{4},$ (109)
$\displaystyle p_{L}$
$\displaystyle=-\varepsilon-\tau\partial_{\tau}\varepsilon,$ (110)
$\displaystyle p_{T}$
$\displaystyle=\varepsilon+\frac{1}{2}\tau\partial_{\tau}\varepsilon.$ (111)

The Einstein equations are then solved using numerical techniques due to
Chesler and Yaffe (namely, the spectral method) Chesler:2008hg ;
Chesler:2009cy ; Chesler:2013lia . There is also publicly available code in
Wu:2011yd ; paul-attractor .

Various initial conditions at $\tau=\tau_{0}$ were generated in
Romatschke:2017vte  by taking $\Sigma$ to be

$\displaystyle\Sigma^{3}=r^{2}(1+r\tau)+\frac{r^{2}s_{2}+rs_{3}+s_{4}}{r^{4}+\bar{c}},$
(112)

where

$\displaystyle s_{2}$
$\displaystyle=\frac{1}{20}\left(4\varepsilon+3\tau\partial_{\tau}\varepsilon\right),$
(113) $\displaystyle s_{3}$
$\displaystyle=\frac{1}{40}\left(13\partial_{\tau}\varepsilon+5\tau\partial^{2}_{\tau}\varepsilon\right),$
(114) $\displaystyle s_{4}$
$\displaystyle=\frac{1}{560\tau}\left(9\partial_{\tau}\varepsilon+135\tau\partial_{\tau}^{2}\varepsilon+35\tau^{2}\partial_{\tau}^{3}\varepsilon\right).$
(115)

Then choosing the profile $\varepsilon\sim\tau^{\alpha}$, one can generate
initial data by changing $\alpha$ and $\bar{c}.$

Sorting through the initial data, there is a clear late-time attractor as can
be seen in Fig. 5. However, at earlier times the story is not so clear due to
the well-known oscillatory nature of $N=4$ SYM, which arises primarily from
the existence of non-hydrodynamic modes at early times Romatschke:2017vte ;
Spalinski:2017mel .

To identify the holographic attractor, it was posited in Romatschke:2017vte
that the relevant condition should be that the attractor experiences little
oscillatory behavior and whose initial value tends to

$\displaystyle\lim_{\tau\rightarrow
0}\frac{\partial\ln\varepsilon}{\partial\ln\tau}\rightarrow-1.$ (116)

The attractor is then visible in the rightmost panel of Fig. 5 as a solid
black line, while the various initial conditions produce a variety of
oscillatory behavior. Such an ambiguity is not found in relativistic
hydrodynamics and kinetic theory, where the approach to the attractor tends to
be monotonic from a given set of initial conditions.

Figure 5: The attractor for rBRSSS hydrodynamics (left panel), kinetic theory
in the relaxed time approximation (middle panel) and the holographic theory
from Sec. 3.1 (right panel). The gray lines denote solutions with different
initial conditions, while the solid lines are the numerical attractor. The
late time behavior should approach the hydrodynamics, which are denoted by the
dashed lines. Note that the agreement between the attractor and hydrodynamics
at late times occur at different timescales, $\omega=\tau T$, for the
different theories. Figure from Romatschke:2017vte .

In contrast to the Romatschke approach, another study was undertaken in
Kurkela:2019set , where the initial conditions of the anisotropic part of the
metric was modelled instead in two ways (note: $\rho=1/r$)

$\displaystyle
B(\tau_{0},\rho)+\frac{2}{3}\log{(\tau_{0}+\rho)}=\begin{cases}e^{-40\rho
T}32\rho^{5}T^{5}&\text{(UV)},\\\ 32\rho^{5}T^{5}&\text{(IR)}.\end{cases}$
(117)

This particular choice of initial conditions was chosen such that the initial
ratio of the longitudinal pressure to the energy density
$p_{L}(\tau_{0})/\varepsilon(\tau_{0})$ and its first derivative are equal.
The meaningful difference between the two families of initial conditions is in
their extent in the holographic direction, with the UV solution localized
almost entirely on the boundary. Both sets of initial conditions display the
expected oscillatory behavior due to the expected decay of quasinormal modes.
Intriguingly, the more localized UV initial condition seems to approach the
hydrodynamic behavior sooner than the corresponding IR family of initial
conditions. For either family of initial conditions, the approach to the
attractor takes longer than the corresponding set-up in kinetic theory or IS
hydrodynamics, irrespective of initialization time. Thus, results of
Kurkela:2019set  seem to indicate that the attractor only exists in a
meaningful way at late times for holographic set-ups with the early time
dynamics dominated by transients, which is further supported by related work
of higher-order fluid dynamics Chattopadhyay:2018pwe .

####  3.1.1 Borel summation

Whether the condition (116) is the “best” or unambigous identifier for a
holographic attractor is up for debate. A perhaps more rigorous approach is
due to Spalinski:2017mel ; Spalinski:2018mqg . Using the same metric as in
(97), the hydrodynamic gradient expansion of the energy density in terms of
the proper time

$\displaystyle\varepsilon=\frac{3N_{c}^{2}\pi^{2}}{8}\tau^{-4/3}\sum_{i=2}^{\infty}\varepsilon_{i}\tau^{-2(i-2)/3}$
(118)

was computed up to $240^{th}$ order Heller:2013fn . The longitudinal and
transverse pressure in the Bjorken case is given by (110) and (111), so that
one can easily compute the anisotropy parameter (33), which has an expansion
in terms of $\omega=T\tau$:

$\displaystyle A(\omega)=\sum_{k=1}^{\infty}a_{k}\omega^{-k},$ (119)

Thus, the expansion of the energy density in terms of proper time (118) can be
re-expressed as an expansion in the anisotropy parameter in terms of the
dimensionless time, $\varepsilon(\tau)\leftrightarrow A(\omega)$. One then
computes the Borel transform of the original power series, which is given by

$\displaystyle\mathcal{B}A(x)=\sum_{k=1}^{\infty}\frac{a_{k}}{k!}x^{k}$ (120)

The intent of the Borel transform is to provide a method of handling divergent
series. The gradient expansion (119) is a divergent expansion with larger
order terms diverging factorially Heller:2015dha . One can then perform the
Borel resummation via the inverse Borel transform

$\displaystyle A(\omega)=\omega\int_{C}dxe^{-\omega x}\mathcal{B}A(x),$ (121)

where the contour, $C$, connects $x=0$ to $x=\infty$ in the complex plane.
Practically, since only 240 terms are known due to the high technical cost of
computing the coefficients, the analytic continuation is done using Padé
approximants. The Padé approximants are, in some respects, a generalization of
the Taylor series, with the improvement that singularities are explicitly
included in the approximant. This is done by considering a ratio of
polynomials, i.e. we can approximate a function via

$\displaystyle
f(x)\approx\frac{\sum_{i}^{N}a_{i}x^{i}}{1+\sum_{j}^{M}b_{j}x^{j}},$ (122)

where the denominator can better capture the singularity structure of the
function that is being approximated, compared to Taylor series, which does not
contain such terms. For more on Padé approximants, see numrecipe  and
Ellis:1995jv  for application to QCD sum rules.

Does this method provide a good approximation to the attractor? In the case
the attractor is known, such as in the previously mentioned HSJW model (20),
the known attractor can be compared to the result of the integrating the
series truncated at $240^{th}$ order via (121), where there was quantitative
agreement down to about $\omega\approx 0.3,$ indicating the method provides
meaningful results. In the case of holography, where the exact attractor is
not known, it was found in Spalinski:2017mel  that there is good agreement
with Romatschke:2017vte  up to $\omega\sim 0.4,$ whereas the qualitative
behavior is different for smaller $\omega$.

###  3.2 Violation of energy conditions

As the main workhorse for the study of holographic attractors and the go-to
holographic toy model for heavy ion collisions, it is important to understand
the holographic Bjorken model in Sec. 3.1. For this reason, it is interesting
to note the holographic Bjorken flow can evolve to violate the dominant and
weak energy conditions Rougemont:2021qyk , even with initial conditions that
satisfy these bounds. For conformal Bjorken flow, the weak energy condition
was given in the first holographic study of Bjorken flow Janik:2005zt

$\displaystyle T^{\mu\nu}t^{\mu}t^{\nu}\geq 0,$ (123)

where $t^{\mu}$ is any time-like vector, which implies for Bjorken flow

$\displaystyle\varepsilon(\tau)\geq 0,\quad\varepsilon^{\prime}(\tau)\leq
0,\quad\text{and}\quad\tau\varepsilon^{\prime}(\tau)\geq-4\varepsilon(\tau).$
(124)

Note that since we are considering a conformal system, the strong energy
condition

$\displaystyle T_{\mu\nu}t^{\mu}t^{\nu}\geq\frac{T^{\mu}_{\mu}}{2}=0,$ (125)

is equivalent to the weak energy condition. The dominant energy condition
implies that

$\displaystyle\varepsilon(\tau)>0,\quad\text{and}\quad-1<\frac{p_{T}-p_{L}}{\varepsilon}<2.$
(126)

It was found that the violations occur at times around $\tau\approx 0.5\text{
fm/c}$, which is within the prehydrodynamic regime in the QGP Kurkela:2018wud
. This can have implications for the validity of kinetic theory, the typical
effective description used in this regime. Kinetic theory satisfies the energy
conditions: the weak energy condition implies that the energy density is non-
negative and the dominant energy condition is the requirement that matter does
not exceed the speed of light wald_2009 . If these violations were shown to
exist in the energy mometum tensor during nuclear collisions, this could
indicate that this transient regime would lie outside of the domain of
validity of kinetic theory descriptions Rougemont:2021qyk , indicating that
strongly-coupled hydrodynamization as predicted by holography would be the
preferred description.

###  3.3 Gauss-Bonnet attractor

Attractors have been studied in Gauss-Bonnet (GB) holography Casalderrey-
Solana:2017zyh ; Gushterov:2018erj ; Casalderrey-Solana:2019npu . The key
advantage of studying GB holographic systems is that one has access to
intermediate coupling (in addition to infinitely strong and weak coupling)
regimes. The only parameter is $\lambda_{GB}$, which controls the size of the
higher derivative corrections. The dual field theory is unknown, but finite
and negative values of $\lambda_{GB}$ are at least qualitatively similar to
what one would expect from finite coupling corrections to $N=4$ SYM. The
negative values of $\lambda_{GB}$ do lead to violations of causality in the
ultraviolet Camanho:2014apa , but the interest of these works is in the
direction of the infrared dynamics and the system’s approach to hydrodynamics.

Here, we will briefly outline some details of the construction and the lessons
learnt. The GB action is

$\displaystyle
S=\frac{1}{2\kappa_{5}^{2}}\int\text{d}^{5}x\sqrt{-g}\left[R+\frac{12}{L^{2}}+\frac{\lambda_{GB}L^{2}}{2}(R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}-4R_{\mu\nu}R^{\mu\nu}+R^{2})\right],$
(127)

where $\kappa_{5}$ is related to Newton’s constant in 5 dimensions and $L$ is
the AdS radius for the $\lambda_{GB}=0$ theory. Note that the $\lambda_{GB}=0$
case corresponds to the example in Sec. 3.1 with the terms proportional to
$\lambda_{GB}$ representing higher derivatives of the metric. The shear
viscosity of this theory is given by

$\displaystyle\frac{\eta}{s}=\frac{1-4\lambda_{GB}}{4\pi}.$ (128)

Using a similar metric ansatz of the Eddington-Finkelstein form (97), one then
solves the Einstein equations as outlined in Sec. 3.1.

For the values of $\lambda_{GB}$ studied in Casalderrey-Solana:2017zyh , the
attractor is estimated using Borel summation as outlined in Sec. 3.1.1. It was
found that the relaxation time of the non-hydrodynamic modes occurs when the
gradients of the system are still large, i.e. for anisotropy parameter $>1$,
with first order viscous hydrodynamics comparable to the estimated attractor.
Interestingly, it is found that by varying $\lambda_{GB}$ there is a
qualitative interpolation in the singularities found in the complex plane
between the holographic infinite coupling limit and the weakly coupled limit
in kinetic theory, which can be inferred from Fig. 6.

Figure 6: Positions of poles of the Padé approximant in the complex $\chi$
plane (which we call $x$ in Sec. 3.1.1) for different values of the coupling
$\lambda_{GB}$, compared to the RTA kinetic theory result. Note that
$\lambda_{GB}=0$ represents the typical pole structure seen in holography.
Figure from Casalderrey-Solana:2017zyh .

###  3.4 $D$-dimensional Bjorken flow

The attractor in arbitrary $D=n+1$ dimensions was studied in Casalderrey-
Solana:2018uag . For a review of large $D$ gravity, see Emparan:2020inr . The
essential technical difference is that the metric (97) has the same form and
dependence, but there are more transverse directions, i.e.

$\displaystyle ds^{2}=2drd\tau-
Ad\tau^{2}+\Sigma^{2}\left(e^{B}dx_{\perp}^{2}+e^{-(n-2)B}d\xi^{2}\right),$
(129)

which is chosen such that the determinant of the term in brackets is equal to
one. Operationally, the $D-$dimensional Einstein equations are rewritten in
the same way as in (98)-(102). The equations are then expanded in a large
$n=D-1$ limit up to $n^{-3}$ order. The key ingredient to compare holographic
theories in different dimensions is the requirement that the late time
behavior is the same, which is enforced by introducing the late time scale,
$\Lambda,$ which arises from considering the late time behavior of the ideal
Bjorken expansion

$\displaystyle\varepsilon_{\rm
ideal}=\frac{\Lambda^{n}}{(\Lambda\tau)^{n/(n-1)}}.$ (130)

Interestingly in the case studied in Casalderrey-Solana:2018uag , at the time
scale $\tau\sim 1/\Lambda$ the decoupling of non-hydrodynamic modes from the
hydrodynamic modes occurs in an average sense: these modes are rapidly
oscillating around the hydrodynamized evolution, averaging to zero over the
typical expansion time scale. Coupled with the fact that the non-hydrodynamic
excitations have amplitudes which are not parametrically suppressed, it seems
like the typical approach to the attractor would be different to the usual
$D=5$ case.

##  4 Attractor in kinetic theory

In heavy ion collisions, the conceptual picture is that kinetic theory tracks
the evolution of (quasi-)particles. The applicability of the theory is early
times post collision up to the time when the build-up of the soft gluon bath
is significant Busza:2018rrf . This corresponds to roughly the
hydrodynamization time, $\tau\sim 1\text{ fm/c}$. The bridging between regimes
where the dominant effective description is kinetic theory to hydrodynamics
has been developed in KøMPøST Kurkela:2018wud ; Kurkela:2018vqr .

As a workable theory when the heavy ion system is quite far-from-equilibrium,
it is interesting to study which role attractors play in a kinetic
description. To that end, we present the reader with an example kinetic
computation in Sec. 4.1, before describing further progress in the following
subsections.

###  4.1 Kinetic theory primer

Here we provide a short recap of key concepts regarding kinetic theory,
following the discussion in Gorbunov:2011zzc . Kinetic theory concerns itself
with the microscopic detail of particles directly. The fundamental object of
interest in kinetic theory, keeping track of the particles’ positions and
momenta, is the one particle distribution function, $f=f(x^{\mu},p_{\mu})$,
where $x^{\mu}$ are the spacetime coordinates and $p_{\mu}$ is the four-
momentum. Note that the distribution function should not be confused with the
dimensionless quantity (30), used to characterize the hydrodynamic attractor.
The evolution of the distribution function in phase space follows Liouville’s
theorem landau , namely

$\displaystyle
0=\frac{df}{d\lambda}=\frac{dx^{\mu}}{d\lambda}\partial_{\mu}f+\frac{dp_{\mu}}{d\lambda}\frac{\partial
f}{\partial p_{\mu}}$ (131)

where $\lambda$ is an affine parameter. This equation is commonly known as the
collisionless Boltzmann equation. It is common to consider the Boltzmann
equation on-shell, e.g. $p_{\mu}p^{\mu}=-m^{2}$ in our sign convention. This
means that the one of the components of the four-momentum is determined in
terms of the other momenta with a common choice being $p_{0}=p_{0}(p_{i}).$
For a massive gas in a Minkowski background, this is simply
$p_{0}=\sqrt{p_{i}\delta^{ij}p_{j}+m^{2}}.$

Macroscopic quantities can be computed straightforwardly by considering the
moments of the distrbution functions, i.e. integrating the distribution
function with increasing powers of the four-momentum over momentum. For
instance, the particle number current and energy momentum tensor are given by
Andreasson:2011ng ; Gorbunov:2011zzc

$\displaystyle n^{\mu}(x^{\mu})$
$\displaystyle=\int\frac{d^{3}p}{(2\pi)^{3}\sqrt{-g}}\frac{p^{\mu}}{p^{0}}f(x^{\mu},p_{\mu}),$
(132) $\displaystyle T^{\mu\nu}(x^{\mu})$
$\displaystyle=\int\frac{d^{3}p}{(2\pi)^{3}\sqrt{-g}}\frac{p^{\mu}p^{\nu}}{p^{0}}f(x^{\mu},p_{\mu}),$
(133)

respectively.

Noting that the four-momentum is $p^{\mu}=\frac{dx^{\mu}}{d\lambda}$ and the
four-force is $F_{\mu}=\frac{dp_{\mu}}{d\lambda}$, we find that the
distribution function evolves according to the Boltzmann equation in a curved
spacetime in the absence of external forces

$\displaystyle
p^{\mu}\partial_{\mu}f+p_{\mu}p^{\beta}\Gamma^{\mu}_{\phantom{\mu}\alpha\beta}\frac{\partial
f}{\partial p_{\alpha}}=C[f],$ (134)

where $C[f]$ is the collision kernel, which encompasses all order scattering
processes. In practice, it is more expedient to work with a simpler, truncated
collision kernel. For analytic study, a helpful choice of the collision kernel
is the relaxation time approximation (RTA) (also known as the Bhatnagar-Gross-
Krook (BGK) approximation Bhatnagar:1954zz ), where the collision term is
given by

$\displaystyle C[f]=\frac{p^{\mu}u_{\mu}}{\tau_{R}}(f-f_{eq}),$ (135)

where $\tau_{R}$ is the relaxation time, which in principle can depend on the
proper time, $\tau_{R}=\tau_{R}(\tau)$. Here, we will take the equilibrium
distribution function to be given by

$\displaystyle f_{eq}=e^{p^{\mu}u_{\mu}/T}$ (136)

where $T$ is the temperature and will be a function of proper time, $\tau$.
The interpretation of the RTA collision term is that the system relaxes back
to its equilibrium value, which can be easily seen in the following example.
Assume that we are in flat Minkowski spacetime with only time dependence in
the distribution function. Taking $f_{\rm eq}$ to be a constant and shifting
$f-f_{\rm eq}\rightarrow f$ to write the Boltzmann equation (134) in the RTA
(135) as

$\displaystyle\partial_{0}f=-\frac{1}{\tau_{R}}f,$ (137)

which has solution after shifting back as

$\displaystyle f=f_{\rm eq}+e^{-t/\tau_{R}}.$ (138)

Thus, we see that the distribution function exponentially rapidly approaches
its equilibrium value in a timescale, $\tau_{R}.$

An important ingredient in the RTA is what is known as the matching condition.
Since the energy momentum tensor (133) is conserved,
$\nabla_{\mu}T^{\mu\nu}=0,$ the first moment of the Boltzmann equation must
vanish. This sets a condition on the collision kernel, namely that the moment
of the collision kernel should vanish, leading to the dynamical matching
conditions. In other words, the energy density, $\varepsilon$, should match
the equilibrium energy density, $\varepsilon_{\rm eq}$, i.e.

$\displaystyle T^{00}$ $\displaystyle=T^{00}_{\rm eq},$ (139)
$\displaystyle\Rightarrow\int\frac{d^{3}p}{(2\pi)^{3}}\frac{(p^{0})^{2}}{p^{0}}(f-f_{\rm
eq})=0.$ (140)

####  4.1.1 Attractor via the method of moments

A straightforward way of seeing the attractor in kinetic theory is using the
method of moments. In the context of attractors, the tower of moments was
studied in the RTA Kurkela:2019set ; Strickland:2019hff  as well as with more
realistic collision kernels Almaalol:2020rnu , following prior work in
studying the divergence of the Chapman-Enskog expansion Denicol:2016bjh .
Moreover, for an application of studying the fixed points, e.g. free streaming
and hydrodynamic, found in kinetic theory via the method of moments, see
Blaizot:2017ucy ; Blaizot:2019scw ; Blaizot:2021cdv . Here we outline the
basics of the construction.

Instead of working with the distribution function directly, it is helpful to
consider the generalized moments, $\rho$, of the distribution function
Denicol:2016bjh

$\displaystyle\rho_{n,l}=\int\frac{d^{3}k}{(2\pi)^{3}\tau}(k^{0})^{n}\left(\frac{k_{\eta}}{k^{0}\tau}\right)^{2l}f(\tau,k_{0},k_{\eta}).$
(141)

The moments contain all of the microscopic information of the distribution
function. For instance, it is clear by inspection that the lower moments have
a straightforward physical interpretation, namely the energy density is given
by $\varepsilon=\rho_{1,0}$ and the longitudinal momentum is given by
$p_{L}=\rho_{1,1}$.

The evolution equation for the moments is found by integrating the Boltzmann
equation. One sees that

$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}\tau}(k^{0})^{n}\left(\frac{k_{\eta}}{k^{0}\tau}\right)^{2l}\partial_{\tau}f$
$\displaystyle=\partial_{\tau}\rho_{n,l}-\int\frac{d^{3}k}{(2\pi)^{3}}k_{\eta}^{2l}\partial_{\tau}\left[\tau^{-1-2l}(k^{0})^{n-2l}\right]f,$
(142)
$\displaystyle=\partial_{\tau}\rho_{n,l}+\frac{1+2l}{\tau}\rho_{n,l}+\frac{n-2l}{\tau}\rho_{n,l+1},$
(143)

where we integrated by parts and used that in the Bjorken massless case
$k_{0}=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{\eta}^{2}/\tau^{2}}$. Integrating over the
equilibrium distribution function leads to

$\displaystyle\rho^{\rm
eq}_{n,l}=\int\frac{d^{3}k}{(2\pi)^{3}}(k^{0})^{n}\left(\frac{k_{\eta}}{\tau
k^{0}}\right)^{2m}f_{\rm eq}.$ (144)

Putting together the pieces, we find that the general equation for moments is
given by

$\displaystyle\partial_{\tau}\rho_{n,l}+\frac{1+2l}{\tau}\rho_{n,l}+\frac{n-2l}{\tau}\rho_{n,l+1}=-\frac{\rho_{n,l}-\rho^{\rm
eq}_{nm}}{\tau_{R}}.$ (145)

In the case considered by Kurkela:2019set , namely the set of lowest lying
equations $(n,l)=(1,0)$ and $(n,l)=(1,1)$, one finds the following set of
equations

$\displaystyle 0$
$\displaystyle=\partial_{\tau}\varepsilon+\frac{1}{\tau}(\varepsilon+p_{1}),$
(146) $\displaystyle 0$
$\displaystyle=\partial_{\tau}p_{1}+\frac{1}{\tau}(3p_{1}-p_{2})+\frac{p_{1}-\varepsilon/3}{\tau_{R}}.$
(147)

As a reminder, the relaxation time is in principle a function of proper time.
As above, we will consider the conformal case, where
$\tau_{R}^{-1}\propto\varepsilon^{1/4}.$ We will use the rescaled time
$t\equiv\tau/\tau_{R}$ (not to be confused with the Minkowski time, $t$). A
helpful identity to be used below is that

$\displaystyle\frac{1}{\tau_{R}}\frac{dt}{d\tau}=\frac{1}{\varepsilon^{1/4}}\frac{d(\tau\varepsilon^{1/4})}{d\tau}=1+\frac{\tau}{4}\partial_{\tau}\varepsilon=\frac{1}{4}(3-x),$
(148)

where we introduced $x\equiv p_{1}/\varepsilon$. Finally, calling $y(t)\equiv
p_{2}/\varepsilon$, we see that (146) and (147) can be combined in one
equation as Kurkela:2019set

$\displaystyle\frac{3}{4}\left(3-x\right)tx^{\prime}(t)=3x^{2}-6x+3y+t(1-3x).$
(149)

In analyzing (149), a useful starting point is to consider is at $t=0.$
Requiring the derivative of $x$ is finite at $t=0,$ we see that we have two
cases

$\displaystyle x(0)=1\pm\sqrt{1-y(0)}.$ (150)

At late times, it is clear that (149) becomes

$\displaystyle\frac{3}{4}(3-x)x^{\prime}(t)=1-3x,$ (151)

from which it follows that the generic late-time behavior is
$\lim_{t\rightarrow\infty}x(t)\rightarrow\frac{1}{3}.$

Depending on the initialization time, the approach to the attractor, which we
call $x_{A},$ has two qualitatively different regimes, as can be seen in Fig.
7. When the initial time is chosen to be $t_{0}\gg 1$, the attractor is
approached exponentially quickly, while for $t_{0}\ll 1$ the approach to the
attractor follows approximately a power law decay, $x-x_{A}\sim t^{-8/3}$.
Similar behavior, namely a transition from power law to exponential behavior,
is also seen in the case of when an approximate analytic attractor is known
Jaiswal:2019cju .

It is worth bearing in mind that the number of moments is infinite. It is
beneficial to understand how the attractor might appear at higher moments in
the distribution function. Work in the Bjorken flow in the RTA
Strickland:2018ayk ; Strickland:2019hff  indicates that the higher moments all
approach the attractor at late times, although higher moments with $l=0$ in
(141) approach the exact attractor solution much slower than the other modes.

Figure 7: The RTA attractor. The ratio of the longitudinal pressure to energy
density is plotted for various initializations of (149). Figure from
Kurkela:2019set .

###  4.2 Exact solution in the RTA

Another method to study the attractor in kinetic theory is the use of the
exact solution of the Boltzmann RTA, following e.g. Strickland:2018ayk . Note
that an analytic solution in terms of confluent hypergeometric functions was
found in the RTA with massless particles with a relaxation time related to the
proper time via $\tau_{R}\propto\tau^{1-\Delta}$, where $\Delta$ is a constant
Blaizot:2020gql ; Blaizot:2021cdv . As in the hydrodynamic case, we will limit
our initial consideration to $0+1$ dimensional Bjorken expansion. The RTA
attractor has been further examined in Strickland:2019hff .

We begin by taking the flow velocity to be

$\displaystyle z^{\mu}=(z\tau^{-1},0,0,t\tau^{-1}),$ (152)

which is clearly spacelike
$z^{\mu}\eta_{\mu\nu}z^{\nu}=({t^{2}-z^{2}})/{\tau^{2}}=1$. Due to our
assumption of boost invariance, the distribution function will depend only on
the proper time, $\tau$, the transverse momenta, $\vec{p}_{T}$, and the boost-
invariant variable PhysRevD.30.2371 ; Bialas:1987en

$\displaystyle\omega=tp_{L}-zE,$ (153)

where $E=p^{0}.$ One can then define the quantity

$\displaystyle v(\tau,\omega,\vec{p}_{T})=Et-
p_{L}z=\sqrt{\omega^{2}+(m^{2}+\vec{p}_{T}^{2})\tau^{2}},$ (154)

We will consider massless particles, setting $m=0$ in what follows. The above
relations lead to a straightforward identification of the energy and
longitudinal momentum of a particle as

$\displaystyle E\tau^{2}$ $\displaystyle=vt+\omega z,$ (155) $\displaystyle
p_{L}\tau^{2}$ $\displaystyle=\omega t+vz.$ (156)

Note that in this case, the integration measure reads

$\displaystyle dP$
$\displaystyle=\frac{d^{4}p}{(2\pi)^{4}}2\pi\delta(p^{2}-m^{2})2\theta(E),$
(157) $\displaystyle=\frac{dp_{L}d^{2}p_{T}}{(2\pi)^{3}E}=\frac{d\omega
d^{2}p_{T}}{(2\pi)^{3}v}$ (158)

With the simplifying assumptions and parameterizations outlined above, the RTA
Boltzmann equation takes the simple form

$\displaystyle\partial_{\tau}f=\frac{f_{(eq)}-f}{\tau_{R}(\tau)},$ (159)

where we note that the relaxation time $\tau_{R}$ is a function of proper
time. The RTA Boltzmann equation has a formal solution BAYM198418 ;
Florkowski:2013lza  given by

$\displaystyle
f(\tau,\omega,p_{T})=D(\tau,\tau_{0})f_{0}(\omega,p_{T})+\int^{\tau}_{\tau_{0}}\frac{\text{d}\tau^{\prime}}{\tau_{\rm
eq}(\tau^{\prime})}D(\tau,\tau^{\prime})f_{\rm
eq}(\tau^{\prime},\omega,p_{T}),$ (160)

where the damping function is given by

$\displaystyle
D(a,b)=\exp\left[-\int^{a}_{b}\frac{\text{d}\tau^{\prime}}{\tau_{\rm
eq}(\tau^{\prime})}\right]$ (161)

and $f_{0}$ is an arbitrary initial distribution function. Note that the first
term in (160) carries much of the information of the initial conditions at
early times, but at times greater than $\tau_{R}$, this first term decays
exponentially. Thus the second term describes the late-time behavior of the
system. Finally, note that for conformal RTA, we have Denicol:2010xn ;
Denicol:2011fa

$\displaystyle\tau_{\rm eq}=\frac{5{\eta}}{sT(\tau)}.$ (162)

To introduce anisotropy, one can use the Romatschke-Strickland distribution
function Romatschke:2003ms  with anisotropy parameter, $\xi$,

$\displaystyle f=f_{eq}\left(\sqrt{p^{2}_{\perp}+(1+\xi)p_{z}^{2}}\right),$
(163)

where without loss of generality, we aligned the anisotropy in the $z$
direction. This was studied in the RTA and with a $\lambda\phi^{4}$ scalar
collision kernel Almaalol:2018ynx , which was extended further by considering
dynamical fugacity in the RTA Strickland:2019hff . Interestingly, models with
number-conserving kernels were found to approach the attractor more rapidly,
especially the lowest lying mode, although that was not the case for the
aHydro attractor Almaalol:2018jmz .

One proceeds to solve the equations with the integral iterative method, first
outlined in this context in Florkowski:2013lza ; Florkowski:2013lya . It
essentially boils down to integrating both sides of (160) over the momentum
squared to arrive at an expression for the energy density. Using the matching
conditions, one knows that $\varepsilon=\varepsilon_{\rm eq}$. Then one can
read off the temperature from $\varepsilon\sim T^{4}$. With the temperature,
the equilibrium distribution function is known at all times and so the
integration can be done in (160) to find the distribution function, $f$.

###  4.3 Adiabatic hydrodynamization

Another interesting approach was to study a system with a gapped Hamiltonian,
where the RTA Boltzmann equation of the momentum-weighted distribution
function (integrating out the modulus $p=\sqrt{p^{2}_{T}+p^{2}_{z}}$, but not
performing the angular integration) is recast as a time dependent Schrödinger
equation Brewer:2019oha . The lowest eigenmode is identified as the pre-
hydrodynamic mode, which eventually evolves into the hydrodynamic mode,
remaining gapped from the faster modes throughout the evolution. In this set-
up, the approach to the attractor is mainly characterized by such pre-
hydrodynamic modes, which the authors of Brewer:2019oha  refer to as adiabatic
hydrodynamization. This is due to adiabaticity, i.e. the excitations of the
pre-hydrodynamic modes are suppressed. Future work involves relaxing the
assumption of the RTA to a more realistic case, e.g. a Fokker-Planck collision
kernel. To make better contact with the QGP, it is useful to go beyond the
Bjorken expansion by considering transverse modes. A possible experimental
signature could be found in the flow coefficients $v_{n}$ of small systems,
where the pre-hydrodynamic modes might not have the chance to fully
hydrodynamize.

###  4.4 Hydrodynamic generators

An intriguing notion that has been raised is the existence of hydrodynamic
generators McNelis:2020jrn ; McNelis:2021zsp , where the exact solution of the
RTA Boltzmann equation at late times generates the Borel-resummed Chapman-
Enskog expansion. The Chapman-Enskog expansion is an expansion in the Knudsen
number, $K_{N}$, which is a dimensionless parameter given by the ratio of the
mean free path and a characteristic macroscopic length scale, particular to
the system at hand. For the RTA Boltzmann equation (159), the Chapman-Enskog
expansion is McNelis:2020jrn

$\displaystyle
f_{CE}(\tau,p)=\sum_{n=0}^{\infty}(-\tau_{R}(\tau)\partial_{\tau})^{n}f_{eq}(\tau,p).$
(164)

This series was found to be factorially divergent PhysRevLett.56.1571 ;
Denicol:2016bjh  even for small Knudsen number. One can try to Borel resum the
Chapman-Enskog series as outlined in Sec. 3.1.1 to find

$\displaystyle
f_{CE}(\tau,p)=\int_{0}^{\infty}dze^{-z}\sum_{n=0}^{\infty}\frac{z^{n}}{n!}(-\tau_{R}(\tau)\partial_{\tau})^{n}f_{eq}(\tau,p).$
(165)

Using the exact solution (160) in the case the relaxation time is kept
constant, after some algebra one finds

$\displaystyle
f(\tau,p)=e^{-z_{0}}f_{0}(\tau_{0},p)+\int_{0}^{z_{0}}dze^{-z}\sum_{n=0}^{\infty}\frac{(-z\tau_{R})^{n}}{n!}\partial_{\tau}^{n}f_{eq}(\tau,p),$
(166)

where $z_{0}=({\tau-\tau_{0}})/{\tau_{R}}.$ Clearly, expanding (166) for large
proper time reduces to (165). From here, the authors conjecture that even in
the case the relaxation time is proper time dependent, the hydrodynamic
generator is given by McNelis:2020jrn

$\displaystyle
f_{G}(\tau,p)=\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}D(\tau,\tau^{\prime})f_{\rm
eq}(\tau^{\prime},p)}{\tau_{r}(\tau^{\prime})},$ (167)

where $D(\tau,\tau^{\prime})$ is the damping factor (161). The hydrodynamic
generator is a nice representation of the Chapman-Enskog expansion as it is
finite, is amenable to generalizations to $(3+1)$ dimensions and to systems
without Bjorken flow, and satisfies the RTA Boltzmann equation

$\displaystyle\partial_{\tau}f_{G}=\frac{f_{\rm eq}-f}{\tau_{R}(\tau)}.$ (168)

Interestingly, the RTA hydrodynamic generator contains the Chapman-Enskog
expansion with coupled non-hydrodynamic modes with different decay times. It
is found that the higher order gradient corrections are suppressed during
hydrodynamization, which indicates that these higher corrections do not
necessarily play a large role, even with large gradients. This story seems to
resemble the modes found via resurgence in Heller:2016rtz ; Heller:2018qvh .

Finally, it should be pointed out that the RTA has its limitations. The
derivation above rested on the knowledge of an exact solution in the RTA, so
incorporating more realistic collision kernels will require the damping
function (161) to be replaced more generally with a Green’s function, see e.g.
Kamata:2020mka .

##  5 An eye to phenomenology

We have seen that hydrodynamic attractors are found in many systems. The study
of these attractors offer plenty of theoretical insights, e.g. to the
applicability of hydrodynamics, resurgence and transseries, but what about
about real world applications? In this section, we highlight some of the
phenomenology in heavy ion physics that attractors played an important part.

###  5.1 KøMPøST

Figure 8: Evolution of the background energy density divided by the ideal
hydrodynamic energy density (177) as a function of proper time (normalized by
the kinetic relaxation time $\tau_{R}\equiv(\eta/s)/T_{\rm id})$. Note that
attractor curve goes smoothly from the early time free streaming limit to
viscous hydrodynamics at late times. Figure from Kurkela:2018wud .

KøMPøST is a framework that matches the early time, nonequilibrium QCD kinetic
theory to a later time when a hydrodynamic description is more appropriate
Kurkela:2018wud ; Kurkela:2018vqr . The initial energy momentum tensor is
described by kinetic theory beginning at the time scale, $\tau_{\rm EKT}$.
This is then propagated to later times when hydrodynamics is applicable,
$\tau_{\rm hydro}$, via nonequilibrium Green’s functions,
$G^{\mu\nu}_{\alpha\beta}$, which are computed in QCD kinetic theory. The
linearized perturbations of the stress tensor

$\displaystyle T^{\mu\nu}(\tau_{\rm
EKT},\textbf{x})=\bar{T}^{\mu\nu}(\tau_{\rm EKT})+\delta T^{\mu\nu}(\tau_{\rm
EKT},\textbf{x}),$ (169)

evolve according to

$\displaystyle\delta T^{\mu\nu}(\tau_{\rm hydro},\textbf{x})=\int
d^{2}x^{\prime}G^{\mu\nu}_{\alpha\beta}(\textbf{x},\textbf{x}^{\prime},\tau_{\rm
hydro},\tau_{\rm EKT})\delta T^{\alpha\beta}(\tau_{\rm
EKT},x^{\prime})\frac{\bar{T}^{\tau\tau}(\tau_{\rm
hydro})}{\bar{T}^{\tau\tau}(\tau_{\rm EKT})}.$ (170)

The barred quantities denote the background energy momentum tensor, which is
computed by taking a spatial average over the causal circle. As observed in
Keegan:2015avk ; Heller:2016rtz ; Kurkela:2018wud , the dynamics in these
timescales becomes essentially independent of the microscopic kinetic coupling
when measured in (macroscopic) dimensionless time

$\displaystyle\bar{\omega}=\frac{\tau T}{4\pi\eta/s}.$ (171)

Even far from equilibrium for $\bar{\omega}<1$, it is observed that the energy
momentum tensor is a function of $\bar{\omega}$ only with hydrodynamization
occuring after $\bar{\omega}>1.$ As such, the attractor can be viewed as a
bridge between the kinetic and hydrodynamic regimes. The universal attractor
behavior can be seen in Fig. 8, where two values of the kinetic theory
coupling is plotted for a range of $\eta/s$. The nonequilibrium evolution of
the system goes from early time free streaming to late time viscous
hydrodynamics. The universal scaling function is due to the attractor
behavior.

KøMPøST has been incorporated in a number of contexts, e.g. to pre-
hydrodynamic evolution and its signatures in final-state heavy-ion observables
NunesdaSilva:2020bfs  and more recently, to describe photon observables in
small collision systems Gale:2021zlc .

###  5.2 Entropy production

The fact that hydrodynamic attractors can serve as a bridge between late-time
and earlier, far-from-equilibrium behavior was recently exploited to extract
phenomenologically relevant information in the QGP context Giacalone:2019ldn .
The central idea is that the multiplicity of final state particles measured in
detectors is intimately linked to the entropy production during the pre-
equilibrium phase. The total entropy per unit spacetime rapidity, $\eta_{s},$
during the near-equilbrium transverse expansion of the QGP is

$\displaystyle\frac{dS}{d\eta_{s}}=A_{\perp}(s\tau)_{\rm hydro},$ (172)

which is related to the charged particle multiplicity at late times via

$\displaystyle\frac{dN_{\rm
ch}}{d\eta_{s}}\frac{dS}{dy}\approx\frac{1}{J}A_{\perp}(s\tau)_{\rm
hydro}\frac{dN_{\rm ch}}{dy},$ (173)

where $J$ denotes the Jacobian factor, $y$ is the rapidity, $dy/d\eta\approx
1.1$, $A_{\perp}=\pi R^{2}$ and $2R$ is the transverse extent of the system.
Then the initial energy density can be linked to the experimentally observed
particle multiplicities via the neat pocket formula Giacalone:2019ldn

$\displaystyle(s\tau)_{\rm
hydro}=\frac{4}{3}C_{\infty}^{3/4}\left(4\pi\frac{\eta}{s}\right)^{1/3}\left(\frac{\pi^{2}}{30}\nu_{\rm
eff}\right)^{1/3}(\varepsilon\tau)_{0}^{2/3},$ (174)

where $C_{\infty}$ is a constant of order unity and $\nu_{\rm eff}$ is the
number of effective degrees of freedom. For an ideal gas of gluons, $\nu_{\rm
eff}=16$, while for a mixture of quarks and gluons, $\nu_{\rm eff}=47.5$.

To derive this result, the knowledge that there is an attractor can be put to
work to find

$\displaystyle\frac{\varepsilon(\tau)\tau^{4/3}}{\varepsilon_{\rm
hydro}\tau_{\rm hydro}^{4/3}}=\mathcal{E}(\bar{\omega}),$ (175)

where the quantity on the right hand side is the attractor curve. To see this,
we follow the derivation of the link between the initial state energy density,
$\varepsilon_{0}$, and the near thermal system energy density,
$\varepsilon_{\rm hydro}$, beginning from the boost-invariant conservation law

$\displaystyle\partial_{\tau}\varepsilon=-\frac{\varepsilon+P_{L}}{\tau}.$
(176)

The energy density is related to the effective temperature via

$\displaystyle\varepsilon(\tau)=\frac{\pi^{2}}{30}\nu_{\rm eff}T^{4}_{\rm
eff}(\tau),$ (177)

It is then convenient to introduce the dimensionless scaling variable
$\bar{\omega}$ given by (171). Finally, calling the anisotropy
$f(\bar{\omega})\equiv P_{L}/\varepsilon$, we see that we can rewrite (176) as

$\displaystyle\frac{d\bar{\omega}}{d\tau}\partial_{\bar{\omega}}\log\varepsilon=-\frac{1+f}{\tau}.$
(178)

Using the definition of $\bar{\omega}$ and the relation between the energy
density and the effective temperature (177), we see that

$\displaystyle\frac{4\pi\eta}{s}\frac{d\bar{\omega}}{d\tau}$
$\displaystyle=T_{\rm eff}+\tau\partial_{\tau}T_{\rm eff},$ (179)
$\displaystyle=T_{\rm eff}+\frac{1}{4}T_{\rm
eff}(\tau\varepsilon^{-1}\partial_{\tau}\varepsilon),$ (180)
$\displaystyle=T_{\rm eff}\left(\frac{3}{4}-\frac{1}{4}f\right).$ (181)

Thus,

$\displaystyle\frac{\partial}{\partial\bar{\omega}}\log\varepsilon=-\frac{1+f}{\frac{\tau
T_{\rm
eff}}{4\pi\eta/s}\left(\frac{3}{4}-\frac{1}{4}f\right)}=-\frac{1+f}{\bar{\omega}\left(\frac{3}{4}-\frac{1}{4}f\right)}.$
(182)

After integration, we arrive at

$\displaystyle\varepsilon(\tau)=\varepsilon_{0}\exp\left(-\int_{\bar{\omega_{0}}}^{\bar{\omega}_{\rm
hydro}}\frac{d\bar{\omega}}{\bar{\omega}}\frac{1+f}{\frac{3}{4}-\frac{1}{4}f}\right).$
(183)

From here, we arrive straightforwardly at (175). Note that the energy
attractor (175) interpolates smoothly between free streaming and late time
viscous hydrodynamics

$\displaystyle\mathcal{E}(\bar{\omega}\ll
1)=C_{\infty}^{-1}\bar{\omega}^{4/9},$ (184)
$\displaystyle\mathcal{E}(\bar{\omega}\gg 1)=1-\frac{2}{3\pi\bar{\omega}},$
(185)

respectively. The pocket formula (174) and the charged particle multiplicity
(173) can be used to estimate the initial energy density for central lead
collisions at $\sqrt{s_{NN}}=2.76$ TeV to be $\varepsilon\approx 270$ GeV/fm3
at the time $\tau_{0}=0.1$ fm/c Giacalone:2019ldn . Furthermore, the study of
such energy attractors has been used to characterize potential probes of the
QGP at early times, such as intermediate mass dileptons Coquet:2021cuv .
Interestingly, the evolution of a system following an attractor solution has
been shown to result in higher yields (in fact, maximum yields) in thermal
dileption and photon particle production Naik:2021yph , compared to solutions
away from the attractor.

The previous argument was turned on its head in Jankowski:2020itt , where the
assumption of free streaming at early times was explored by positing the
energy density scales with some to-be-determined parameter, $\beta$, namely
$\varepsilon\sim\tau^{-\beta}$ (the free streaming case has $\beta=1$). Using
different initial state models and comparing the final state to experimental
data, they find that the preferred expansion parameter $\beta$ is slightly,
but distinctly different from the free streaming case.

###  5.3 Non-zero chemical potential

Real world QGP has non-zero baryon chemical potential. As such, it is natural
to include the effect of non-zero chemical potential on the attractor story
described so far. This has been explored both in hydrodynamic Dore:2020fiq ;
Dore:2020jye  and kinetic setups Du:2020zqg ; Du:2020dvp . In the hydrodynamic
case, the $0+1d$ DNMR Bjorken case described in Sec. 2.2 is supplemented with
an evolution of the baryon number density and relaxation of the bulk
viscosity. The late-time attractor regime was not reached due to the short
lifetime of the hydrodynamic runs, so it is more apt to talk of potential
attractors, which the system looked like it was approaching. Although
attractors were not systematically studied, there was the indication that near
the critical point, the bulk attractor should have a pronounced effect. Also,
the potential attractor line has a different behavior when comparing MIS to
DNMR near the critical point, indicating care should be taken as there are
phenomenological implications and differing results depending on which set of
equations of motion is used Dore:2020jye .

In the kinetic approach of Du:2020zqg ; Du:2020dvp , quark/anti-quark and
gluonic distribution functions are taken into account, which respectively take
the form

$\displaystyle f_{q_{f}/\bar{q}_{f}}$
$\displaystyle=f_{q_{f}/\bar{q}_{f}}^{0}\frac{\sqrt{p_{T}^{2}+p_{\parallel}^{2}}}{\sqrt{p_{T}^{2}+\xi_{0}^{2}p_{\parallel}^{2}}}e^{-2(p_{T}^{2}+\xi_{0}^{2}p_{\parallel}^{2})/(3Q_{0}^{2})},$
(186) $\displaystyle f_{g}$
$\displaystyle=\frac{{f_{g}^{0}Q_{0}}}{\sqrt{p_{T}^{2}+\xi_{0}^{2}p_{\parallel}^{2}}}e^{-2(p_{T}^{2}+\xi_{0}^{2}p_{\parallel}^{2})/(3Q_{0}^{2})},$
(187)

where $\xi_{0}$ is the initial anisotropy and $Q_{0}$ is the initial energy
scale. The evolution of the anisotropic distribution functions with form
similar to (163) is given by the QCD Boltzmann equation Mueller:1999pi  with a
collision kernel with leading order elastic ($2\leftrightarrow 2$) and
inelastic ($1\leftrightarrow 2$) interactions taken into account, i.e. for a
distribution function of species $a$, $f_{a}=f_{a}(\tau,p_{T},p_{\parallel})$

$\displaystyle\left[\partial_{\tau}-\frac{p_{\parallel}}{\tau}\frac{\partial}{\partial
p_{\parallel}}\right]f_{a}=-C^{2\leftrightarrow 2}[f]-C^{2\leftrightarrow
1}[f].$ (188)

The evolution of the system is given by the conservation of the energy
momentum tensor and the conserved current, $\Delta J^{\mu}_{f}=(\Delta
n_{f},0,0,0)$, via

$\displaystyle 0$
$\displaystyle=\partial_{\tau}\varepsilon+\frac{\varepsilon+p_{L}}{\tau},$
(189) $\displaystyle 0$ $\displaystyle=\partial_{\tau}\Delta
n_{f}+\frac{\Delta n_{f}}{\tau}.$ (190)

The nonequilibrium system has a build up of longitudinal pressure due to
interactions while the longitudinal expansion slows down, before finally at
late times approaching local equilibrium with $p_{L}=p_{T}=\varepsilon/3$.

Like in the case of zero chemical potential discussed previously, the validity
of viscous hydrodynamics as a description of the non-equilibrium evolution in
the case of non-zero chemical potential also begins at timescales roughly
given by $4\pi\eta/(sT)$ Du:2020zqg ; Du:2020dvp . However, at early times
non-zero chemical potential can have a large impact on the pressure for times
earlier than $4\pi\eta/(sT)$, although the late time approach is found to be
affected by around $~{}10\%$. This is to be contrasted with results from
kinetic theory computations involving a single component distribution
function, as discussed in Sec. 4, where the rapid memory loss of the initial
conditions characteristic to attractor behavior allows for a universal
approach to viscous hydrodynamics. As such, it is clear that the chemical
composition of the QGP has an effect on the approach to the universal approach
to viscous hydrodynamics at late time.

###  5.4 Improved freeze-out procedures

It is worthwhile to keep in mind that for all the discussion of hydrodynamic
variables as a means to describe heavy ion collisions, we have left out a
crucial ingredient. Experiments themselves measure the distribution of
particles and not the fluid-dynamic variables directly as discussed
previously. Following the collision and initial evolution of the quark gluon
plasma, the distribution of particles “freeze-out” as color neutral hadrons
and free stream towards the detectors at time scales around $\tau\sim 10\text{
fm/c}$. The translation of distribution of particles to hydrodynamic variables
is known as the freeze-out procedure.

Recently, it was pointed out in Almaalol:2020rnu  that possible improvements
to freeze-out procedures can be made in light of the developments surrounding
the attractor. Using Arnold-Moore-Yaffe (AMY) effective kinetic theory in one-
dimensional Bjorken flow, where the collision kernel includes gluon scattering
and inelastic gluon splitting using the pure glue code from
PhysRevLett.115.182301 , a comparison was made between aHydro and linearized
freeze-out prescriptions. It was found that the freeze-out prescriptions
studied tend to reproduce the low-order moments of the distribution at late
times $\tau\gg 5\text{ fm/c}.$ However, for higher moments, it was found that
the linearized prescriptions do poorly compared to aHydro. The last fact
indicates that the aHydro freeze-out prescription can be of use for the study
of smaller systems, like $pA$ or peripheral $AA$ collisions.

###  5.5 Beyond conformal symmetry

The majority of the work discussed so far centered on conformal systems, which
provides the advantage of theoretical clarity and sometimes the feasibility of
analytical computations. However, in more realistic systems, one should not
expect full conformal symmetry to hold, especially with a full three
dimensional expansion and a realistic equation of state. There are a number of
works that go beyond the status of the field, including the exploratory work
of attractors with a non-conformal equation of state relevant for the QCD
critical point considered in Dore:2020fiq ; Dore:2020jye , see Sec. 5.3 for
more details.

Here we discuss an initial survey that was made in Romatschke:2017acs  for a
$0+1$-dimensional Bjorken kinetic theory of a gas of particles of mass $m$ and
for inhomogeneous $2+1$ dimensional BRSSS hydrodynamics. Two quantities were
proposed to investigate whether attractor solutions are present, namely

$\displaystyle A_{1}$
$\displaystyle=\frac{u^{\mu}\nabla_{\mu}\varepsilon}{(\varepsilon+P)\nabla_{\lambda}u^{\lambda}},$
(191) $\displaystyle A_{2}$
$\displaystyle=\frac{\text{Tr}T^{\mu\nu}+\varepsilon-3P}{\zeta T}.$ (192)

The first quantity represents a simple rewriting of the evolution of the
energy density from the Navier-Stokes equations, while the second denotes the
deviation of the equation of state from the conformal limit. When expressed as
a function of the inverse gradient strength $\Gamma$,

$\displaystyle\Gamma=\left[\frac{\eta}{2s}\frac{\sigma^{\mu\nu}\sigma_{\mu\nu}}{T\nabla_{\lambda}u^{\lambda}}+\frac{\zeta}{s}\frac{\nabla_{\lambda}u^{\lambda}}{T}\right]^{-1},$
(193)

both $A_{1}$ and $A_{2}$ exhibit attractor behavior. However, the approach of
different solutions is slower for $A_{2}$ than $A_{1}$, which may indicate
that at early times, the attractor is harder to discern in non-conformal
systems. In the case of the $2+1$ dimensional BRSSS hydrodynamics, the
apparent attractor was seen by studying $A_{1}(\Gamma)$, which seemed to
greatly depart from the expected Navier-Stokes limit, although such a
simulation proved too computationally costly to extend for much later times.

In many respects, although the exploration of the non-conformal limit is just
at beginning, there are a number of lessons that have been learned already. As
mentioned above, moving away from the conformal limit also can provide insight
into the nature of the early time attractor. Recent work involved studying a
$0+1$ dimensional gas of massive particles evolving under the Boltzmann
equation in the RTA, which was compared to non-conformal hydrodynamics
Chattopadhyay:2021ive ; Jaiswal:2021uvv . It is found that only the
longitudinal pressure has an early time attractor, while that there is no
early time hydrodynamic attractor for the bulk and shear viscous stresses,
adding evidence to the related indication in Kurkela:2019set . Moreover, the
comparison of the quantity (191) in Jaiswal:2021uvv  indicates that while for
late times $A_{1}$ does exhibit universal attractor behavior, at early times
there is a noticeable difference between the solutions depending on the choice
of different masses. Another interesting observation in non-conformal Bjorken
flow is that it is not a guarentee one even reaches local thermal equilibrium
at late times. For instance, in the case of variable relaxation time scaling,
namely $\tau_{R}\propto T^{-\Delta}$ ($\Delta=1$ corresponds to the conformal
limit), the system can evolve to a late time regime out of local thermal
equilibrium for $\Delta\neq 1$ Heller:2018qvh , in contrast to conformal
Bjorken flow.

##  6 Conclusion and Outlook

The story of hydrodynamic attractors in heavy ion collisions is still being
written. We can take this point to reflect on where the field stands and where
it might go. At this stage, it seems like attractors for Bjorken flow are
well-established for a variety of effective theories relevant for heavy ion
collisions, such as relativistic hydrodynamics, holography and kinetic theory.
Attractors provide a conceptual mechanism to understand the applicability of
hydrodynamics even in far-from-equilibrium systems, as well as providing some
phenomenological insight. Notwithstanding the current successes, there is
still much work to be done. Some open questions remain:

  * •

What are the theoretical properties of the QCD attractor and under which
conditions does QCD exhibit attractor behaviour?

  * •

How do attractors behave in systems with relaxed symmetries? Studies of
attractors have begun to go beyond the Bjorken and Gubser cases described
here. Ultimately, the goal should be to consider the full symmetries in heavy
ion collisions, moving away from the conformal limit as in Romatschke:2017acs
; Chattopadhyay:2021ive ; Jaiswal:2021uvv ; Chen:2021wwh  and including
transverse expansion in $3+1$ dimensions, as well as the quark degrees of
freedom as in Kurkela:2018xxd ; Kurkela:2018oqw ; Dore:2020jye ; Du:2020zqg ;
Du:2020dvp . There is some promising work to take transverse dynamics into
account Ambrus:2021sjg . Recently, the attractor was studied in a Hubble
expanding background Du:2021fok , where the scale factor was assumed to follow
a power law form $a(t)\sim t^{\alpha}$ instead of solving the Einstein
equations to find $a(t)$, an assumption informed by cosmology.

  * •

Beyond numerical techniques, how can we extend the determination of attractors
to earlier times? Especially in the more complicated setting of relaxed
symmetries, it is not clear if there is an unambiguous early time attractor.
Also, in the holographic context, it seems like there is only the late-time
attractor Kurkela:2019set  (although the next question also can be relevant
for the early time attractor).

  * •

How would finite coupling corrections affect the story of holographic
attractors? The Gauss-Bonnet setup described in Sec. 3.3 can act as a
laboratory to study higher derivative corrections. The observation that
varying a single parameter allows the interpolation from the strongly coupled
to the expectation from weakly coupled kinetic theory should be viewed as an
encouraging step in this direction. It can be speculated that finite coupling
corrections thus might provide an unambiguous method of determining the
possible early time behavior of the attractor.

  * •

How do hydro attractors fit in with non-thermal fixed points? It would be
interesting to see a study of the evolution of a kinetic model from a non-
thermal fixed point to a late time hydrodynamic attractor.

  * •

Are there any signatures of attractors in smaller systems, e.g. $pp$
collisions? In $pp$ collisions, the number of particles produced is regularly
quite low compared to typical heavy ion collisions. It is not clear if
hydrodynamization would occur in such short-lived settings Wiedemann:2021wju .
If so, there could still be the approach to the hydrodynamic attractor for
generic initial conditions without reaching the attractor at all, potentially
indicating some level of pre-hydrodynamization in the smaller systems.

  * •

Which observables are sensitive to hydrodynamic attractors? Attractors seem to
be a useful tool in using final state information to probe earlier times in
heavy ion collisions, see Sec. 5. Furthermore, one can envisage using them to
determine which final state observables depend directly on the initial state.

  * •

How universal are hydrodynamic attractors? The hydrodynamic limit exists in a
number of neighboring fields, including cosmology and condensed matter
systems. The lessons learned about attractors in the high energy community can
surely be applied to expanding systems in a variety of fields.

### Acknowledgements

I would like to extend my gratitude to Kirill Boguslavski, Eduardo Grossi,
Michal Heller, Aleksi Kurkela, Mauricio Martinez, Aleksas Mazeliauskas, Ayan
Mukhopadhyay, David Müller, Jorge Noronha, Anton Rebhan, Paul Romatschke,
Alexandre Serantes, Soeren Schlichtling, Michał Spaliński, Michael Strickland,
Victor Svensson, Derek Teaney, Urs Wiedemann and Yi Yin for helpful
discussions. Many thanks to Eduardo Grossi, Ayan Mukhopadhyay and David Müller
for comments on the manuscript. This work is supported by the Erwin
Schrödinger Fellowship of the Austrian Science Fund (FWF), project no. J4406.

##  Appendix A List of abbreviations

3+1D

3+1 dimensions / dimensional

AdS/CFT

Anti-de Sitter/Conformal Field Theory

AMY

Arnold-Moore-Yaffe

BGK

Bhatnagar-Gross-Krook

BRSSS

Baier-Romatschke-Son-Starinets-Stephanov

DNMR

Denicol-Niemi-Molnar-Rischke

GB

Gauss-Bonnet

HJSW

Heller-Janik-Spalinski-Witaszczyk

KøMPøST

Kurkela-Mazeliauskas-Paquet-Schlichtling-Teaney

LHC

Large Hadron Collider

MIS

Müller-Israel-Stewart

QCD

quantum chromodynamics

QGP

quark-gluon plasma

RHIC

Relativistic Heavy Ion Collider

RTA

relaxed time approximation

SYM

Super Yang-Mills

## References

  * (1) P. Kovtun, _Lectures on hydrodynamic fluctuations in relativistic theories_ , _J. Phys. A_ 45 (2012) 473001, [1205.5040].
  * (2) P. Romatschke and U. Romatschke, _Relativistic Fluid Dynamics In and Out of Equilibrium_.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 5, 2019, 10.1017/9781108651998.
  * (3) C. Shen and L. Yan, _Recent development of hydrodynamic modeling in heavy-ion collisions_ , _Nucl. Sci. Tech._ 31 (10, 2020) 122, [2010.12377].
  * (4) U. Heinz and R. Snellings, _Collective flow and viscosity in relativistic heavy-ion collisions_ , _Ann. Rev. Nucl. Part. Sci._ 63 (2013) 123–151, [1301.2826].
  * (5) P. Chaikin and T. Lubensky, _Principles of Condensed Matter Physics_.  Cambridge University Press, 2000.
  * (6) S. N. Shore, _An Introduction to Astrophysical Hydrodynamics_.  Academic Press.
  * (7) S. Weinberg, _Cosmology_.  Oxford University Press, 2008.
  * (8) D. A. Teaney, _Viscous Hydrodynamics and the Quark Gluon Plasma_ , pp. 207–266.  2010\.  0905.2433.  10.1142/9789814293297_0004.
  * (9) R. Snellings, _Elliptic Flow: A Brief Review_ , _New J. Phys._ 13 (2011) 055008, [1102.3010].
  * (10) P. Romatschke and U. Romatschke, _Viscosity Information from Relativistic Nuclear Collisions: How Perfect is the Fluid Observed at RHIC?_ , _Phys. Rev. Lett._ 99 (2007) 172301, [0706.1522].
  * (11) G. Policastro, D. T. Son and A. O. Starinets, _The Shear viscosity of strongly coupled N=4 supersymmetric Yang-Mills plasma_ , _Phys. Rev. Lett._ 87 (2001) 081601, [hep-th/0104066].
  * (12) P. Kovtun, D. T. Son and A. O. Starinets, _Holography and hydrodynamics: Diffusion on stretched horizons_ , _JHEP_ 10 (2003) 064, [hep-th/0309213].
  * (13) P. Kovtun, D. T. Son and A. O. Starinets, _Viscosity in strongly interacting quantum field theories from black hole physics_ , _Phys. Rev. Lett._ 94 (2005) 111601, [hep-th/0405231].
  * (14) J. Berges, M. P. Heller, A. Mazeliauskas and R. Venugopalan, _QCD thermalization: Ab initio approaches and interdisciplinary connections_ , _Rev. Mod. Phys._ 93 (2021) 035003, [2005.12299].
  * (15) W. Florkowski, M. P. Heller and M. Spalinski, _New theories of relativistic hydrodynamics in the LHC era_ , _Rept. Prog. Phys._ 81 (2018) 046001, [1707.02282].
  * (16) M. P. Heller and M. Spalinski, _Hydrodynamics Beyond the Gradient Expansion: Resurgence and Resummation_ , _Phys. Rev. Lett._ 115 (2015) 072501, [1503.07514].
  * (17) M. P. Heller, R. A. Janik and P. Witaszczyk, _The characteristics of thermalization of boost-invariant plasma from holography_ , _Phys. Rev. Lett._ 108 (2012) 201602, [1103.3452].
  * (18) M. P. Heller, R. A. Janik and P. Witaszczyk, _Hydrodynamic Gradient Expansion in Gauge Theory Plasmas_ , _Phys. Rev. Lett._ 110 (2013) 211602, [1302.0697].
  * (19) W. Florkowski, E. Maksymiuk and R. Ryblewski, _Coupled kinetic equations for fermions and bosons in the relaxation-time approximation_ , _Phys. Rev. C_ 97 (2018) 024915, [1710.07095].
  * (20) P. Romatschke, _Relativistic Hydrodynamic Attractors with Broken Symmetries: Non-Conformal and Non-Homogeneous_ , _JHEP_ 12 (2017) 079, [1710.03234].
  * (21) A. Buchel, M. P. Heller and J. Noronha, _Entropy Production, Hydrodynamics, and Resurgence in the Primordial Quark-Gluon Plasma from Holography_ , _Phys. Rev. D_ 94 (2016) 106011, [1603.05344].
  * (22) J. Casalderrey-Solana, N. I. Gushterov and B. Meiring, _Resurgence and Hydrodynamic Attractors in Gauss-Bonnet Holography_ , _JHEP_ 04 (2018) 042, [1712.02772].
  * (23) M. Rangamani, _Gravity and Hydrodynamics: Lectures on the fluid-gravity correspondence_ , _Class. Quant. Grav._ 26 (2009) 224003, [0905.4352].
  * (24) G. Giacalone, A. Mazeliauskas and S. Schlichting, _Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions_ , _Phys. Rev. Lett._ 123 (2019) 262301, [1908.02866].
  * (25) J. Jankowski, S. Kamata, M. Martinez and M. Spaliński, _Constraining the initial stages of ultrarelativistic nuclear collisions_ , _Phys. Rev. D_ 104 (2021) 074012, [2012.02184].
  * (26) T. Dore, E. McLaughlin and J. Noronha-Hostler, _Far From Equilibrium Hydrodynamics and the Beam Energy Scan_ , _J. Phys. Conf. Ser._ 1602 (2020) 012017, [2006.04206].
  * (27) T. Dore, J. Noronha-Hostler and E. McLaughlin, _Far-from-equilibrium search for the QCD critical point_ , _Phys. Rev. D_ 102 (2020) 074017, [2007.15083].
  * (28) A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting and D. Teaney, _Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory_ , _Phys. Rev. Lett._ 122 (2019) 122302, [1805.01604].
  * (29) A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting and D. Teaney, _Effective kinetic description of event-by-event pre-equilibrium dynamics in high-energy heavy-ion collisions_ , _Phys. Rev. C_ 99 (2019) 034910, [1805.00961].
  * (30) G. Basar and G. V. Dunne, _Hydrodynamics, resurgence, and transasymptotics_ , _Phys. Rev. D_ 92 (2015) 125011, [1509.05046].
  * (31) I. Aniceto, G. Basar and R. Schiappa, _A Primer on Resurgent Transseries and Their Asymptotics_ , _Phys. Rept._ 809 (2019) 1–135, [1802.10441].
  * (32) A. Behtash, C. N. Cruz-Camacho and M. Martinez, _Far-from-equilibrium attractors and nonlinear dynamical systems approach to the Gubser flow_ , _Phys. Rev. D_ 97 (2018) 044041, [1711.01745].
  * (33) A. Behtash, S. Kamata, M. Martinez and H. Shi, _Dynamical systems and nonlinear transient rheology of the far-from-equilibrium Bjorken flow_ , _Phys. Rev. D_ 99 (2019) 116012, [1901.08632].
  * (34) A. Behtash, S. Kamata, M. Martinez and H. Shi, _Global flow structure and exact formal transseries of the Gubser flow in kinetic theory_ , _JHEP_ 07 (2020) 226, [1911.06406].
  * (35) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Turbulent thermalization process in heavy-ion collisions at ultrarelativistic energies_ , _Phys. Rev. D_ 89 (2014) 074011, [1303.5650].
  * (36) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Universal attractor in a highly occupied non-Abelian plasma_ , _Phys. Rev. D_ 89 (2014) 114007, [1311.3005].
  * (37) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Basin of attraction for turbulent thermalization and the range of validity of classical-statistical simulations_ , _JHEP_ 05 (2014) 054, [1312.5216].
  * (38) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Universality far from equilibrium: From superfluid Bose gases to heavy-ion collisions_ , _Phys. Rev. Lett._ 114 (2015) 061601, [1408.1670].
  * (39) K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron, _Highly occupied gauge theories in 2+1 dimensions: A self-similar attractor_ , _Phys. Rev. D_ 100 (2019) 094022, [1907.05892].
  * (40) A. Kurkela, W. van der Schee, U. A. Wiedemann and B. Wu, _Early- and Late-Time Behavior of Attractors in Heavy-Ion Collisions_ , _Phys. Rev. Lett._ 124 (2020) 102301, [1907.08101].
  * (41) W. A. Hiscock and L. Lindblom, _Generic instabilities in first-order dissipative relativistic fluid theories_ , _Phys. Rev. D_ 31 (1985) 725–733.
  * (42) P. Kostadt and M. Liu, _Causality and stability of the relativistic diffusion equation_ , _Phys. Rev. D_ 62 (2000) 023003, [cond-mat/0010276].
  * (43) I. Müller, _Zum Paradoxon der Wärmeleitungstheorie_ , _Z. Phys._ 198 (1967) 329–344.
  * (44) W. Israel, _Nonstationary irreversible thermodynamics: A Causal relativistic theory_ , _Annals Phys._ 100 (1976) 310–331.
  * (45) W. Israel and J. Stewart, _Transient relativistic thermodynamics and kinetic theory_ , _Annals Phys._ 118 (1979) 341–372.
  * (46) R. Baier, P. Romatschke, D. T. Son, A. O. Starinets and M. A. Stephanov, _Relativistic viscous hydrodynamics, conformal invariance, and holography_ , _JHEP_ 04 (2008) 100, [0712.2451].
  * (47) G. S. Denicol, H. Niemi, E. Molnar and D. H. Rischke, _Derivation of transient relativistic fluid dynamics from the Boltzmann equation_ , _Phys. Rev. D_ 85 (2012) 114047, [1202.4551].
  * (48) M. Strickland, J. Noronha and G. Denicol, _Anisotropic nonequilibrium hydrodynamic attractor_ , _Phys. Rev. D_ 97 (2018) 036020, [1709.06644].
  * (49) A. Muronga, _Causal theories of dissipative relativistic fluid dynamics for nuclear collisions_ , _Phys. Rev. C_ 69 (2004) 034903, [nucl-th/0309055].
  * (50) M. Strickland, _Anisotropic Hydrodynamics: Three lectures_ , _Acta Phys. Polon. B_ 45 (2014) 2355–2394, [1410.5786].
  * (51) M. McNelis, D. Bazow and U. Heinz, _(3+1)-dimensional anisotropic fluid dynamics with a lattice QCD equation of state_ , _Phys. Rev. C_ 97 (2018) 054912, [1803.01810].
  * (52) D. Almaalol and M. Strickland, _Anisotropic hydrodynamics with a scalar collisional kernel_ , _Phys. Rev. C_ 97 (2018) 044911, [1801.10173].
  * (53) H. Alalawi and M. Strickland, _Improved anisotropic hydrodynamics ansatz_ , _Phys. Rev. C_ 102 (2020) 064904, [2006.13834].
  * (54) P. Romatschke, _Relativistic (Lattice) Boltzmann Equation with Non-Ideal Equation of State_ , _Phys. Rev. D_ 85 (2012) 065012, [1108.5561].
  * (55) U. Heinz, D. Bazow, G. S. Denicol, M. Martinez, M. Nopoush, J. Noronha et al., _Exact solutions of the Boltzmann equation and optimized hydrodynamic approaches for relativistic heavy-ion collisions_ , _Nucl. Part. Phys. Proc._ 276-278 (2016) 193–196, [1509.05818].
  * (56) M. P. Heller, R. A. Janik, M. Spaliński and P. Witaszczyk, _Coupling hydrodynamics to nonequilibrium degrees of freedom in strongly interacting quark-gluon plasma_ , _Phys. Rev. Lett._ 113 (2014) 261601, [1409.5087].
  * (57) I. Aniceto and M. Spaliński, _Resurgence in Extended Hydrodynamics_ , _Phys. Rev. D_ 93 (2016) 085008, [1511.06358].
  * (58) M. P. Heller, R. Jefferson, M. Spaliński and V. Svensson, _Hydrodynamic Attractors in Phase Space_ , _Phys. Rev. Lett._ 125 (2020) 132301, [2003.07368].
  * (59) A. Nunez and A. O. Starinets, _AdS / CFT correspondence, quasinormal modes, and thermal correlators in N=4 SYM_ , _Phys. Rev. D_ 67 (2003) 124013, [hep-th/0302026].
  * (60) J. D. Bjorken, _Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region_ , _Phys. Rev. D_ 27 (1983) 140–151.
  * (61) A. R. Liddle, P. Parsons and J. D. Barrow, _Formalizing the slow roll approximation in inflation_ , _Phys. Rev. D_ 50 (1994) 7222–7232, [astro-ph/9408015].
  * (62) G. S. Denicol and J. Noronha, _Analytical attractor and the divergence of the slow-roll expansion in relativistic hydrodynamics_ , _Phys. Rev. D_ 97 (2018) 056021, [1711.01657].
  * (63) D. Dorigoni, _An Introduction to Resurgence, Trans-Series and Alien Calculus_ , _Annals Phys._ 409 (2019) 167914, [1411.3585].
  * (64) G. A. Edgar, _Transseries for Beginners_ , _Real Analysis Exchange_ 35 (2009) 253 – 310.
  * (65) I. Aniceto, B. Meiring, J. Jankowski and M. Spaliński, _The large proper-time expansion of Yang-Mills plasma as a resurgent transseries_ , _JHEP_ 02 (2019) 073, [1810.07130].
  * (66) T. Schäfer and E. V. Shuryak, _Instantons in QCD_ , _Rev. Mod. Phys._ 70 (1998) 323–426, [hep-ph/9610451].
  * (67) O. Costin, _Exponential asymptotics, transseries, and generalized Borel summation for analytic rank one systems of ODE’s_ , _arXiv Mathematics e-prints_ (Aug., 2006) math/0608414, [math/0608414].
  * (68) O. Costin and R. D. Costin, _On the formation of singularities of solutions of nonlinear differential systems in antistokes directions_ , _Inventiones Mathematicae_ 145 (Sept., 2001) 425–485, [math/0202234].
  * (69) I. Aniceto and R. Schiappa, _Nonperturbative Ambiguities and the Reality of Resurgent Transseries_ , _Commun. Math. Phys._ 335 (2015) 183–245, [1308.1115].
  * (70) S. S. Gubser, _Symmetry constraints on generalizations of Bjorken flow_ , _Phys. Rev. D_ 82 (2010) 085027, [1006.0006].
  * (71) F. S. Bemfica, M. M. Disconzi and J. Noronha, _Causality and existence of solutions of relativistic viscous fluid dynamics with gravity_ , _Phys. Rev. D_ 98 (2018) 104064, [1708.06255].
  * (72) G. S. Denicol and J. Noronha, _Hydrodynamic attractor and the fate of perturbative expansions in Gubser flow_ , _Phys. Rev. D_ 99 (2019) 116004, [1804.04771].