aainstitutetext: The Institute of Mathematical Sciences, Chennai 600113,
Indiabbinstitutetext: Homi Bhabha National Institute, Training School Complex,
AnushaktiNagar, Mumbai 400094, Indiaccinstitutetext: Center for Strings,
Gravitation and Cosmology, Indian Institute of Technology Madras, Chennai
600036, Indiaddinstitutetext: Institut für Theoretische Physik, Technische
Universität Wien, Wiedner Hauptstr. 8-10, A-1040 Vienna, Austria

# Hydrodynamization in hybrid Bjorken flow attractors

Toshali Mitra c Sukrut Mondkar c Ayan Mukhopadhyay d Anton Rebhan d and
Alexander Soloviev<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Hybrid fluid models, consisting of two sectors with more weakly and more
strongly self-interacting degrees of freedom coupled consistently as in the
semi-holographic framework, have been shown to exhibit an attractor surface
for Bjorken flow. Retaining only the simple viscid fluid descriptions of both
sectors, we can establish analytically the universality of a feature of the
bottom-up thermalization scenario, namely that the energy density of the more
weakly coupled (perturbative) sector at sufficiently early times dominates
over that of the strongly coupled (holographic) one with an exponent
determined by the respective viscosities and relaxation times. On the
attractor surface, the hydrodynamization times of both subsectors decrease
following a conformal scaling behavior with the total energy densities at the
respective hydrodynamization (times a pre-factor), reach their minimum values
and subsequently increase rapidly. The minimum values are obtained when the
respective energy densities are of the order of the inverse of the
dimensionful inter-system coupling. The pre-factor in the conformal window is
universal and identical to that in usual conformal attractors for the dominant
sector, while that of the other sector increases with the ratio of the
dominant to subdominant sector’s energy-densities at sufficiently early time.
Studying initial conditions motivated by saturation scenarios in $p$-$p$ and
Pb-Pb collisions, we find that the perturbative sector hydrodynamizes much
later and the holographic sector hydrodynamizes earlier in $p$-$p$ collisions,
since the total energy densities at the hydrodynamization times fall inside
and outside the conformal window for the two sectors respectively.

###### Keywords:

Hydrodynamic attractor, quark-gluon plasma, hydrodynamization

##  1 Introduction

One of the most remarkable discoveries in non-equilibrium dynamics recently
has been that quantum many-body systems can hydrodynamize far away from
equilibrium Chesler:2009cy ; Chesler:2010bi ; Heller:2011ju ; Chesler:2015bba
; Attems:2016tby ; Attems:2017zam ; Romatschke:2017vte ; Florkowski:2017olj .
Hydrodynamization refers to the stage of dynamical evolution in which the
expectation values of the energy-momentum tensor and conserved currents start
following the constitutive relations of hydrodynamics to a very good
approximation. Remarkably, hydrodynamization can happen in both weakly coupled
and strongly coupled systems even when they are far away from equilibrium and
have large pressure anisotropies Keegan:2015avk . Theoretically,
hydrodynamization can be formalized in terms of a hydrodynamic attractor in
phase space which is built out of a resummation of the divergent series of the
late time effective hydrodynamic expansion such that the time-dependent
energy-momentum tensor approaches that of the attractor for any arbitrary
initial condition Heller:2011ju ; Heinz:2013th ; Heller:2013fn ;
Heller:2015dha ; Heller:2016gbp ; Spalinski:2017mel ; Romatschke:2017vte ;
Romatschke:2017acs ; Denicol:2017lxn ; Behtash:2017wqg ; Casalderrey-
Solana:2017zyh ; Denicol:2018pak ; Strickland:2019hff ; Jaiswal:2019cju ;
Kurkela:2019set ; Giacalone:2019ldn ; Dash:2020zqx ; Heller:2020anv ;
Blaizot:2020gql ; Almaalol:2020rnu ; Ambrus:2021sjg ; Du:2021fok ;
Blaizot:2021cdv ; Chattopadhyay:2021ive ; Naik:2021yph ; Chen:2021wwh ;
Du:2022bel ; Almaalol:2022ijz ; Chen:2022ryi ; Alalawi:2022pmg , see
Soloviev:2021lhs  for a recent review. The hydrodynamic attractor in diverse
conformal systems, irrespective of whether they are weakly or strongly coupled
(e.g., a kinetic system or a holographic strongly interacting field theory),
are quasi-universal in the sense that the evolution is almost identical when
the phase space is described by suitable variables (to be introduced in due
course of time). The study of the hydrodynamic attractor has been extensively
carried out in the context of the Bjorken flow to learn about the
hydrodynamization of the quark-gluon plasma (QGP) produced in heavy-ion
collisions (see e.g. Berges:2020fwq ; Busza:2018rrf  for recent reviews),
although such a simplified set-up captures the longitudinal expansion only and
not the transverse flow.

To understand the hydrodynamization of the QGP, one needs to follow its
evolution from the earliest stages, where a perturbative description in the
forms of the glasma effective theory Gelis:2010nm  and subsequently kinetic
theory Arnold:2002zm ; Kurkela:2015qoa  is applicable, to later epochs when
the soft gluons form a strongly coupled thermal bath. The semi-holographic
framework, briefly reviewed in Sec. 2, allows for simple phenomenological
constructions where the evolution of both perturbative (weakly self-
interacting) and non-perturbative (strongly self-interacting) degrees of
freedom can be described simultaneously in a consistent way Faulkner:2010tq ;
Mukhopadhyay:2013dqa ; Iancu:2014ava ; Mukhopadhyay:2015smb ; Banerjee:2017ozx
; Kurkela:2018dku ; Ecker:2018ucc . In the simplest model, one retains only
two effective metric couplings, and reduces both sectors to viscous fluids
each following the Müller-Israel-Stewart (MIS) Muller:1967zza ; Israel:1979wp
description with different amounts of viscosities. The effective metric
couplings permit the exchange of energy and momentum, while the total energy-
momentum tensor, which can be constructed readily from the subsector energy-
momentum tensors, is conserved in the physical background (Minkowski) metric.
For further simplicity, the fluids describing both sectors are assumed to be
conformal. The full system deviates from conformality, however, since the
effective metric couplings (denoted by $\gamma$, $\tilde{\gamma}=r\gamma$)
have mass dimension $-4$. We assume them to be given in terms of a saturation
scale that sets the scale for the energy densities of the QGP when
perturbative descriptions tend to break down.111At this point, we note that
although we are calling the full non-perturbative description semi-
holographic, we are only using an effective fluid description for the strongly
coupled non-perturbative sector with the expectation that our simplification
should be sufficient to learn about hydrodynamization. The explicit dual
holographic description will be necessary to understand the role of
irreversible energy transfer to a dynamical black hole which is expected to be
a slow process determined by the coupling between the two sectors based on
results obtained in Ecker:2018ucc .

The existence of hydrodynamic attractors in such hybrid fluid models has been
established in Mitra:2020mei . These attractor curves form a two-dimensional
surface in four-dimensional phase space such that for any initial condition
the evolution in phase space approaches a curve on this attractor surface
after some time. At late time, the full system can be described as a single
fluid with equation of state and transport coefficients dependent on two
constant parameters which label the curves on the attractor surface. A robust
feature of the attractor surface is an initial evolution reminiscent of the
so-called bottom-up thermalization scenario of Ref. Baier:2000sb , for at
sufficiently early time, the energy density of the perturbative sector always
dominates over the strongly self-interacting non-perturbative sector.
Remarkably, the universality of this feature of the bottom-up thermalization
scenario thus emerges simply from just effective fluid descriptions within the
non-perturbative semi-holographic framework. The initial dominance of the
weakly coupled sector is determined by an exponential factor with an exponent
determined only by the viscosities and relaxation times of the two sectors.
Note, however, that our set-up certainly does not apply at arbitrarily early
times, where it will be at best a crude approximation to a full semi-
holographic treatment.

In this paper, we systematically study hydrodynamization in these hybrid
hydrodynamic attractors following preliminary observations in Mitra:2020mei
indicating that one can obtain significant new insights into small vs large
system collisions (high-multiplicity $p$-$p$ or $A$-$p$ vs. $A$-$A$ heavy ion
collisions Loizides:2016tew ; Schlichting:2016sqo  by detailed exploration of
the attractor surface.

The general characterization of the hydrodynamization times of the two
components on the hybrid attractor surface is as follows. Let us denote the
hydrodynamization times of the perturbative and non-perturbative (holographic)
sectors as $\tau_{hd}$ and $\widetilde{\tau}_{hd}$ respectively. When the
total energy density at the hydrodynamization times is in the regime

$0.0001\gamma^{-1}\lessapprox\mathcal{E}(\tau_{hd}),\mathcal{E}(\widetilde{\tau}_{hd})\lessapprox\gamma^{-1},$

we find that the hydrodynamization times of both sectors follow conformal
behaviors, i.e.

$\tau_{hd}\approx p\times(4\pi C_{\eta})\times\mathcal{E}(\tau_{hd})^{-1/4}$

and

$\widetilde{\tau}_{hd}\approx\widetilde{p}\times(4\pi\widetilde{C}_{\eta})\times\mathcal{E}(\widetilde{\tau}_{hd})^{-1/4}$

with $p$ and $\widetilde{p}$ being constants. Generically, due to the feature
of bottom-up thermalization discussed above, the energy density of the
perturbative sector is dominant at some reference time
$\tau_{0}\ll\gamma^{1/4}$. For such cases, we find that

$p\approx 0.63,$

which is approximately the same value as in usual conformal attractors (which
we obtain in the decoupling limit). Therefore, the hydrodynamization time of
the perturbative sector in the conformal regime is both universal and
identical to that of conformal attractors studied earlier in the literature.
Furthermore, $\widetilde{p}$ for the strongly self-interacting sector depends
only on the ratio of energy densities of the two subsectors at early time. As
this ratio of the energy densities of the perturbative to holographic sectors
at time (say) $\tau_{0}=0.001\gamma^{1/4}$ is increased from $1:1$ to
$1000:1$, the constant $\widetilde{p}$ increases from $0.64$ to $2.1$. The
hydrodynamization time of the relatively soft sector is thus sensitive to
initial conditions, and can be significantly larger than in usual conformal
attractors even in the conformal window.

Remarkably, for fine-tuned initial conditions on the attractor surface, where
the holographic sector has dominant energy at a fixed early initialization
time, the roles in the conformal window is reversed. The $\widetilde{p}$ of
the holographic sector is then universally $\approx 0.63$ as in the usual
conformal attractors, while $p$ of the perturbative sector is determined by
the ratio of the energy densities at the initialization time. We can thus
deduce that the hydrodynamization time of the sub-dominant sector increases
because it is forced to share its energy with the dominant sector prior to
hydrodynamization, while the dominant sector is dynamically driven to behave
in a universal manner.

On the attractor surface, the hydrodynamization times of both subsectors take
their minimal values when the respective energy densities are close to
$\gamma^{-1}$, i.e.

$\mathcal{E}({\tau}_{hd}),\mathcal{E}(\widetilde{\tau}_{hd})\sim\gamma^{-1}.$

For larger $\mathcal{E}({\tau}_{hd})$, $\mathcal{E}(\widetilde{\tau}_{hd})$
the hydrodynamization times of both sectors increase rapidly and do not follow
any scaling behavior like that in the conformal window.

For pondering potential phenomenological consequences of these general
features of the hybrid attractor, we interpret solutions with lower total
energy densities at the reference time $\tau\sim\gamma^{1/4}$ as corresponding
to small systems (high-multiplicity $p$-$p$ or $A$-$p$ collisions)
Loizides:2016tew ; Schlichting:2016sqo  as opposed to $A$-$A$ heavy-ion
collisions. To quantify our results, we identify the energy scale of the
intersystem coupling $\gamma^{-1/4}$ with the saturation scale
$Q_{s}^{p\mbox{-}p}$ of $p$-$p$ collisions, and set various initial conditions
such that the total energy density $\propto Q_{s}^{4}$ at $\tau\approx
Q_{s}^{-1}$, with $Q_{s}$ in $A$-$A$ systems chosen from typical values in the
IPsat and bCGC models Kowalski:2007rw . Generically, we find:

  * •

The total energy densities at hydrodynamization time of the perturbative
sector ($\mathcal{E}({\tau}_{hd})$) are within the conformal window
approximately for both $p$-$p$ and Pb-Pb collisions. Since
$\mathcal{E}({\tau}_{hd})$ is smaller in $p$-$p$ collisions, naturally the
hydrodynamization time of the perturbative sector is significantly larger in
$p$-$p$ collisions as readily follows from the scaling behavior of $\tau_{hd}$
in the conformal window.

  * •

The total energy densities at hydrodynamization time of the holographic sector
($\mathcal{E}(\widetilde{\tau}_{hd})$) are outside the conformal window for
both $p$-$p$ and Pb-Pb collisions, i.e.
$\mathcal{E}(\widetilde{\tau}_{hd})>\gamma^{-1}$. Since
$\mathcal{E}({\tau}_{hd})$ is smaller in $p$-$p$ collisions, the
hydrodynamization time of the softer sector is smaller in $p$-$p$ collisions
than in the case of Pb-Pb collisions. This follows from the growth of
$\widetilde{\tau}_{hd}$ with $\mathcal{E}(\widetilde{\tau}_{hd})$.

Our simple hybrid model therefore provides a new perspective on how collective
flow can emerge from the softer components even in small system collisions.

The organization of the paper is as follows. In Sec. 2, we review our set-up
of two fluids undergoing Bjorken flow and coupled mutually via the democratic
effective metric coupling Banerjee:2017ozx ; Kurkela:2018dku . We further
discuss the method to discern the attractor surface in this hybrid system. In
Sec. 3, we provide more analytic details of early time behavior on the
attractor surface, before turning our attention to the late-time hydrodynamic
behavior. In both cases, we expand the discussion beyond our earlier results
in Mitra:2020mei . Finally in Sec. 4, we systematically study the
hydrodynamization on the attractor surface and connect to the phenomenology of
$p$-$p$ and Pb-Pb collisions.

##  2 Set-up

###  2.1 The hybrid action

A minimalist set-up required to obtain our hybrid fluid model involves only
two leading democratic effective metric couplings Banerjee:2017ozx ;
Kurkela:2018dku . The action formalism for the full theory describing the
perturbative and non-perturbative degrees of freedom coupled via the mutual
democratic effective metric coupling is Kurkela:2018dku

$\displaystyle S[\psi,\tilde{\psi},G_{\mu\nu}]$
$\displaystyle=S[\psi,{g}_{\mu\nu}]+\tilde{S}[\tilde{\psi},{\tilde{g}}_{\mu\nu}]+S_{\rm
int}[{g}_{\mu\nu},{\tilde{g}}_{\mu\nu},G_{\mu\nu}],$ $\displaystyle S_{\rm
int}$ $\displaystyle=\frac{1}{2\gamma}\int
d^{d}x\sqrt{-G}(g_{\mu\alpha}-G_{\mu\alpha})G^{{\alpha\beta}}({\tilde{g}}_{\nu\beta}-G_{\nu\beta})G^{\mu\nu}$
$\displaystyle+\frac{1}{2\gamma}\frac{\gamma^{\prime}}{d\gamma^{\prime}-\gamma}\int
d^{d}x\sqrt{-G}(g_{\mu\nu}G^{\mu\nu}-d)({\tilde{g}}_{\alpha\beta}G^{\alpha\beta}-d).$
(1)

Above $S$ and $\tilde{S}$ denote the Wilsonian effective actions for the
perturbative and non-perturbative sectors respectively, with respective
degrees of freedom $\psi$ and $\tilde{\psi}$, and living in respective
effective metrics ${g}_{\mu\nu}$ and ${\tilde{g}}_{\mu\nu}$, while the
physical metric of the full theory is $G_{\mu\nu}$ (to be set to Minkowski
metric eventually). Both of these sectors inhabit the same topological space,
and they are superficially invisible to each other as they experience
different background metrics, but those are determined by the respective
opposite system. If $\tilde{S}$ is strongly self-interacting and admits a
holographic description, then it is the on-shell action of an extra-
dimensional gravitational theory whose boundary metric should be identified
with ${\tilde{g}}$. In what follows, we will not need any explicit description
of both $S$ and $\tilde{S}$.

Since both $S$ and $\tilde{S}$ are individually diffeomorphism-invariant in
the respective background metrics ${g}$ and ${\tilde{g}}$, we obtain the two
Ward identities

$\displaystyle\nabla_{\mu}t^{\mu\nu}=0\quad\text{and}\quad\tilde{\nabla}_{\mu}\tilde{t}^{\mu\nu}=0,$
(2)

where the energy momentum tensors are

$\displaystyle t^{\mu\nu}=\frac{2}{\sqrt{-g}}\frac{\delta S}{\delta
g_{\mu\nu}}\quad\text{and}\quad\tilde{t}^{\mu\nu}=\frac{2}{\sqrt{-{\tilde{g}}}}\frac{\delta\tilde{S}}{\delta{\tilde{g}}_{\mu\nu}},$
(3)

and $\nabla$ and $\tilde{\nabla}$ are built out of ${g}$ and ${\tilde{g}}$,
respectively.

The interaction term $S_{int}$ in (2.1) relates these two different background
metrics in terms of the actual physical background metric $G$ ultralocally, as
follows. We readily see that the two effective metrics appear in the full
action as auxiliary fields, and extremizing the action with respect to them
leads to imposing

$\frac{\delta S}{\delta{g}_{\mu\nu}}=0,\quad\frac{\delta
S}{\delta{g}_{\mu\nu}}=0.$ (4)

The above yield the following coupling equations

$\displaystyle{g}_{\mu\nu}$
$\displaystyle=G_{\mu\nu}+\frac{\sqrt{-{\tilde{g}}}}{\sqrt{-G}}\left[\gamma\,\tilde{t}^{\alpha\beta}G_{\alpha\mu}G_{\beta\nu}+\gamma^{\prime}\,t^{\alpha\beta}G_{\alpha\beta}G_{\mu\nu}\right],$
$\displaystyle{\tilde{g}}_{\mu\nu}$
$\displaystyle=G_{\mu\nu}+\frac{\sqrt{-{g}}}{\sqrt{-G}}\left[\gamma\,t^{\alpha\beta}G_{\alpha\mu}G_{\beta\nu}+\gamma^{\prime}\,t^{\alpha\beta}G_{\alpha\beta}G_{\mu\nu}\right],$
(5)

where $\gamma$ and $\gamma^{\prime}$ are the two dimensionful couplings with
mass dimension $-4$. It is expected that
$\gamma,\gamma^{\prime}\sim\mathcal{O}(\Lambda_{int}^{-4})$, where
$\Lambda_{int}$ is a suitable scale separating perturbative and non-
perturbative phenomena (and is also larger than the relevant scale of
observation). In the context of heavy-ion collisions, $\Lambda_{int}$ can be
identified with a saturation scale $Q_{s}$ determining the initial conditions.
It is evident from (2.1) that higher order effective metric couplings will be
suppressed by powers of $(T/Q_{s})^{4}$, where $T$ is an effective temperature
determining the scales of the stress tensors.

We can compute the energy momentum tensor of the full system by varying the
action (2.1) with respect to the physical background metric $G_{\mu\nu}$,
which yields

$\displaystyle T^{\mu}_{\phantom{\mu}\nu}$
$\displaystyle=\frac{1}{2}\left(\frac{\sqrt{-g}}{\sqrt{-G}}(t^{\mu}_{\phantom{\mu}\nu}+t_{\nu}^{\phantom{\nu}\mu})+\frac{\sqrt{-{\tilde{g}}}}{\sqrt{-G}}(\tilde{t}^{\mu}_{\phantom{\mu}\nu}+\tilde{t}_{\nu}^{\phantom{\nu}\mu})\right)+\Delta
K\delta^{\mu}_{\nu},$ $\displaystyle\Delta K$
$\displaystyle=-\frac{1}{2}\frac{\sqrt{-g}}{\sqrt{-G}}\frac{\sqrt{-{\tilde{g}}}}{\sqrt{-G}}\Big{[}\gamma\,t^{\mu\nu}G_{\mu\alpha}t^{\alpha\beta}G_{\beta\nu}+\gamma^{\prime}\,t^{\alpha\beta}G_{\alpha\beta}\tilde{t}^{\rho\sigma}G_{\rho\sigma}\Big{]}.$
(6)

In the first equation above, all indices for subsystem variables are lowered
using the respective background metrics, while that of the full system (on the
left hand side) is lowered using the physical background metric. It can be
explicitly checked that the coupling equations (2.1) and the two Ward
identities (2) indeed imply the local conservation of the full energy-momentum
tensor in the physical background metric, i.e.

$\displaystyle\nabla_{\mu}^{(B)}T^{\mu}_{\phantom{\mu}\nu}=0.$ (7)

The mutual effective metric coupling thus leads to a full energy-momentum
tensor (6) which can be obtained ultra-locally from the subsystem energy-
momentum tensors alone, and does not require any further microscopic input.
This ensures that a coarse-grained hydrodynamic description of the full
energy-momentum tensor in terms of those of the subsystems. It implements the
hypothetical Wilsonian perspective that the coarse-grained perturbative
description should determine the complete non-perturbative description at any
scale. Therefore, assuming that $\tilde{S}$ and the inter-system couplings can
be derived from $S$, we can state that the coarse-grained descriptions of the
two subsystems are related,222In Banerjee:2017ozx , a derivation of the semi-
holographic procedure was proposed by noting that the perturbative series has
renormalons which forbid Borel resummation, and this can be cured by
considering non-perturbative contributions. A physical observable, such as the
full energy-momentum tensor, which receives contributions from both sectors,
should be a function of the perturbative coupling with a finite radius of
convergence, provided that the intersystem couplings are appropriate functions
of the perturbative couplings. This was illustrated in a toy bi-holographic
set-up. and the full hydrodynamic description can be constructed in terms of
the subsystem hydrodynamic variables. A detailed description will follow soon.

In most set-ups, the full system is solved in a self-consistent iterative
procedure Mukhopadhyay:2015smb ; Ecker:2018ucc ; Mondkar:2021qsf . At the
first step of iteration, the effective metrics are identified with the
background metric $G_{\mu\nu}$, and the two systems are solved independently
with given initial data. In the next step, the effective metrics of the two
subsystems are updated via the coupling equations (2.1) using the respective
subsystem stress tensors obtained from the previous iteration. The subsystems
are solved once again with the same initial conditions. The iterations are
thus continued with fixed initial conditions until they converge implying that
the full energy-momentum tensor is conserved in the physical background
metric.

This iterative procedure was first proposed in the context of heavy-ion
collisions with glasma-like initial conditions for the perturbative sector and
an empty $AdS$-like initial conditions for the holographic sector
Iancu:2014ava . Subsequent works Mukhopadhyay:2015smb ; Ecker:2018ucc ;
Mondkar:2021qsf  have shown that the convergence of the semi-holographic
iterative procedure can not only be achieved, but to a very good accuracy over
three to four iterations only. For example, the convergence of the iterative
procedure has been demonstrated with non-trivial initial conditions for a
scalar field in the boundary coupled to the bulk Ecker:2018ucc ;
Mondkar:2021qsf . Note that since both systems are marginally deformed via
their effective metrics in each step of the iteration, the iterative procedure
is well defined in the full quantum theory, however we will not discuss this
further here.

Here, we will restrict ourselves to a large $N$ approximation in both sub-
sectors which ensures factorization of expectation values of multi-trace
operators into those of the single trace operators like the energy-momentum
tensor, and thus ignore stochastic hydrodynamic phenomena. Furthermore, we
will see that the iterative procedure is also not necessary in the
hydrodynamic limit – the coupling equations and the dynamical equations can be
evolved simultaneously ensuring conservation of the full energy-momentum
tensor in the physical background metric automatically.

###  2.2 Hybrid fluid model

Although the above discussion included an action, it is not necessary to
describe the evolution of the system, which is governed by the effective
conservation laws (2) and the algebraic coupling equations (2.1) in the
hydrodynamic limit. The only required inputs are the energy-momentum tensors
of each sector, $t^{\mu\nu}$ and $\tilde{t}^{\mu\nu}$, determined via the
constitutive relations, and the physical background metric $G_{\mu\nu}$. Since
the hydrodynamic limit is found in many theories (we have in mind a strongly
coupled holographic theory for one subsystem and a weakly coupled kinetic
theory for the other), we can consider the simpler case in which the two
subsystems are two different relativistic fluids.

Hence, consider two conformal fluids with energy momentum tensors in their
respective metrics, given by

$\displaystyle t^{\mu\nu}$
$\displaystyle=(\varepsilon+P)u^{\mu}u^{\nu}+P{g}^{\mu\nu}+\Pi^{\mu\nu},$ (8)
$\displaystyle\tilde{t}^{\mu\nu}$
$\displaystyle=(\tilde{\varepsilon}+\tilde{P})\tilde{u}^{\mu}\tilde{u}^{\nu}+\tilde{P}{\tilde{g}}^{\mu\nu}+\tilde{\Pi}^{\mu\nu},$
(9)

where $t^{\mu\nu}{g}_{\mu\nu}=-\varepsilon+3P=0$,
$\tilde{t}^{\mu\nu}{\tilde{g}}_{\mu\nu}=-\tilde{\varepsilon}+3\tilde{P}=0$.
The energy momentum tensors are conserved in their respective backgrounds (2).
The four-velocities are properly normalized with respect to their metrics,
i.e.

$\displaystyle
u^{\mu}{g}_{\mu\nu}u^{\nu}=-1\quad\text{and}\quad\tilde{u}^{\mu}{\tilde{g}}_{\mu\nu}\tilde{u}^{\nu}=-1.$
(10)

The dissipation tensors, $\Pi^{\mu\nu}$ and $\tilde{\Pi}^{\mu\nu}$, are taken
to be orthogonal to their respective four velocities

$\displaystyle\Pi^{\mu\nu}u_{\mu}=0\quad\text{and}\quad\tilde{\Pi}^{\mu\nu}\tilde{u}_{\mu}=0.$
(11)

The two fluids are assumed to have dissipative dynamics dictated by Müller-
Israel-Stewart (MIS) theory Muller:1967zza ; Israel:1979wp ,333The more weakly
and the more strongly coupled sectors are also simply assumed to have higher
and lower specific viscosity, respectively, as well as correspondingly shorter
and longer relaxation time scales. While not pursued here, one could also
consider including second-order derivatives of the shear tensor of the
strongly coupled sector as in Noronha:2011fi ; Heller:2014wfa  such as to make
the dynamics of the latter more similar to the holographic dynamics, which is
governed by quasinormal modes. (We thank Michal Heller for this remark.)
namely

$\displaystyle(\tau_{\pi}u^{\mu}\nabla_{\mu}+1)\Pi^{\mu\nu}$
$\displaystyle=-\eta\sigma^{\mu\nu},$
$\displaystyle(\tilde{\tau}_{\pi}\tilde{u}^{\mu}\tilde{\nabla}_{\mu}+1)\tilde{\Pi}^{\mu\nu}$
$\displaystyle=-\tilde{\eta}\tilde{\sigma}^{\mu\nu},$ (12)

where $\eta$ and $\tilde{\eta}$ are the subsystem viscosities. Note that we
have truncated at first order, and since we are working with conformal
subsystems, we have set the bulk viscosities to zero.

It is worthwhile to underline that although the two subsystems are conformal,
the complete system (6) will in general not be conformal,
$T^{\mu\nu}g^{(B)}_{\mu\nu}\neq 0$. We will see that there is an emergent
conformality for large proper time. It is worthwhile to mention that the full
system also has emergent conformality at thermal equilibrium at high
temperature (with respect to the scale of the intersystem coupling)
Kurkela:2018dku .

###  2.3 Explicit equations for boost-invariant Bjorken flow

We now provide the explicit form of the complete set of equations, (2), (2.1)
and (12). We will work in the Milne coordinates: $\tau=\sqrt{t^{2}-z^{2}}$,
$x$, $y$, and $\zeta={\rm arctanh}(z/t)$, where $x$ and $y$ are the transverse
coordinates and $z$ is the longitudinal coordinate for the expanding system.
In these coordinates, the Minkowski metric assumes the form

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

The above is thus the physical background metric for the full system. The
Bjorken flow is boost-invariant, and therefore in the Milne coordinates, all
physical variables depend only on the proper time $\tau$.

Motivated by the form of the background metric, we make the following boost-
invariant ansatz for the effective metrics:

$\displaystyle{g}_{\mu\nu}(\tau)=\text{diag}(-a(\tau)^{2},b(\tau)^{2},b(\tau)^{2},c(\tau)^{2}),$
$\displaystyle{\tilde{g}}_{\mu\nu}(\tau)=\text{diag}(-\tilde{a}(\tau)^{2},\tilde{b}(\tau)^{2},\tilde{b}(\tau)^{2},\tilde{c}(\tau)^{2}).$
(14)

The four velocities are then

$\displaystyle
u^{\mu}=\left(\frac{1}{a(\tau)},0,0,0\right)\quad\text{and}\quad
u^{\mu}=\left(\frac{1}{\tilde{a}(\tau)},0,0,0\right).$ (15)

Here we consider the dissipation tensors, $\Pi^{\mu\nu}$ and
$\tilde{\Pi}^{\mu\nu}$, to be symmetric, transverse and traceless i.e.,
${g}_{\mu\nu}\Pi^{\mu\nu}=0$ and ${\tilde{g}}_{\mu\nu}\tilde{\Pi}^{\mu\nu}$.
These conditions imply that there is only one independent component of
$\Pi^{\mu\nu}$ for each sector, which we call $\Pi^{\eta}_{\,\,\eta}=-\phi$
(the longitudinal non-equilibrium pressure) following the convention in
Heller:2015dha , which means that

$\displaystyle\Pi^{\mu\nu}$
$\displaystyle=\text{diag}\left(0,\frac{\phi(\tau)}{2},\frac{\phi(\tau)}{2},\frac{-\phi(\tau)}{\tau^{2}}\right),$
(16) $\displaystyle\tilde{\Pi}^{\mu\nu}$
$\displaystyle=\text{diag}\left(0,\frac{\tilde{\phi}(\tau)}{2},\frac{\tilde{\phi}(\tau)}{2},\frac{-\tilde{\phi}(\tau)}{\tau^{2}}\right).$
(17)

Then explicitly, the only non-vanishing components of the conservation
equations (2) are

$\displaystyle 0$
$\displaystyle=\partial_{\tau}\varepsilon+\frac{4}{3}\varepsilon\partial_{\tau}\log(b^{2}c)+\phi\partial_{\tau}\log(\frac{b}{c}),$
(18) $\displaystyle 0$
$\displaystyle=\partial_{\tau}\tilde{\varepsilon}+\frac{4}{3}\tilde{\varepsilon}\partial_{\tau}\log(\tilde{b}^{2}\tilde{c})+\tilde{\phi}\partial_{\tau}\log(\frac{\tilde{b}}{\tilde{c}}),$
(19)

the MIS equations are (12)

$\displaystyle 0$
$\displaystyle=\tau_{\pi}\partial_{\tau}\phi+\frac{4}{3}\eta\partial_{\tau}\log\frac{b}{c}+\left[a+\frac{4}{3}\tau_{\pi}\log(b^{2}c)\right]\phi,$
(20) $\displaystyle 0$
$\displaystyle=\tilde{\tau}_{\pi}\partial_{\tau}\tilde{\phi}+\frac{4}{3}\tilde{\eta}\partial_{\tau}\log\frac{\tilde{b}}{\tilde{c}}+\left[\tilde{a}+\frac{4}{3}\tilde{\tau}_{\pi}\log(\tilde{b}^{2}\tilde{c})\right]\tilde{\phi},$
(21)

and finally, the algebraic coupling equations (14) are

$\displaystyle 1-a^{2}$
$\displaystyle=\gamma\frac{\tilde{a}\tilde{b}^{2}\tilde{c}}{\tau}\Big{[}\frac{\tilde{\varepsilon}}{\tilde{a}^{2}}+r(-\frac{\tilde{\varepsilon}}{\tilde{a}^{2}}+\frac{\tilde{\phi}}{\tilde{b}^{2}}+\frac{2\tilde{\varepsilon}}{3\tilde{b}^{2}}+\tau^{2}(\frac{\tilde{\varepsilon}}{3\tilde{c}^{2}}-\frac{\tilde{\phi}}{\tilde{c}^{2}}))\Big{]},$
(22) $\displaystyle b^{2}-1$
$\displaystyle=\gamma\frac{\tilde{a}\tilde{b}^{2}\tilde{c}}{\tau}\Big{[}\frac{(\tilde{\phi}+\frac{2\tilde{\varepsilon}}{3})}{2\tilde{b}^{2}}-r(-\frac{\tilde{\varepsilon}}{\tilde{a}^{2}}+\frac{\tilde{\phi}}{\tilde{b}^{2}}+\frac{2\tilde{\varepsilon}}{3\tilde{b}^{2}}+\tau^{2}(\frac{\tilde{\varepsilon}}{3\tilde{c}^{2}}-\frac{\tilde{\phi}}{\tilde{c}^{2}}))\Big{]},$
(23) $\displaystyle c^{2}-\tau^{2}$
$\displaystyle=\gamma\frac{\tilde{a}\tilde{b}^{2}\tilde{c}}{\tau}\Big{[}\frac{\tau^{4}}{\tilde{c}^{2}}(\frac{\tilde{\varepsilon}}{3}-\tilde{\phi})-r\tau^{2}(-\frac{\tilde{\varepsilon}}{\tilde{a}^{2}}+\frac{\tilde{\phi}}{\tilde{b}^{2}}+\frac{2\tilde{\varepsilon}}{3\tilde{b}^{2}}+\tau^{2}(\frac{\tilde{\varepsilon}}{3\tilde{c}^{2}}-\frac{\tilde{\phi}}{\tilde{c}^{2}}))\Big{]},$
(24) $\displaystyle 1-\tilde{a}^{2}$
$\displaystyle=\gamma\frac{ab^{2}c}{\tau}\Big{[}\frac{\varepsilon}{a^{2}}+r(-\frac{\varepsilon}{a^{2}}+\frac{\phi}{b^{2}}+\frac{2\varepsilon}{3b^{2}}+\tau^{2}(\frac{\varepsilon}{3c^{2}}-\frac{\phi}{c^{2}}))\Big{]},$
(25) $\displaystyle\tilde{b}^{2}-1$
$\displaystyle=\gamma\frac{ab^{2}c}{\tau}\Big{[}\frac{(\phi+\frac{2\varepsilon}{3})}{2b^{2}}-r(-\frac{\varepsilon}{a^{2}}+\frac{\phi}{b^{2}}+\frac{2\varepsilon}{3b^{2}}+\tau^{2}(\frac{\varepsilon}{3c^{2}}-\frac{\phi}{c^{2}}))\Big{]},$
(26) $\displaystyle\tilde{c}^{2}-\tau^{2}$
$\displaystyle=\gamma\frac{ab^{2}c}{\tau}\Big{[}\frac{\tau^{4}}{c^{2}}(\frac{\varepsilon}{3}-\phi)-r\tau^{2}(-\frac{\epsilon(\tau)}{a^{2}}+\frac{\phi}{b^{2}}+\frac{2\varepsilon}{3b^{2}}+\tau^{2}(\frac{\varepsilon}{3c^{2}}-\frac{\phi}{c^{2}}))\Big{]},$
(27)

where we introduced $r\equiv-\gamma^{\prime}/\gamma$. We see that in the
decoupling limit, namely $\gamma=0$, the right hand side vanishes and the
effective metrics reduce to the usual background Milne metric (13) in which
the two fluids evolve independently.

It is clear that the above system actually has four dynamical variables,
namely the energy densities and the non-equilibrium pressures
$(\varepsilon,\tilde{\varepsilon},\phi,\tilde{\phi})$ which follow the first-
order equations (18), (19), (20) and (21). The six effective metric variables,
$a,\tilde{a},b,\tilde{b},c,\tilde{c}$, can be eliminated by determining them
as functions of the four dynamical variables by solving the six algebraic
equations (22)-(27). The evolution of the full system can thus be readily
determined numerically for given initial values of
$(\varepsilon,\tilde{\varepsilon},\phi,\tilde{\phi})$.

In what follows, we parametrize the transport coefficients in the two
conformal subsystems, and the respective MIS relaxation times via
dimensionless quantities $C_{\eta}$, $\tilde{C}_{\eta}$, $C_{\tau}$ and
$\tilde{C}_{\tau}$ defined via

$\displaystyle\eta=C_{\eta}\frac{\varepsilon+P}{\varepsilon^{1/4}},\quad$
$\displaystyle\text{and}\quad\tau_{\pi}=C_{\tau}\varepsilon^{-1/4},$ (28)
$\displaystyle\tilde{\eta}=\tilde{C}_{\eta}\frac{\tilde{\varepsilon}+\tilde{P}}{\tilde{\varepsilon}^{1/4}},\quad$
$\displaystyle\text{and}\quad\tilde{\tau}_{\pi}=\tilde{C}_{\tau}\tilde{\varepsilon}^{-1/4}.$
(29)

In the simulations that follow, we will usually choose

$\displaystyle\tilde{C}_{\eta}=\frac{1}{4\pi},\quad\tilde{C}_{\tau}=\frac{2-\log
2}{2\pi},\quad C_{\eta}=10\tilde{C}_{\eta},\quad\text{and}\quad
C_{\tau}=5C_{\eta}.$ (30)

motivated by the values obtained in holographic descriptions and in kinetic
theories, respectively. This marks the tilded sector as more strongly coupled
than the untilded one.

We shall set

$\displaystyle r\equiv-\gamma^{\prime}/\gamma=2$ (31)

unless noted otherwise. In the inviscid Minkowski case Kurkela:2018dku , the
condition of UV completeness restricts the value of $r$ to $r>1$. Moreover, it
was found in the inviscid case that there is a second order phase transition
for $1<r\lessapprox 1.11$, which we discuss in Appendix B. Note that the
coupling $\gamma$ is the only dimensionful scale in our framework. Further, we
will work with the dimensionless anisotropy parameter,
$\chi=\phi/(\epsilon+P)$. Thus the conservation and the MIS equations read

$\displaystyle\frac{4}{3}\chi\left(\frac{b^{\prime}}{b}-\frac{c^{\prime}}{c}\right)+\frac{4}{3}\left(\frac{2b^{\prime}}{b}+\frac{c^{\prime}}{c}\right)+\frac{\varepsilon^{\prime}}{\varepsilon}$
$\displaystyle=0,$ (32)
$\displaystyle\frac{4}{3}\tilde{\chi}\left(\frac{\tilde{b}^{\prime}}{\tilde{b}}-\frac{\tilde{c}^{\prime}}{\tilde{c}}\right)+\frac{4}{3}\left(\frac{2\tilde{b}^{\prime}}{\tilde{b}}+\frac{\tilde{c}^{\prime}}{\tilde{c}}\right)+\frac{\tilde{\varepsilon}^{\prime}}{\tilde{\varepsilon}}$
$\displaystyle=0,$ (33)
$\displaystyle-\frac{a\chi{\varepsilon}^{1/4}}{C_{\text{$\tau$}}}-\Big{(}\frac{4C_{\text{$\eta$}}}{3C_{\text{$\tau$}}}-\frac{4\chi{}^{2}}{3}\Big{)}\Big{(}\frac{b^{\prime}}{b}-\frac{c^{\prime}}{c}\Big{)}-\chi^{\prime}$
$\displaystyle=0,$ (34)
$\displaystyle-\frac{\tilde{a}\tilde{\chi}{\tilde{\varepsilon}}^{1/4}}{\tilde{C}_{\tau}}-\Big{(}\frac{4\tilde{C}_{\eta}}{3\tilde{C}_{\text{$\tau$}}}-\frac{4\tilde{\chi}{}^{2}}{3}\Big{)}\Big{(}\frac{\tilde{b}^{\prime}}{\tilde{b}}-\frac{\tilde{c}^{\prime}}{\tilde{c}}\Big{)}-\tilde{\chi}^{\prime}$
$\displaystyle=0,$ (35)

where the prime denotes a derivative with respect to the proper time, e.g.
$\varepsilon^{\prime}=\partial_{\tau}\varepsilon.$ The full energy momentum
tensor (6) is given by

$\displaystyle\mathcal{T}^{\mu\nu}=\text{diag}(\mathcal{E},P_{T},P_{T},P_{L}),$
(36)

where the total energy density, total transverse pressure and longitudinal
pressure are

$\displaystyle\mathcal{E}=-\frac{ab^{2}c\varepsilon+\tilde{a}\tilde{b}^{2}\tilde{c}\tilde{\varepsilon}}{\tau}+\frac{\gamma}{36a\tilde{a}c\tilde{c}\tau^{2}}\Big{(}a^{2}\tilde{a}^{2}\Big{[}2b^{2}\tau^{2}(\varepsilon-3\phi)\left(\tilde{b}^{2}(r-1)\tau^{2}(\tilde{\varepsilon}-3\tilde{\phi})+\tilde{c}^{2}r(2\tilde{\varepsilon}+3\tilde{\phi})\right)$
$\displaystyle+c^{2}(2\varepsilon+3\phi)\left(2\tilde{b}^{2}r\tau^{2}(\tilde{\varepsilon}-3\tilde{\phi})+\tilde{c}^{2}(2r-1)(2\tilde{\varepsilon}+3\tilde{\phi})\right)\Big{]}$
$\displaystyle-6a^{2}\tilde{b}^{2}\tilde{c}^{2}r\tilde{\varepsilon}\left(b^{2}\tau^{2}(\epsilon-3\phi)+c^{2}(2\varepsilon+3\phi)\right)$
$\displaystyle-6b^{2}c^{2}\varepsilon\left(\tilde{b}^{2}\left(\tilde{a}^{2}r\tau^{2}(\tilde{\varepsilon}-3\tilde{\phi})-3\tilde{c}^{2}(r-1)\tilde{\varepsilon}\right)+\tilde{a}^{2}\tilde{c}^{2}r(2\tilde{\varepsilon}+3\tilde{\phi})\right)\Big{)},$
(37) $\displaystyle
P_{T}=\frac{ab^{2}c(2\varepsilon+3\phi)+\tilde{a}\tilde{b}^{2}\tilde{c}(2\tilde{\varepsilon}+3\tilde{\phi})}{6\tau}$
$\displaystyle+\frac{\gamma}{36ac\tilde{a}\tilde{c}\tau^{2}}\Big{[}ca\tilde{a}^{2}\Big{(}2b^{2}\tau^{2}(\varepsilon-3\phi)\left(\tilde{b}^{2}(r-1)\tau^{2}(\tilde{\varepsilon}-3\tilde{\phi})+\tilde{c}^{2}r(2\tilde{\varepsilon}+3\tilde{\phi})\right)$
$\displaystyle+c^{3}(2\varepsilon+3\phi)\left(2\tilde{b}^{2}r\tau^{2}(\tilde{\varepsilon}-3\tilde{\phi})+\tilde{c}^{2}(2r-1)(2\tilde{\varepsilon}+3\tilde{\phi})\right)\Big{)}$
$\displaystyle-6b^{2}c^{2}\varepsilon\left(\tilde{b}^{2}\left(\tilde{a}^{2}r\tau^{2}(\tilde{\varepsilon}-3\text{$\phi$t})-3\tilde{c}^{2}(r-1)\tilde{\varepsilon}\right)+\tilde{a}^{2}\tilde{c}^{2}r(2\tilde{\varepsilon}+3\tilde{\phi})\right)$
$\displaystyle-6ac\tilde{b}^{2}\tilde{c}^{2}r\tilde{\varepsilon}\left(b^{2}\tau^{2}(\varepsilon-3\phi)+c^{2}(2\epsilon+3\phi)\right)\Big{]},$
(38) $\displaystyle
P_{L}=P_{T}-\frac{3\left(ab^{2}c\phi+\tilde{a}\tilde{b}^{2}\tilde{c}\tilde{\phi}\right)}{2\tau}.$
(39)

###  2.4 Determining the attractor surface

As discussed above, the hybrid fluid system has a four-dimensional phase space
with the energy densities and non-equilibrium pressures,
$(\varepsilon,\tilde{\varepsilon},\phi,\tilde{\phi})$, as the dynamical
variables. It was shown in Mitra:2020mei  that the full system has a two-
dimensional attractor surface. For any initial condition
$(\varepsilon,\tilde{\varepsilon},\phi,\tilde{\phi})(\tau_{0})$, the system
evolves towards a particular curve on the attractor surface. The curves which
compose this attractor surface are characterized by dynamical evolutions where
the dimensionless anisotropies, $\chi$ and $\tilde{\chi}$ do not diverge as we
trace back in time to the moment of collision $\tau=0$. Furthermore,

$\lim_{\tau\rightarrow 0}\chi=\sigma_{1},\quad\lim_{\tau\rightarrow
0}\tilde{\chi}=\sigma_{2},$ (40)

where

$\displaystyle\sigma_{1}=\sqrt{\frac{C_{\eta}}{C_{\tau}}},\quad\sigma_{2}=\sqrt{\frac{\tilde{C}_{\eta}}{\tilde{C}_{\tau}}}.$
(41)

See Denicol:2017lxn  for similar analytic behavior of the attractor. In Fig. 1
we show the evolution of the anisotropies in the different sectors and the
total system for one particular attractor solution and neighboring
trajectories (from Mitra:2020mei , where more plots are given for this
particular solution).

Figure 1: Anisotropies $\chi$ of one particular attractor solution (thick
lines) and four neighboring trajectories (thin lines). The less viscous
(strongly coupled) system is displayed by blue curves, the more viscous (more
weakly coupled) system by red curves, the total system by orange curves. From
Mitra:2020mei , where more plots are given for this particular solution.

In order to determine the attractor curves on the attractor surface we proceed
as follows. We initialize at a time $\tau_{0}$ and choose our initial energy
densities and non-equilibrium pressures,
$(\varepsilon,\tilde{\varepsilon},\phi,\tilde{\phi})(\tau_{0})$. We then solve
the coupling equations (22)-(27) at the initial time to determine the metric
components at the initial time. Subsequently, we evolve the system backwards
in time via (18), (19), (20) and (21), and tune the values of
$(\phi,\tilde{\phi})(\tau_{0})$ with fixed
$(\varepsilon,\tilde{\varepsilon})(\tau_{0})$ such that $\chi$ and
$\tilde{\chi}$ tend to constants for early times following (40). Having thus
found an attractor curve to a satisfying level of precision (in this paper,
$10^{-15})$, we then evolve the system forward in time to study the late time
approach to hydrodynamics.

With parameters chosen as in (30) and (31), we have $\sigma_{1}=\sqrt{1/5}$
and $\sigma_{2}=\sqrt{1/(2(2-\ln 2))}$. We will initialize at a time
$\gamma^{-1/4}\tau_{0}=1$ and evolve backwards in time to discern the early
time behavior. Now let us find the rigid parameter by comparing different
attractor solutions.

For illustrative purposes, we will discuss one approximate attractor solution,
which can be found by choosing

$\displaystyle\gamma\varepsilon(\tau_{0})$
$\displaystyle=\gamma\tilde{\varepsilon}(\tau_{0})=0.18,$
$\displaystyle\gamma\phi_{1}(\tau_{0})$
$\displaystyle=0.10002105\quad\text{and}\quad\gamma\phi_{2}(\tau_{0})=0.0501539602.$
(42)

and we evolve the system backward in time.

Figure 2: Illustrative behaviour of an attractor curve – see text in Sec. 2.4
for the choice of parameters. Top: $\chi$ and $\tilde{\chi}$ approach
$\sigma_{1}$ and $\sigma_{2}$ respectively as $\tau\rightarrow 0$. Bottom:
$\varepsilon$ and $\tilde{\varepsilon}$ satisfy the scaling laws (3.1) at
early time. We also see that $\gamma\varepsilon$ tends towards
$\sqrt{(r-1)/r}=1/\sqrt{2}$ (orange dotted line) as claimed in (3.1).

We see that in the left panel of Fig. 2 that with these initial conditions,
$\chi\rightarrow\sigma_{1}$ and $\tilde{\chi}\rightarrow\sigma_{2}$ for early
times indicating that we are close to the attractor. Fig. 2 further shows that
the scalings (3.1) for $\varepsilon$ and $\tilde{\varepsilon}$ (see below)
come out to be as predicted in the limit $\tau\rightarrow 0$.

##  3 Mapping the attractor surface at early and late times

###  3.1 Early time behavior: dominance of the weakly coupled sector

The attractor surface has a universal scaling behavior at $\tau=0$ which is
determined in terms of these two dimensionless viscous parameters,
$\sigma_{1}$ and $\sigma_{2}$ defined in (41).

Crucially, we need

$0<\sigma_{1}<1\quad{\rm
and}\quad\frac{3}{4}\sigma_{1}+\frac{1}{4}<\sigma_{2}<1.$ (43)

The relations $\sigma_{1,2}<1$ are imposed by causality Romatschke:2017ejr .
The two relations above imply that $\sigma_{2}>\sigma_{1}$ as should be the
case if the second subsystem is more strongly coupled than the first. The
stricter lower bound on $\sigma_{2}$ is necessary for reasons to be described
below.

We first note that with the above conditions, the anisotropies, $\chi$ and
$\tilde{\chi}$, go to $\sigma_{1}$ and $\sigma_{2}$ respectively at $\tau=0$
on the attractor surface, as we see in Fig. 2. Moreover, we find numerically
that $\tilde{\varepsilon}$ goes to zero while $\varepsilon$ goes to a
constant, namely $\gamma^{-1}\sqrt{r-1}/\sqrt{r}$. With these numerical
inputs, assuming (43), and by expanding the conservation, MIS and coupling
equations in $\tau$ around $\tau=0$, we readily determine the following
universal early time scaling behavior of the energy densities and the
anisotropies:

$\displaystyle\gamma\varepsilon$
$\displaystyle\rightarrow\sqrt{\frac{r-1}{r}},\quad\quad\gamma\tilde{\varepsilon}\rightarrow{{g_{2}}}\tau^{\frac{4}{3}(\sigma_{2}-\sigma_{1})},$
$\displaystyle\chi$
$\displaystyle\rightarrow\sigma_{1}+k_{1}\tau^{\frac{1}{3}(\sigma_{1}+2)},\quad\quad\tilde{\chi}\rightarrow\sigma_{2}+k_{2}\tau^{\frac{5}{3}\sigma_{2}-\frac{4}{3}\sigma_{1}+\frac{2}{3}},$
(44)

and the following scaling behaviors of the components of the effective
metrics:

$\displaystyle a$ $\displaystyle\rightarrow
a_{10}\tau^{\frac{1}{3}(\sigma_{1}-1)},\quad b\rightarrow
b_{10}\tau^{\frac{1}{3}(\sigma_{1}-1)},\quad c\rightarrow
c_{10}\tau^{\frac{1}{3}(\sigma_{1}+2)},$ $\displaystyle\tilde{a}$
$\displaystyle\rightarrow
a_{20}\tau^{\frac{4}{3}\sigma_{2}-\sigma_{1}-\frac{1}{3}},\quad\tilde{b}\rightarrow
b_{20}\tau^{\frac{1}{3}(\sigma_{1}-1)},\quad\tilde{c}\rightarrow
c_{20}\tau^{\frac{1}{3}(\sigma_{1}+2)}.$ (45)

Since the attractor surface is two-dimensional, there should be only two
independent coefficients of the above early time expansions, which can be
chosen to be $k_{2}$ and $g_{2}$. The latter determine the remaining
coefficients as follows:

$\displaystyle
b_{10}=k_{2}g_{2}^{-5/4}\gamma^{{1}/{4}}\tilde{C}_{\tau}\frac{\sqrt{3}(13\sigma_{2}-4\sigma_{1}+2)\sqrt{r}}{\sigma_{2}(1-4r+2(1-r)\sigma_{1})\sqrt{1-4r-4(1-r)\sigma_{1}}},$
$\displaystyle a_{10}=\sqrt{\frac{r-1}{r}}b_{10},$ $\displaystyle
c_{10}=b_{10},$ $\displaystyle
a_{20}=-k_{2}g_{2}^{-1/4}\gamma^{{1}/{4}}\tilde{C}_{\tau}\frac{13\sigma_{2}-4\sigma_{1}+2}{3\sigma_{2}},$
$\displaystyle b_{20}=\sqrt{\frac{4r-1+2(r-1)\sigma_{1}}{3{r(r-1)}}}b_{10},$
$\displaystyle c_{20}=\sqrt{\frac{4r-1-4(r-1)\sigma_{1}}{3{r(r-1)}}}b_{10},$
$\displaystyle
k_{1}=-C_{\tau}^{-1}\gamma^{-\frac{1}{4}}r^{-5/8}(r-1)^{5/8}\frac{3\sigma_{1}}{2+9\sigma_{1}}b_{10}.$
(46)

As a result, in the limit $\tau\rightarrow 0$, the following five identities
should hold

$\displaystyle\frac{b}{\tau^{\frac{1}{3}(\sigma_{1}-1)}}=\frac{c}{\tau^{\frac{1}{3}(\sigma_{1}+2)}}=\sqrt{\frac{r}{r-1}}\frac{a}{\tau^{\frac{1}{3}(\sigma_{1}-1)}}=\frac{\tilde{b}}{\tau^{\frac{1}{3}(\sigma_{1}-1)}}\sqrt{r\gamma\tilde{\varepsilon}\frac{\tilde{c}}{\tau\tilde{a}}}$
(47)
$\displaystyle=\sqrt{\frac{3\sqrt{r(r-1)}}{\gamma\varepsilon(4r-1-2(r-1)\chi)}}\frac{\tilde{b}}{\tau^{\frac{1}{3}(\sigma_{1}-1)}}=\sqrt{\frac{3\sqrt{r(r-1)}}{\gamma\varepsilon(4r-1+4(r-1)\chi)}}\frac{\tilde{c}}{\tau^{\frac{1}{3}(\sigma_{1}+2)}}$

on any curve on the attractor surface. Fig. 3 demonstrates the numerical check
of the above identities on an attractor curve.

The stricter lower bound on $\sigma_{2}$ in Eq. (43) is necessary because the
above scaling relations follow from solving the coupling equations, etc., with
the assumption that $\tilde{a}$ vanishes as $\tau\rightarrow 0$. From the
explicit scaling of $\tilde{a}$ in (3.1), it is clear that the lower bound is
then self-consistent. Furthermore, it is also clear from (3.1) that
$\tilde{\chi}$ also does not diverge as

$\sigma_{2}-\frac{4}{5}\sigma_{1}+\frac{2}{5}>0$

as a consequence of the lower bound on $\sigma_{2}$ and $\sigma_{1}<1$.

If instead of Eq. (43) we impose

$0<\sigma_{2}<1\quad{\rm and}\quad{\rm
max}\left(0,\frac{3}{4}\sigma_{2}-\frac{1}{2}\right)<\sigma_{1}<\sigma_{2},$
(48)

we obtain that $\tilde{a}$ diverges as
$\tau^{(5/3)\sigma_{1}-(4/3)\sigma_{2}-(1/3)}$ as $\tau\rightarrow 0$ and all
other scaling relations discussed above are also altered. Nevertheless, $\chi$
and $\tilde{\chi}$ still reaches $\sigma_{1}$ and $\sigma_{2}$ as
$\tau\rightarrow 0$ on the attractor surface. However, the above condition
involves a lower bound on $\sigma_{1}$ which cannot be motivated otherwise,
while (43) can be readily satisfied with reasonable weak coupling and strong
coupling parameters. Therefore, we do not further investigate the above
condition.

Figure 3:  Check of the identities (47) at early times. Similar identities
given in (141) also holds for the disperser surface, and can be verified to a
very high accuracy numerically as evident from Fig. 19.

It is easy to see from the above that the physical energy densities of the two
subsystems behave as follows at early time

$\mathcal{E}_{1}=\varepsilon\frac{ab^{2}c}{\tau}\sim\tau^{\frac{4}{3}\sigma_{1}-\frac{4}{3}},\quad\mathcal{E}_{2}=\varepsilon\frac{\tilde{a}\tilde{b}^{2}\tilde{c}}{\tau}\sim\tau^{\frac{8}{3}\sigma_{2}-\frac{4}{3}\sigma_{1}-\frac{4}{3}},$
(49)

and therefore

$\frac{\mathcal{E}_{2}}{\mathcal{E}_{1}}\sim\tau^{\frac{8}{3}(\sigma_{2}-\sigma_{1})}.$
(50)

Since $\sigma_{2}>\sigma_{1}$, it follows that the energy density of weakly
self-interacting sector is always dominant at sufficiently early times,
leading to the universal feature of an energy distribution akin to the so-
called bottom-up thermalization scenario Baier:2000sb  on the attractor
surface. Also note that $\mathcal{E}_{1}$ always diverges as $\tau\rightarrow
0$ since $\sigma_{1}<1$. However, $\mathcal{E}_{2}$ may or may not diverge in
this limit, without contradicting any of the two conditions in (43), depending
on the choice of parameters.

While all initial conditions evolve towards specific curves on the attractor
surface, they also get repelled from a disperser surface in which $\chi$ and
$\tilde{\chi}$ approach $-\sigma_{1}$ and $-\sigma_{2}$ respectively in the
limit $\tau\rightarrow 0$. On the disperser surface, the opposite of the
energy distribution in bottom-up thermalization takes place, since both
$\mathcal{E}_{2}$ and $\mathcal{E}_{2}/\mathcal{E}_{1}$ diverge as
$\tau\rightarrow 0$. In contrast to (3.1), $\gamma\tilde{\varepsilon}$ goes to
the constant $\sqrt{r-1}/\sqrt{r}$ while $\gamma\varepsilon$ vanishes as
$\tau\rightarrow 0$ on the disperser surface. For more details, see Appendix
C. One can also find two other such half-disperser surfaces, in which $\chi$
and $\tilde{\chi}$ approach $\sigma_{1}$ and $-\sigma_{2}$ respectively, or
$\chi$ and $\tilde{\chi}$ approach $-\sigma_{1}$ and $\sigma_{2}$ respectively
as $\tau\rightarrow 0$. These can be analyzed in a similar manner. All these
three surfaces are not relevant for understanding hydrodynamization with
generic initial conditions.

###  3.2 Initialization and bounds on initial conditions

We initialize the evolution at $\tau_{0}>0$ on the two-dimensional attractor
surface. Therefore, the two parameters $k_{2}$ and $g_{2}$ determining the
early time expansions about $\tau=0$ discussed above are not convenient.
Instead we choose $\gamma\varepsilon(\tau_{0})$ and
$\gamma\tilde{\varepsilon}(\tau_{0})$, or rather $\gamma\varepsilon(\tau_{0})$
and the ratio $\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ as our
parameters. Clearly $\phi(\tau_{0})$ and $\tilde{\phi}(\tau_{0})$ should be
fine-tuned such that $\chi$ and $\tilde{\chi}$ reach $\sigma_{1}$ and
$\sigma_{2}$ as $\tau\rightarrow 0$.

Crucially, $\gamma\varepsilon(\tau_{0})$ is bounded from above for any choice
of $\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$, as otherwise we do
not obtain physical solutions. However, the physical energy densities and the
total energy density can be arbitrarily large. This is similar to the
equilibrium solutions studied in Kurkela:2018dku  where $\varepsilon$ and
$\tilde{\varepsilon}$ were bounded from above, but the total energy density
and pressure were unbounded.

As shown on the left panel of Fig. 4, the upper bound on
$\gamma\varepsilon(\tau_{0})$ is a function of the ratio
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ and it goes to $1$ as
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})\rightarrow\infty$. As
shown on the right panel of Fig. 4, the total energy density diverges as
$\gamma\varepsilon(\tau_{0})$ reaches its upper bound for any fixed ratio
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$.

To initialize on the attractor surface and find the upper bound on
$\gamma\varepsilon(\tau_{0})$ at the initial time
$\tau_{0}=10^{-3}\gamma^{1/4}$ we can proceed as follows. For each
$\gamma\varepsilon(\tau_{0})$ and
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$, we fix the values of
$\phi$ and $\tilde{\phi}$ by setting $\chi=\sigma_{1}$ and
$\tilde{\chi}=\sigma_{2}$ at $\tau=\tau_{0}$. Then with chosen precision
(approx. $10^{-3}$) we determine the upper limit of
$\gamma\varepsilon(\tau_{0})$ such that we obtain physical solutions. We then
tune $\phi$ and $\tilde{\phi}$ again until we find that $\chi=\sigma_{1}$ and
$\tilde{\chi}=\sigma_{2}$ at $\tau=10^{-6}\gamma^{1/4}$ up to a chosen
precision – this implies that our evolution is closer to the attractor
surface. We also recompute the upper bound on $\gamma\varepsilon(\tau_{0})$
for each $\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$. The new $\phi$
and $\tilde{\phi}$ for each $\gamma\varepsilon(\tau_{0})$ and
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$, and also the new upper
bound on $\gamma\varepsilon(\tau_{0})$ for each
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ differ from the previous
ones at a subpercent level giving estimate of our errors.

We can systematically correct this error by iterative fine-tuning. At each
step of the iteration, we fine tune $\phi$ and $\tilde{\phi}$ for each
$\gamma\varepsilon(\tau_{0})$ and
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ such that
$\chi=\sigma_{1}$ and $\tilde{\chi}=\sigma_{2}$ at progressively earlier
times. Practically, the very first step of iteration gives a very precise
estimate of initial values of $\phi$ and $\tilde{\phi}$ on the attractor
surface for each $\gamma\varepsilon(\tau_{0})$ and
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$, and also the upper
bound on $\gamma\varepsilon(\tau_{0})$ for each
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ if we initialize at
$\tau_{0}=10^{-3}\gamma^{1/4}$ with our choices of parameters. We also find
that the estimates of hydrodynamization times and other results reported below
change only sub-percent if we correct for our deviation from the attractor
surface systematically in the iterative procedure.

Figure 4: Left: Plot of upper bound on $\gamma\varepsilon(\tau_{0})$ for
different initial ratios
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$. The parameters used are
found in (30) and (31), and $\tau_{0}=10^{-3}\gamma^{1/4}$. Right:
$\gamma\varepsilon(\tau_{0})$ vs the total energy density, for the ratio
$\varepsilon/\tilde{\varepsilon}=2$ at $\tau_{0}$. The subsystem energy,
$\gamma\varepsilon(\tau_{0})$, saturates to its upper bound of $\sim 0.482$
for large total energy density, $\gamma E$.

###  3.3 Hydrodynamic expansion

At late time, the full system can be mapped to a hydrodynamic derivative
expansion with powers of $\tau^{-2/3}$ counting number of derivatives. The
hydrodynamic expansion is determined only by two dimensionless parameters,
which could be chosen for instance as

$\lim_{\tau\rightarrow\infty}\gamma^{2/3}\varepsilon\tau^{4/3}\quad{\rm
and}\quad\lim_{\tau\rightarrow\infty}\gamma^{2/3}\tilde{\varepsilon}\tau^{4/3}.$
(51)

(For reasons to be mentioned below, we will choose a different pair of
parameters.) Each pair of these parameters determines a unique curve on the
attractor surface. Thus any initial condition in the four-dimensional phase
space evolves to a specific curve on the attractor surface labelled by a
specific pair of parameters determining the late-time hydrodynamic expansion.
On the attractor surface itself, the pair of parameters $k_{2}$ and $g_{2}$
determining the early time expansion should relate to the pair of parameters
$\gamma\varepsilon(\tau_{0})$ and
$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ giving the
initialization at time $\tau_{0}$. In turn the latter pair of parameters
should match uniquely to a pair of parameters giving the late-time
hydrodynamic expansion. We postpone the discussion of the latter matching and
turn now to the hydrodynamic expansion itself.

The entire hydrodynamic expansion can be expressed in terms of the perfect
fluid components which give the leading $\tau$ dependence of $\varepsilon$ and
$\tilde{\varepsilon}$ at late time as follows

$\varepsilon=\epsilon_{1,pf}+\cdots,\quad\tilde{\varepsilon}=\epsilon_{2,pf}+\cdots,\quad\epsilon_{1,pf}=\epsilon_{10}\left(\frac{\tau}{\tau_{0}}\right)^{-4/3},\quad\epsilon_{2,pf}=\epsilon_{20}\left(\frac{\tau}{\tau_{0}}\right)^{-4/3},$
(52)

with $\cdots$ denoting viscous corrections and $\tau_{0}$ is an arbitrary
moment of time. Note that at late time the two fluids decouple. The expansion
dilutes the energy densities and therefore the mutual coupling via the
effective metrics. As a result, the leading order perfect fluid expansions
assume the usual forms for both systems. It is clear that the choice of
$\tau_{0}$ above is arbitrary. A different choice $\tilde{\tau}_{0}$ simply
rescales $\epsilon_{10}$ and $\epsilon_{20}$ as follows:

$\epsilon_{10}=\tilde{\epsilon}_{10}(\tilde{\tau}_{0}/\tau_{0})^{4/3},\quad\epsilon_{20}=\tilde{\epsilon}_{20}(\tilde{\tau}_{0}/\tau_{0})^{4/3}.$
(53)

We want to label hydrodynamic curves by variables which are dimensionless and
invariant under reprametrizations $\tau_{0}\rightarrow\tilde{\tau}_{0}$. It is
straightforward to see that the relevant variables are

$\displaystyle\alpha:=$
$\displaystyle\gamma(\epsilon_{20}\tau_{0}^{4/3})\sqrt{\epsilon_{10}\tau_{0}^{4/3}}=\gamma\epsilon_{20}\sqrt{\epsilon_{10}}\tau_{0}^{2},$
(54) $\displaystyle\beta:=$
$\displaystyle\gamma(\epsilon_{10}\tau_{0}^{4/3})\sqrt{\epsilon_{20}\tau_{0}^{4/3}}=\gamma\epsilon_{10}\sqrt{\epsilon_{20}}\tau_{0}^{2}.$
(55)

These two parameters uniquely label each attractor (hydrodynamic) curve as
shown below. It is also useful to define

$\mathfrak{r}:=(\epsilon_{20}/\epsilon_{10})^{1/4}=\sqrt{\alpha/\beta}=\lim_{\tau\rightarrow\infty}{(\tilde{\varepsilon}/\varepsilon)^{1/4}}.$
(56)

For the sake of convenience, we define dimensionless time:

$\omega_{1}:=\epsilon_{1,pf}^{1/4}\tau,\quad\omega_{2}:=\epsilon_{2,pf}^{1/4}\tau.$
(57)

Note however that $\omega_{1}$ and $\omega_{2}$ are not independent time
variables. In fact

$\omega_{2}=\mathfrak{r}\omega_{1}.$ (58)

The hydrodynamic expansion in $\tau$, obtained by solving all coupling
equations, the two MIS equations and the two conservation equations, can be
readily stated in terms of these variables as follows

$\displaystyle\varepsilon^{1/4}\tau$
$\displaystyle=\omega_{1}-\frac{2}{3}C_{\eta}+\frac{3\alpha-4C_{\eta}C_{\tau}}{18}\frac{1}{\omega_{1}}+\mathcal{O}(\omega_{1}^{-2}),$
(59) $\displaystyle\chi$
$\displaystyle=\frac{4C_{\eta}}{3\omega_{1}}+\frac{8C_{\eta}C_{\tau}}{9\omega_{1}^{2}}+\mathcal{O}(\omega_{1}^{-3}),$
(60) $\displaystyle\tilde{\varepsilon}^{1/4}\tau$
$\displaystyle=\omega_{2}-\frac{2}{3}\tilde{C}_{\eta}+\frac{3\beta-4\tilde{C}_{\eta}\tilde{C}_{\tau}}{18}\frac{1}{\omega_{2}}+\mathcal{O}(\omega_{2}^{-2}),$
(61) $\displaystyle\tilde{\chi}$
$\displaystyle=\frac{4\tilde{C}_{\eta}}{3\omega_{2}}+\frac{8\tilde{C}_{\eta}\tilde{C}_{\tau}}{9\omega_{2}^{2}}+\mathcal{O}(\omega_{2}^{-3}),$
(62)

The above equations show all terms up to second order in derivatives. Defining

$w_{1}:=\varepsilon^{1/4}\tau,\quad w_{2}:=\tilde{\varepsilon}^{1/4}\tau$ (63)

we can invert (59)

$\omega_{1}=w_{1}+\frac{2}{3}C_{\eta}-\frac{3\alpha-4C_{\eta}C_{\tau}}{18}\frac{1}{w_{1}}+\mathcal{O}(\omega_{1}^{-2}).$
(64)

In terms of $w_{1}$ we can readily express the hydrodynamic expansion as a
phase space curve

$\displaystyle\chi$
$\displaystyle=\frac{4C_{\eta}}{3w_{1}}+\frac{8C_{\eta}C_{\tau}}{9w_{1}^{2}}+\mathcal{O}(w_{1}^{-3}),$
(65) $\displaystyle w_{2}$
$\displaystyle=\sqrt{\frac{\alpha}{\beta}}w_{1}-\frac{2}{3}\left(\tilde{C}_{\eta}-\sqrt{\frac{\alpha}{\beta}}C_{\eta}\right)$
(66)
$\displaystyle+\left(\sqrt{\frac{\beta}{\alpha}}\frac{3\beta-4\tilde{C}_{\eta}\tilde{C}_{\tau}-24\tilde{C}_{\eta}^{2}}{18}-\sqrt{\frac{\alpha}{\beta}}\frac{3\alpha-
C_{\eta}C_{\tau}-24C_{\eta}^{2}}{18}\right)\frac{1}{w_{1}}+\mathcal{O}(w_{1}^{-2}),$
$\displaystyle\tilde{\chi}$
$\displaystyle=\sqrt{\frac{\beta}{\alpha}}\frac{4\tilde{C}_{\eta}}{3w_{1}}+\frac{\beta}{\alpha}\frac{8\tilde{C}_{\eta}\left(\tilde{C}_{\eta}+C_{\eta}\sqrt{\frac{\beta}{\alpha}}+\tilde{C}_{\tau}\right)}{9w_{1}^{2}}+\mathcal{O}(w_{1}^{-3}).$
(67)

With $w_{1}$, $w_{2}$, $\chi$ and $\tilde{\chi}$ as the phase space variables,
the phase space curves representing the hydrodynamic expansion (or the
attractor) are then given by the above expressions, $\chi(w_{1})$,
$w_{2}(w_{1})$ and $\tilde{\chi}(w_{1})$. As manifest from above, each curve
is labelled by two dimensionless, time-reparametrisation invariant parameters,
$\alpha$ and $\beta$.

The hydrodynamic late-time expansion of the components of the effective
metrics are:

$\displaystyle a$
$\displaystyle=1-\frac{1}{2}\alpha\frac{1}{\omega_{1}^{2}}+\frac{4}{3}\tilde{C}_{\eta}\sqrt{\alpha\beta}\frac{1}{\omega_{1}^{3}}+\mathcal{O}(\omega_{1}^{-4}),$
(68) $\displaystyle b$
$\displaystyle=1+\frac{1}{6}\alpha\frac{1}{\omega_{1}^{2}}+\mathcal{O}(\omega_{1}^{-4}),$
$\displaystyle c$
$\displaystyle=1+\frac{1}{6}\alpha\frac{1}{\omega_{1}^{2}}-\frac{4}{3}\tilde{C}_{\eta}\sqrt{\alpha\beta}\frac{1}{\omega_{1}^{3}}+\mathcal{O}(w_{1}^{-4}),$
(69) $\displaystyle\tilde{a}$
$\displaystyle=1-\frac{1}{2}\frac{\beta^{2}}{\alpha}\frac{1}{\omega_{1}^{2}}+\frac{4}{3}C_{\eta}\frac{\beta^{2}}{\alpha}\frac{1}{\omega_{1}^{3}}+\mathcal{O}(\omega_{1}^{-4}),$
(70) $\displaystyle\tilde{b}$
$\displaystyle=1+\frac{1}{6}\frac{\beta^{2}}{\alpha}\frac{1}{\omega_{1}^{2}}+\mathcal{O}(\omega_{1}^{-4}),$
$\displaystyle\tilde{c}$
$\displaystyle=1+\frac{1}{6}\frac{\beta^{2}}{\alpha}\frac{1}{\omega_{1}^{2}}-\frac{4}{3}C_{\eta}\frac{\beta^{2}}{\alpha}\frac{1}{\omega_{1}^{3}}+\mathcal{O}(\omega_{1}^{-4}).$
(71)

Remarkably, up to $\mathcal{O}(w_{1}^{-4})$ the volume factors are:

$ab^{2}c=1+\mathcal{O}(w_{1}^{-4}),\quad\tilde{a}\tilde{b}^{2}\tilde{c}=1+\mathcal{O}(w_{1}^{-4}).$
(72)

We can readily compute the full system energy density which behaves as

$\mathcal{E}^{1/4}\tau=\omega-\frac{2}{3}\frac{C_{\eta_{1}}\epsilon_{10}^{3/4}+C_{\eta_{2}}\epsilon_{20}^{3/4}}{(\epsilon_{10}+\epsilon_{20})^{3/4}}+\mathcal{O}(\omega^{-1}),$
(73)

where

$\omega:=\mathcal{E}_{pf}^{1/4}\tau=(\epsilon_{1,pf}+\epsilon_{2,pf})^{1/4}\tau.$
(74)

We find that the full energy-momentum tensor can be mapped to the standard
hydrodynamic Bjorken flow expansion (in the usual Milne metric) of a single
fluid, but with effective equation of state and effective transport
coefficients dependent on the parameters $\alpha$ and $\beta$.

We can explicitly map the full system to a single expanding fluid at late time
and compare to the above discussion. For a single fluid, the evolution of the
energy density in the Milne background metric takes the generic form (up to
first order)

$\mathcal{E}^{\prime}=-\frac{1}{\tau}\left(\mathcal{E}+\mathcal{P}(\mathcal{E})\right)+\frac{1}{\tau^{2}}\left(\zeta(\mathcal{E})+\frac{2}{3}\eta(\mathcal{E})\right),$
(75)

where $\zeta$ and $\eta$ are the bulk and shear viscosity, respectively.
Furthermore, the trace of the (first-order) hydrodynamic energy-momentum
tensor takes the form

$T^{\mu}_{\,\mu}=3\mathcal{P}(\mathcal{E})-\mathcal{E}-\frac{3}{\tau}\zeta(\mathcal{E}).$
(76)

We can assume that

$\displaystyle\mathcal{P}(\mathcal{E})$
$\displaystyle=\frac{\mathcal{E}}{3}+k\mathcal{E}^{2}+\cdots$ (77)
$\displaystyle\eta(\mathcal{E})$
$\displaystyle=\kappa\mathcal{E}^{3/4}+\cdots,\quad\zeta(\mathcal{E})=\tilde{\kappa}\mathcal{E}^{3/4}+\cdots.$
(78)

with $k$, $\kappa$ and $\tilde{\kappa}$ constants. Solving (75) with inputs
from (77) in the large proper time expansion, we obtain

$\displaystyle\mathcal{E}=$
$\displaystyle\mathcal{E}_{0}\left(\frac{\tau_{0}}{\tau}\right)^{4/3}-\frac{\tau_{0}}{2\tau^{2}}(2\kappa+3\tilde{\kappa})\mathcal{E}_{0}^{3/4}$
$\displaystyle+\frac{3}{32\tau^{8/3}}\left(\mathcal{E}_{0}^{1/2}\tau_{0}^{2/3}(4\kappa^{2}+12\kappa\tilde{\kappa}+9\tilde{\kappa}^{2})-8\mathcal{E}_{0}^{2}\tau_{0}^{8/3}k\right)+\cdots,$
(79)

and also (76)

${\rm
Tr}(T)=-\frac{3}{\tau^{2}}\mathcal{E}_{0}^{3/4}\tau_{0}\tilde{\kappa}+\frac{3}{8\tau^{8/3}}\left(\mathcal{E}_{0}^{1/2}\tau_{0}^{2/3}(6\kappa\tilde{\kappa}+9\tilde{\kappa}^{2})+8\mathcal{E}_{0}^{2}\tau_{0}^{8/3}k\right)+\cdots.$
(80)

We can determine the constants $k$, $\kappa$ and $\tilde{\kappa}$ simply by
matching the above expansions to those obtained by solving our system of
equations. We readily see that these constants depend on the specific curve on
the attractor surface. Comparing Eq. (3.3) with Eq. (73), we explicitly find
that the effective viscosity of the full system is

$C_{\eta,eff}=\frac{C_{\eta}\epsilon_{10}^{3/4}+\tilde{C}_{\eta}\epsilon_{20}^{3/4}}{(\epsilon_{10}+\epsilon_{20})^{3/4}}=\frac{C_{\eta}\beta^{3/4}+\tilde{C}_{\eta}\alpha^{3/4}}{(\alpha^{2}+\beta^{2})^{3/4}}$
(81)

which depends manifestly on which curve on the attractor surface the full
system evolves to at late time.In fact $\kappa$ can be obtained from (81):

$\kappa=\frac{4}{3}C_{\eta,eff}$ (82)

Similarly, the hydrodynamic expansion of our hybrid system yields that

${\rm
Tr}(T)=-\frac{8\tau_{0}^{8/3}}{3\tau^{8/3}}\gamma\epsilon_{10}\epsilon_{20}+\cdots.$
(83)

We can immediately conclude that

  * •

$\mathcal{E}_{0}=\epsilon_{10}+\epsilon_{20}.$

  * •

$k=-\frac{1}{3}\gamma\frac{\epsilon_{10}\epsilon_{20}}{\mathcal{E}_{0}^{2}}=-\frac{1}{3}\gamma\frac{\epsilon_{10}\epsilon_{20}}{(\epsilon_{10}+\epsilon_{20})^{2}}=-\frac{1}{3}\gamma\frac{\alpha^{2}\beta^{2}}{(\alpha^{2}+\beta^{2})^{2}},$

  * •

and $\tilde{\kappa}=0,$

i.e. the effective bulk viscosity vanishes since the $\tau^{-2}$ term vanishes
in (83).

Since $k$ and $\kappa$ both depend on $\alpha$ and $\beta$, it’s clear that
the effective equation of state and the transport coefficients of the single
fluid description of the full system depends on the specific curve on the
attractor surface.

In order to interpolate to hydrodynamics, it is useful to do a doubly
logarithmic plot for $\chi$, $\tilde{\chi}$, $\varepsilon$ and
$\tilde{\varepsilon}$. Perfect hydrodynamics begins when the logarithmic plots
of $\varepsilon(\tau)$ and $\tilde{\varepsilon}(\tau)$ become parallel to each
other with slope $-4/3$, i.e. when they both reach the perfect fluid
hydrodynamic tail $1/\tau^{4/3}$, which is evident in Fig. 5 for the first-
order hydrodynamic tails of $\varepsilon$ and $\tilde{\varepsilon}$ for the
first and second attractors, respectively.

Figure 5: The hydrodynamic tails for the attractor with
$(\varepsilon,\tilde{\varepsilon},{\phi},{\tilde{\phi}})(\tau_{0})$ given by
$(0.18,0.18,0.1000210481,0.05015396018505)$

At late time

$\ln\varepsilon\approx-(4/3)\ln\tau+\ln(\epsilon_{10}\tau_{0}^{4/3}),\quad\ln\tilde{\varepsilon}\approx-(4/3)\ln\tau+\ln(\epsilon_{20}\tau_{0}^{4/3}).$
(84)

Therefore the intercepts of the energy densities give us the time-
reparametrisation invariant $\epsilon_{10}\tau_{0}^{4/3}$ and
$\epsilon_{10}\tau_{0}^{4/3}$ with the slopes fitting well to $-4/3$. We
obtain for the attractor

$\epsilon_{10}\tau_{0}^{4/3}\approx
1.346,\,\,\epsilon_{20}\tau_{0}^{4/3}\approx
0.370,\,\,\alpha=0.424,\,\,\beta=0.801.$ (85)

First-order hydrodynamics begins when the logarithmic plots of $\chi(\tau)$
and $\tilde{\chi}(\tau)$ become parallel to each other with slope $-2$, i.e.
when they both reach the hydrodynamic tail $1/\tau^{2}$ given by first-order
hydrodynamics (note in the perfect fluid limit $\chi$ vanishes). This is shown
in Fig. 6 for the first-order hydrodynamic tails of $\chi$ and $\tilde{\chi}$
for the first and second attractors, respectively.

Figure 6: The hydrodynamic tail showing the onset of first-order hydrodynamics
with the same values as in Fig. 5.

We can also match using leading order hydrodynamic expansion of $\chi$
($\tilde{\chi}$):

$\displaystyle\ln{\chi}\approx-2\ln\tau+\ln\frac{16C_{\eta}}{9}+3/4\ln(\epsilon_{10}\tau_{0}^{4/3}),$
$\displaystyle\ln{\tilde{\chi}}\approx-2\ln\tau+\ln\frac{16\tilde{C}_{\eta}}{9}+3/4\ln(\epsilon_{20}\tau_{0}^{4/3}).$
(86)

As a consistency check, we see that the intercepts reproduce
$\epsilon_{10}\tau_{0}^{4/3}$ and $\epsilon_{20}\tau_{0}^{4/3}$ to a high
accuracy.

###  3.4 Late time approach to attractor: fluctuation analysis

The approach to the attractor at late time can be simply understood by
perturbing around the late-time hydrodynamic expansion computed above.

The leading order linear fluctuations of energy densities, anisotropies and
effective metric variables about the hydrodynamic expansion, i.e. the
attractor solution at late time are defined as follows:

$\displaystyle\varepsilon(\hat{\tau})$
$\displaystyle={\epsilon_{10}}{\hat{\tau}}^{-4/3}+{\hat{\tau}}^{-4/3}\left({\epsilon_{11}}\hat{\tau}^{-2/3}+{\epsilon_{12}}\hat{\tau}^{-4/3}\right)+\mathcal{O}\left(\hat{\tau}^{-10/3}\right)+\delta\varepsilon(\hat{\tau}),$
(87) $\displaystyle\phi(\hat{\tau})$
$\displaystyle=f_{10}+\left(f_{11}\hat{\tau}^{-2/3}+f_{12}\hat{\tau}^{-4/3}+f_{13}\hat{\tau}^{-2}++f_{14}\hat{\tau}^{-8/3}\right)+\mathcal{O}\left(\hat{\tau}^{-10/3}\right)+\delta\phi(\hat{\tau}),$
(88) $\displaystyle a(\hat{\tau})$
$\displaystyle=1+\hat{\tau}^{-2/3}\left(A_{11}\hat{\tau}^{-2/3}+A_{12}\hat{\tau}^{-4/3}\right)+\mathcal{O}\left(\hat{\tau}^{-10/3}\right)+\delta
a(\hat{\tau}),$ (89) $\displaystyle b(\hat{\tau})$
$\displaystyle=1+\hat{\tau}^{-2/3}\left(B_{11}\hat{\tau}^{-2/3}+B_{12}\hat{\tau}^{-4/3}\right)+\mathcal{O}\left(\hat{\tau}^{-10/3}\right)+\delta
b(\hat{\tau}),$ (90) $\displaystyle c(\hat{\tau})$
$\displaystyle=\tau_{0}\hat{\tau}+{\hat{\tau}}^{1/3}\left(C_{11}\hat{\tau}^{-2/3}+C_{12}\hat{\tau}^{-4/3}\right)+\mathcal{O}\left(\hat{\tau}^{-10/3}\right)+\delta
c(\hat{\tau}),$ (91)

and similarly for the tilded sector. Above $\hat{\tau}:=\tau/\tau_{0}$ and
$\delta\tilde{\varepsilon}(\hat{\tau})$, $\delta\phi(\hat{\tau})$, etc. denote
the fluctuations about the hydrodynamic expansion. We have truncated to the
first few terms of the hydrodynamic expansions which can be found explicitly
in Appendix A, because only these are sufficient to obtain the leading order
behavior of the fluctuations. The fluctuations have two independent non-
hydrodynamic modes capturing the missing information of the initial conditions
in the four-dimensional phase space, and two hydrodynamic modes which span the
two-dimensional attractor surface. Since the non-hydrodynamic modes, which do
not admit a power series expansion $\hat{\tau}^{-2/3}$, tell us how a generic
evolution relaxes towards the attractor surface, we will focus on finding
these explicitly here.

One can eliminate the fluctuations in six effective metric variables in terms
of other four fluctuations ($\delta\varepsilon(\hat{\tau})$,
$\delta\tilde{\varepsilon}(\hat{\tau})$, $\delta\phi(\hat{\tau})$,
$\delta\tilde{\phi}(\hat{\tau})$) by using six metric coupling equations. This
leaves us with four first-order coupled ordinary differential equations in
four variables, which are given in Appendix A. One can further solve these
equations until one is left with a pair of decoupled fourth-order differential
equations for the anisotropy fluctuations. Depending on the relative values of
parameters, the fluctuations, $\delta\phi(\hat{\tau})$ and
$\delta\tilde{\phi}(\hat{\tau})$, have different leading order forms, which we
detail below.

In the case when
$C_{\tau}\>{\epsilon_{20}}^{1/4}\neq\tilde{C}_{\tau}\>{\epsilon_{10}}^{1/4}$,
the two independent non-hydrodynamic solutions (at leading order) are:

$\displaystyle\delta\phi(\hat{\tau})$
$\displaystyle=m_{1}\;e^{-\frac{3{\epsilon_{10}}^{1/4}\;\tau_{0}}{2C_{\tau}}\hat{\tau}^{2/3}}\hat{\tau}^{-\frac{4}{3}+\frac{2C_{\eta}}{3C_{\tau}}}+m_{2}\;e^{-\frac{3{\epsilon_{20}}^{1/4}\;\tau_{0}}{2\tilde{C}_{\tau}}\hat{\tau}^{2/3}}\hat{\tau}^{-\frac{10}{3}+\frac{2\tilde{C}_{\eta}}{3\tilde{C}_{\tau}}},$
(92) $\displaystyle\delta\tilde{\phi}(\hat{\tau})$
$\displaystyle=n_{1}\;e^{-\frac{3{\epsilon_{10}}^{1/4}\;\tau_{0}}{2C_{\tau}}\hat{\tau}^{2/3}}\hat{\tau}^{-\frac{10}{3}+\frac{2C_{\eta}}{3C_{\tau}}}+n_{2}\;e^{-\frac{3{\epsilon_{20}}^{1/4}\;\tau_{0}}{2C_{\tau_{2}}}\hat{\tau}^{2/3}}\hat{\tau}^{-\frac{4}{3}+\frac{2\tilde{C}_{\eta}}{3\tilde{C}_{\tau}}}$
(93)

where $m_{1},m_{2},n_{1},n_{2}$ are constants, of which only two are
independent:

$\displaystyle m_{1}$ $\displaystyle=$ $\displaystyle
n_{1}\frac{3C_{\eta}\left(\tilde{C}_{\tau}\;{\epsilon_{10}}^{1/4}-C_{\tau}\;{\epsilon_{20}}^{1/4}\right)}{2\tilde{C}_{\eta}C_{\tau}^{2}{\epsilon_{20}}\gamma}\tau_{0},$
$\displaystyle n_{2}$ $\displaystyle=$ $\displaystyle
m_{2}\frac{3\tilde{C}_{\eta}\left(-\tilde{C}_{\tau}\;{\epsilon_{10}}^{1/4}+C_{\tau}\;{\epsilon_{20}}^{1/4}\right)}{2C_{\eta}\tilde{C}_{\tau}^{2}{\epsilon_{10}}\gamma}\tau_{0}.$
(94)

If $C_{\tau}\>{\epsilon_{20}}^{1/4}>\tilde{C}_{\tau}\>{\epsilon_{10}}^{1/4}$,
i.e. if

$\mathfrak{r}=\sqrt{\alpha/\beta}>\tilde{C}_{\tau}/C_{\tau}$ (95)

then

$\exp{\left(-\frac{3{\epsilon_{10}}^{1/4}\;\tau_{0}}{2C_{\tau}}\hat{\tau}^{2/3}\right)}=\exp{\left(-\frac{3{\epsilon_{10}}^{1/4}\;\tau_{0}^{1/3}}{2C_{\tau}}\tau^{2/3}\right)}=\exp{\left(-\frac{3}{2C_{\tau}}\left(\frac{\beta}{\sqrt{\alpha\gamma}}\right)^{1/3}\tau^{2/3}\right)}$
(96)

will be the dominant exponential and give the relaxation towards the
attractor. The powers of $\hat{\tau}$ accompanying this exponential in (92)
and (93) imply that $\delta\tilde{\phi}(\hat{\tau})$ will decay faster towards
the attractor surface than $\delta\phi(\hat{\tau})$ in this case. This implies
that the strongly coupled sector will hydrodynamize faster than the weakly
coupled one. If
$C_{\tau}\>{\epsilon_{20}}^{1/4}<\tilde{C}_{\tau}\>{\epsilon_{10}}^{1/4}$,
i.e. if

$\mathfrak{r}=\sqrt{\alpha/\beta}<\tilde{C}_{\tau}/C_{\tau},$ (97)

then

$\exp{\left(-\frac{3{\epsilon_{20}}^{1/4}\;\tau_{0}}{2\tilde{C}_{\tau}}\hat{\tau}^{2/3}\right)}=\exp{\left(-\frac{3{\epsilon_{20}}^{1/4}\;\tau_{0}^{1/3}}{2\tilde{C}_{\tau}}\tau^{2/3}\right)}=\exp{\left(-\frac{3}{2\tilde{C}_{\tau}}\left(\frac{\alpha}{\sqrt{\beta\gamma}}\right)^{1/3}\tau^{2/3}\right)}$
(98)

is the dominant mode. By an analogous argument, $\delta\phi(\hat{\tau})$ will
decay faster than $\delta\tilde{\phi}(\hat{\tau})$ implying that the weakly
coupled sector will hydrodynamize faster than the strongly coupled one.

Full numerical solutions do confirm that indeed when (95) is realized, the
strongly coupled sector hydrodynamizes faster and the reverse occurs when (97)
holds. However, realizing (97) requires choosing large ratios of
$\tilde{\varepsilon}(\tau_{0})/\varepsilon(\tau_{0})$ such as 100,000 at say
$\tau_{0}={\gamma^{1/4}}$. This is also unnatural because this requires
extreme fine-tuning at early time, because of the bottom-up behavior on the
attractor surface in which $\varepsilon$ goes to a constant while
$\tilde{\varepsilon}$ vanishes as $\tau\rightarrow 0$.

In the case when
$C_{\tau}\;{\epsilon_{20}}^{1/4}=\tilde{C}_{\tau}\;{\epsilon_{10}}^{1/4}$,
i.e.

$\mathfrak{r}=\sqrt{\alpha/\beta}=\tilde{C}_{\tau}/C_{\tau},$ (99)

the two independent solutions to leading order are:

$\displaystyle\delta\phi(\hat{\tau})$
$\displaystyle=e^{-\frac{3{\epsilon_{10}}^{1/4}\tau_{0}}{2C_{\tau}}\hat{\tau}^{2/3}}\left(m_{1}\;\hat{\tau}^{-\frac{4}{3}+\frac{2C_{\eta}}{3C_{\tau}}}+m_{2}\;\hat{\tau}^{-\frac{8}{3}+\frac{2\tilde{C}_{\eta}}{3\tilde{C}_{\tau}}}\right),$
(100) $\displaystyle\delta\tilde{\phi}(\hat{\tau})$
$\displaystyle=e^{-\frac{3{\epsilon_{10}}^{1/4}\tau_{0}}{2C_{\tau}}\hat{\tau}^{2/3}}\left(n_{1}\;\hat{\tau}^{-\frac{8}{3}+\frac{2C_{\eta}}{3C_{\tau}}}+n_{2}\;\hat{\tau}^{-\frac{4}{3}+\frac{2\tilde{C}_{\eta}}{3\tilde{C}_{\tau}}}\right),$
(101)

where $m_{1}$, $m_{2}$, $n_{1}$, and $n_{2}$ are arbitrary constants only two
of which are independent:

$\displaystyle
m_{1}=-\frac{C_{\eta}C_{\tau}^{2}\left(-\tilde{C}_{\eta}C_{\tau}+C_{\eta}\tilde{C}_{\tau}-2C_{\tau}\tilde{C}_{\tau}\right)}{\tilde{C}_{\eta}\tilde{C}_{\tau}^{4}{\epsilon_{10}}\gamma}n_{1},$
$\displaystyle
n_{2}=-\frac{\tilde{C}_{\eta}\left(\tilde{C}_{\eta}C_{\tau}-C_{\eta}\tilde{C}_{\tau}-2C_{\tau}\tilde{C}_{\tau}\right)}{C_{\eta}\tilde{C}_{\tau}^{4}{\epsilon_{10}}\gamma}m_{2}.$
(102)

###  3.5 Matching early to late times

Now we are prepared to study how parameters setting initial conditions on the
Bjorken flow attractor surface map to $\alpha$ and $\beta$, which are defined
in (54) and (55), and which determine the hydrodynamic expansion at late time.
In particular, we find that they scale as a power of the total energy density
at initial time $\tau_{0}$ on the attractor surface. For $r=2$ and
initialization time $\tau_{0}=\gamma^{1/4}$, the scaling of $\alpha$ and
$\beta$ is found to be

$\displaystyle\alpha$ $\displaystyle\sim(\gamma\mathcal{E})^{3/2},$ (103)
$\displaystyle\beta$ $\displaystyle\sim(\gamma\mathcal{E})^{3/2},$ (104)

where $\mathcal{E}$ is the total initial energy density. Using the definitions
of the total entropy density and the individual subsystem entropy densities as

$\displaystyle S_{tot}$ $\displaystyle=s+\tilde{s},$ (105) $\displaystyle s$
$\displaystyle=\frac{4}{3\tau}\varepsilon^{3/4}b^{2}c,$ (106)
$\displaystyle\tilde{s}$
$\displaystyle=\frac{4}{3\tau}\tilde{\varepsilon}^{3/4}\tilde{b}^{2}\tilde{c},$
(107)

we can readily see from the hydrodynamic expansions discussed above that the
individual system entropy densities behave at late time as

$\displaystyle\lim_{\tau\rightarrow\infty}\tau
s=\beta/\sqrt{\alpha},\quad\lim_{\tau\rightarrow\infty}\tau\tilde{s}=\alpha/\sqrt{\beta},$
(108)

which means that the total entropy scales like

$\displaystyle\lim_{\tau\rightarrow\infty}\tau S_{\rm tot}\equiv(\tau S_{\rm
tot})_{\rm hydro}$ $\displaystyle\sim(\gamma\mathcal{E})^{3/4}.$ (109)

We can compare this to the pocket formula in Giacalone:2019ldn  giving the
final state particle multiplicity as a power of the total initial energy
density for single conformal fluid attractors via $(S\tau)_{\rm
hydro}\sim(e\tau)_{0}^{2/3}$ to see that the interactions due to the effective
metric coupling in our hybrid model increases the scaling exponent from $2/3$
to $3/4$ for $\tau_{0}=\gamma^{1/4}$. The pocket formula of Giacalone:2019ldn
states that the total particle multiplicity per unit rapidity444Here, the
rapidity, $\eta$, should not to be confused with the shear viscosity. follows

$\displaystyle\frac{dN}{d\eta}\approx(S\tau)_{\rm hydro}.$ (110)

We would need more microscopic inputs to relate our results to particle
multiplicities, and therefore it is beyond the scope of the present work.
Interestingly, in the case of one fluid in Bjorken flow, the attractor
solution has been shown to provide a maximum for dilepton production (while
the repulsor/disperser represents a minimum) Naik:2021yph . It would be
interesting to do a similar study in our context, however we would then need
to replace the perturbative MIS sector with a kinetic theory.

In Fig. 7, we plot the full system entropy against proper time. We notice that
rather generically, the entropy increases as a function of $\tau$, before
flattening out once interactions have frozen out. Interestingly, we see that
the entropies of the individual subsystem can change in relative size at
earlier $\tau$ as shown in the inset and the right panel of Fig. 7, with the
soft (strongly self-interacting) sector dominating the entropy production at
later time.

Figure 7: The full system and the two subsystem entropies as a function of
proper time, see text for the parameters used. The inset in the left figure
shows the time when the second subsystem has a larger share in entropy density
over the first one, i.e. when the strongly coupled sector begins to dominate
the entropy production. Right: the ratio of the two subsystem entropies. The
parameters used to generate the figure are given by
$\tilde{C}_{\eta}=\frac{1}{4\pi},\tilde{C}_{\tau}=\frac{2-ln(2)}{2\pi},C_{\eta}=2\tilde{C}_{\eta},\tilde{C}_{\tau}=4C_{\tau},\phi(\tau_{0})=0.0220931202,\tilde{\phi}(\tau_{0})=0.0297056408397,\varepsilon(\tau_{0})=0.05,\tilde{\varepsilon}(\tau_{0})=0.1$.

We can also study the trace of the full energy momentum tensor, which encodes
the interaction energy. For late times, we find that the trace falls like the
typical energy density,

$\frac{|{\rm
Tr}T|}{\mathcal{E}}=\frac{\mathcal{E}-P_{L}-2P_{T}}{\mathcal{E}}\sim\tau^{-4/3},$

as follows from (73) and (83). An illustrative example for the full time-
evolution of the interaction energy is shown in Fig. 8. We find that the
interaction energy is non-vanishing only in the intermediate time-scales when
both systems hydrodynamize, however it vanishes both at early and late times.

Figure 8: The trace of the full energy momentum tensor as a function of the
total energy density, which acts as a proxy for the effective temperature.
This plot is for an initial energy density of the hard sector of
$\varepsilon=0.8$ and $\varepsilon/\tilde{\varepsilon}=40$. The green dot
denotes the value at $\tau_{0}=10^{-3}$ and the red point is at $\tau=1000$.

###  3.6 Borel resummation of the hydrodynamic expansion

The full system at late time behaves as a single fluid with the transport
coefficient being determined by which attractor curve on the 2D attractor
surface the system goes to as discussed above. The curves on the attractor
surface are parameterized by the dimensionless parameters, $\alpha$ and
$\beta$ (defined in (54) and (55)) which govern the late time hydrodynamic
expansion and thus the effective transport coefficients of the full system.
Here we aim to understand the hydrodynamic gradient expansion and see how its
Borel resummation can describe the attractor surface in terms of these two
parameters.

We will expand the anisotropies, at late time in powers of inverse proper
time, $x=\tau^{-2/3},$

$\displaystyle\chi=\frac{4C_{\eta}}{3}\Big{(}\frac{\sqrt{\alpha\gamma}}{\beta}\Big{)}^{1/3}x+\frac{8}{9}\Big{(}\frac{\alpha\gamma}{\beta^{2}}\Big{)}^{1/3}(C_{\eta}^{2}+C_{\eta}C_{\tau})x^{2}+\mathcal{O}(x^{3})\,,$
$\displaystyle\tilde{\chi}=\frac{4\tilde{C}_{\eta}}{3}\Big{(}\frac{\sqrt{\beta\gamma}}{\alpha}\Big{)}^{1/3}x+\frac{8}{9}\Big{(}\frac{\beta\gamma}{\alpha^{2}}\Big{)}^{1/3}(\tilde{C}_{\eta}^{2}+\tilde{C}_{\eta}\tilde{C}_{\tau})x^{2}+\mathcal{O}(x^{3})\,.$
(111)

Let’s use the notation $\xi_{i}$ with $i=1,2$ and $\chi_{1}=\chi$ and
$\chi_{2}=\tilde{\chi}$. At arbitrarily high order the coefficients of these
two series expansions $r_{i,n}$ for $\chi_{i}$ show factorial growth following

$\frac{r_{i,n}}{r_{i,n+1}}=\frac{\xi_{i,0}}{n+1}\left(1+\frac{\gamma_{i}+1}{n+1}+\mathcal{O}(n^{-2})\right)$

indicating the series is asymptotic in nature with zero radius of convergence,
which is seen in subsystems in the top panels of Fig. 9. The parameters
$\xi_{i,0}$ and $\gamma_{i}$ for the weak and strong systems depend on
$\alpha$ and $\beta$.

To promote the divergent series to a well-defined function the series needs to
have finite radius of convergence. This is successfully done by adopting the
standard technique of Borel resummation Aniceto:2011nu . The Borel transform
of the series is given by

$\displaystyle\chi_{B,i}(\xi)=\sum_{n=0}^{\infty}\frac{r_{i,n}}{n!}{\xi}^{n}$
(112)

$r_{i,n}$ are the coefficients from the series (3.6). This Borel transformed
function has a finite radius of convergence about the origin in the complex
$\xi$ plane (i.e. the Borel plane) and is related to the original series via
Laplace transform. However the function possesses poles in the complex plane
as a consequence of the finite radius of convergence. Hence to perform the
integral one has to do analytic continuation of the function
$\chi_{B,i}({\xi})$. Since we only have access to a finite number of terms, we
will need to approximate the integrand, which we do via Padé approximant

$\chi_{B,i}({\xi})\approx
P_{N,M,i}({\xi})=\frac{\sum_{i=0}^{N}n_{i}{\xi}^{i}}{1+\sum_{j=1}^{M}m_{j}{\xi}^{j}}$
(113)

where the coefficients of the Padé approximant are fixed to agree with series
coefficients (112). The Padé series can go up to arbitrary orders with a
constraint $N+M=k$, where $k$ is the order up to which Borel transformed
anisotropy is truncated. Most commonly, one considers the symmetric case with
$N=M=k/2$. In our case we have taken the order $k$ to be as high as $500$ with
precision set to thousand digits.

The singularities of the function $\chi_{B,i}({\xi})$ appear as concentrated
poles of the Padé-approximant in the Borel plane and they depend on the
dimensionless parameter $\alpha$ and $\beta$. In Fig. 9, we show the pole
structure of both sectors in the Borel plane. The poles start from a point
nearest to the origin, namely $\xi_{1,0}$ and $\xi_{2,0}$ for $\chi_{B}(\xi)$
and $\tilde{\chi}_{B}(\xi)$ respectively, and accumulate along the real axis
giving an image of branch cut as in the case of decoupled MIS Heller:2015dha .
The starting value of the branch cut is proportional to the inverse of the
slope of the coefficients shown in Fig. 9 and can be identified with the non-
hydrodynamic mode that decays exponentially at late-times. From the analysis
in Sec. 3.4 we readily see that
$\xi_{1,0}=\xi_{2,0}=\frac{3}{2C_{\tau}}\Big{(}\frac{\beta}{\sqrt{\alpha\gamma}}\Big{)}^{1/3}=3\epsilon_{10}^{1/4}\tau_{0}^{1/3}/(2C_{\tau})=\xi_{0}$
when the condition (95) is satisfied, and
$\xi_{1,0}=\xi_{2,0}=\frac{3}{2C_{\tau}}\Big{(}\frac{\alpha}{\sqrt{\beta\gamma}}\Big{)}^{1/3}=3\epsilon_{20}^{1/4}\tau_{0}^{1/3}/(2C_{\tau})=\xi_{0}$
when the condition (97) holds.

With the analytic continuation of the Borel transformed function, we can now
perform the inverse integration Aniceto:2011nu ; Heller:2015dha ; Casalderrey-
Solana:2017zyh

$\displaystyle\chi_{i}(x)=\frac{1}{x}\int_{\mathcal{C}}e^{-\xi/x}P_{N,M,i}({\xi})d\xi$
(114)

where $\mathcal{C}$ is the contour from $0$ to $\infty$. However the presence
of singularity in $\chi_{B,i}(\xi)$ introduces complex ambiguities implying
that the integration depends on the choice of contour $\mathcal{C}$. The
choice of the contour actually determine the Stokes coefficients appearing in
the trans-series that encode the full details of the initial conditions.
Appropriate choices of these Stokes coefficients (contour of integration)
correspond to evolution on the attractor surface. In our case, the contours of
integration can be chosen to be the straight lines given by two different
angles for $\chi$ and $\tilde{\chi}$, respectively, as shown in Fig. 9. We
have two Stokes coefficients corresponding to the two exponential decay modes
discussed in Sec. 3.4 and these give the two missing data (other than the
parameters $\alpha$ and $\beta$ appearing in the hydrodynamic expansion) which
specify the initial conditions fully in the four-dimensional phase space.

Figure 9: The top figures shows the ratio of the absolute values of the
$n+1$-th to $n$-th series coefficients of the hydrodynamic expansion of $\chi$
and $\tilde{\chi}$ as function of $n$, where the left panel is for the strong
system and the right panel is for the weak system. The linear growth of this
ratio implies factorial growth of $r_{n}$ for both weak and strong systems.
The bottom figures show the pole structure of both the subsystems on the Borel
plane. The poles of the strong sector and the weak sector are shown in blue
(on left) and red (on right) respectively. It starts from
$\xi_{1,0}=\xi_{2,0}=\frac{3}{2C_{\tau}}\Big{(}\frac{\beta}{\sqrt{\alpha\gamma}}\Big{)}^{1/3}$,
for $\alpha=1.169$ and $\beta=2.747$. The black dotted lines shows the contour
of integration at an angle $\pm\theta$.

##  4 Hydrodynamization on the attractor surface

Figures 5 and 6 show the existence of hydrodynamic tail of the attractor at a
sufficiently late time, giving evidence of hydrodynamization in hybrid
attractor. This preliminary observation leads us to the systematic study of
characteristics of the hydrodynamization times of the two sectors on the
attractor surface.

###  4.1 Hydrodynamization time

Most commonly, hydrodynamization time is defined as the time when both the
longitudinal and transverse pressure can be described by the first-order
hydrodynamics up to certain accuracy Attems:2017zam . For the hybrid
hydrodynamic attractor, we define the hydrodynamization time of each of the
sectors through

$\displaystyle\frac{|\Delta P_{L}|}{P}:=\frac{|\phi-\phi_{1st}|}{P}<0.1,\hskip
11.38092pt\text{for}\hskip 11.38092pt\tau>\tau_{hd},$ (115)

and similarly for the strongly coupled (tilde) sector, so that $\tau_{hd}$ and
$\tilde{\tau}_{hd}$ denote the hydrodynamization times of the perturbative and
non-perturbative sectors, respectively. In general the hydrodynamization time
of perturbative sector $\tau_{hd}$ is larger than that of the holographic
sector $\tilde{\tau}_{hd}$ as shown in Fig 10, where the blue curve
corresponds to the hard sector and the red one to the holographic sector.

Figure 10: The figure shows a doubly logarithmic plot of the hydrodynamization
time for the hard (in red) and the soft sector (in blue) against the total
energy density at the respective $\tau_{hd}$ and $\tilde{\tau}_{hd}$ for ratio
$5:1$. The hydrodynamization time of the hard sector is longer than the soft
sector. The black and the green dots corresponds to the hydrodynamization
times for initial conditions which are physically realizable in the $p$-$p$
and Pb-Pb collision.

Further, if we analyze the hydrodynamization time of the two components of the
hybrid attractor by classifying the total energy density at the
hydrodynamization time into three regions, one can find the following generic
features of them,

First, when the total energy density at the hydrodynamization time of each
sector lies within the region,

$0.0001\gamma^{-1}\lessapprox\mathcal{E}(\tau_{hd}),\mathcal{E}(\tilde{\tau}_{hd})\lessapprox\gamma^{-1}$

(negative slope in the Fig 10). We have found that the hydrodynamization time
of both the sectors follows conformal behaviours i.e.

$\tau_{hd}\sim p\times(4\pi C_{\eta})\times\mathcal{E}(\tau_{hd})^{-1/4}$

and

$\tilde{\tau}_{hd}\sim\tilde{p}\times(4\pi\tilde{C}_{\eta})\times\mathcal{E}(\tilde{\tau}_{hd})^{-1/4}$

where $p$ and $\tilde{p}$ are constants. The value of $p$, $\tilde{p}$ and the
characteristic of the hydrodynamization time in this conformal window is
determined based on which sector dominates at early time. For instance, with
dominance of the perturbative sector at time $\tau_{0}=10^{-3}\gamma^{1/4}$,
we find

$p\sim 0.63$

This is the same value (within numerical accuracy) that is obtained in the
case of conformal attractors. We find that the hydrodynamization time of the
perturbative sector in this conformal region is universal as shown in Fig 11
and the same as in conformal attractors studied earlier in the literature.
However for the holographic sector as in Fig 12, the value of the $\tilde{p}$
increases from $0.64$ to $2.1$ as the ratio of the energy densities of the
individual sectors at some early time $\tau_{0}=0.001\gamma^{1/4}$ increases
from $1:1$ to $1000:1$. Thus the hydrodynamization time of the holographic
sector shows sensitivity towards the initial ratio of the energy densities of
the subsectors. Hence one finds that even in the conformal window the
hydrodynamization time of the soft sector can be substantially larger than the
conformal attractors.

Figure 11: The doubly logarithmic plot of the scaled hydrodynamization time of
the weakly coupled (hard) sector vs the total energy density at the
${\tau}_{hd}$ shows universality and indistinguishability from other conformal
attractors, when the hard sector is dominant at a reference time (here
$\tau_{0}=0.001\gamma^{1/4}$). The various color labels refer to different
ratios of energy densities of the hard to the strongly coupled (soft) sector
at the initial reference time. We note that the hydrodynamization time of the
dominant hard sector is independent of these ratios.  Figure 12: The figure
shows the doubly logarithmic plot of the scaled hydrodynamization time of the
strongly coupled (softer) sector against the total energy density at the
respective $\tilde{\tau}_{hd}$, when the weakly coupled (hard) sector is
dominant at an initial reference time (here $\tau_{0}=0.001\gamma^{1/4}$). The
various color labels refer to different ratios of energy densities of the hard
to the soft sector at the initial reference time. Values of $\tilde{p}$ for
ratios from $1:1$ to $1000:1$ of the initial energy densities of the
subsectors increases from 0.64, 0.76, 0.93,  1.07, 1.21,  1.76 to 2.1. These
different values of $\tilde{p}$ shows the dependency of $\tilde{\tau}_{hd}$ on
the ratio of initial energy densities of the subsectors.

However, with initial conditions where the energy density of the holographic
sector is dominating on the attractor surface, the roles of conformal window
in the subsector reverses i.e., $\tilde{p}$ assumes the conformal value
$0.63$, thus showing universality as in Fig 13 and coinciding with the usual
conformal attractors. Meanwhile from Fig 14, we find that the perturbative
sector shows dependence on the initial energy densities of the subsectors and
the value of the constant $p$ increases with the increase in the initial
energy density of the perturbative sector. This reversal in the behaviour of
hydrodynamization time in the conformal window can be attributed to the fact
that the sub-dominating sector is forced to share its energy with the
dominating sector even before the hydrodynamization, while the dominating
sector is dynamically driven to act universally.

Figure 13: The figure shows the doubly logarithmic plot of the scaled
hydrodynamization time of the strongly coupled (soft) sector against the total
energy density at the respective $\tilde{\tau}_{hd}$ when the soft sector is
dominant at the reference time $\tau_{0}=0.001\gamma^{1/4}.$ The various color
labels refer to different ratios of energy densities of the hard to the soft
sector at the initial reference time. We note that the hydrodynamization time
of the dominant soft sector is independent of these ratios. The plot shows the
universal behaviour of the hydrodynamization time of the dominant sector in
the conformal window, as the value of $\tilde{p}$ is $0.63$, and is the same
as in Fig. 11 corresponding to the cases where the hard sector was dominant at
initial reference time.  Figure 14: The doubly logarithmic plot of the scaled
hydrodynamization time of the weakly coupled (hard) sector vs the total energy
density at the ${\tau}_{hd}$ when the soft sector is dominant at the reference
time $\tau_{0}=0.001\gamma^{1/4}.$ The various color labels refer to different
ratios of energy densities of the hard to the soft sector at the initial
reference time. We note that the hydrodynamization time of the subdominant
hard sector depends on the ratio of the energy densities. Comparing with Fig.
12, where the hard sector was dominant at initial reference time, we observe
that in both of these cases, the more dominant one sector is, the more delayed
is the hydrodynamization of the sub-dominant sector.

Next, when the energy densities at hydrodynamization times of both the sectors
are close to $\gamma^{-1}$ i.e.,

$\mathcal{E}(\tau_{hd}),\mathcal{E}(\tilde{\tau}_{hd})\sim\gamma^{-1}$

the hydrodynamization times on the attractor surface takes minimum values as
shown in Fig 10.

When $\mathcal{E}(\tau_{hd})$ and $\mathcal{E}(\tilde{\tau}_{hd})$ are very
large, the hydrodynamization time increases rapidly and do not show any
scaling behaviour like that in conformal window. This corresponds to the
positive slope region in Fig 10.

###  4.2 Rough phenomenological match to $p$-$p$ and Pb-Pb collisions

In order to contemplate possible phenomenological implications of these
features, we interpret attractor curves with lower total energy densities at
$\tau\sim\gamma^{1/4}$ as corresponding to small systems (high-multiplicity
$p$-$p$ or $A$-$p$ collisions) as opposed to $A$-$A$ heavy-ion collisions (Pb-
Pb or Au-Au). To relate the two, we identify the energy scale of the
intersystem coupling $\gamma^{-1/4}$ with the saturation scale
$Q_{s}^{p\mbox{-}p}$ of $p$-$p$ collisions, and we assume that
$\mathcal{E}(\tau)\propto Q_{s}^{4}$ at $\tau\sim Q_{s}^{-1}$ with large
$A$-$A$ systems having a larger $Q_{s}$. In the IPsat and bCGC models studied
in Kowalski:2007rw  $Q_{s,A}$ at large nucleon number $A$, impact parameter
$b=0$, and Feynman $x$ parameter around 0.001 was found to be given roughly by

$Q_{s,A}^{2}\sim 0.4(x/0.001)^{-0.3}A^{1/3}Q_{s,p}^{2}.$ (116)

For Pb-Pb collisions, where $A^{1/3}\approx 6$, we take simply

$Q_{s,Pb}^{2}=2.5\,Q_{s,p}^{2}$ (117)

as a typical value and get that at $\tau\sim Q_{s,Pb}^{-1}$ the total physical
energy density is $\mathcal{E}_{\mathrm{Pb\mbox{-}Pb}}(\tau\sim 0.63)\sim
6.25\,\mathcal{E}_{p\mbox{-}p}(\tau\sim 1)$ in units where $\gamma=1$.

The tables 1 and 2 show the hydrodynamization time of accordingly defined
$p$-$p$ and Pb-Pb collisions for various ratios of initial energy densities of
the two subsectors.

$\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ | $\varepsilon(\tau_{0}=10^{-3})$ | $\tau_{hd}$ | $\tilde{\tau}_{hd}$ |  $w_{hd}$/10 | $\tilde{w}_{hd}$ | $\mathcal{E}(\tau_{hd})$ | $\mathcal{E}(\tilde{\tau}_{hd})$
---|---|---|---|---|---|---|---
1 | 0.354767 | 13.3212 | 1.95205 | 0.689669 | 1.66915 | 0.0718448 | 0.534596
3 | 0.556843 | 11.9203 | 2.01584 | 0.660307 | 1.7354 | 0.0941532 | 0.549255
4 | 0.63962 | 11.6363 | 2.00203 | 0.654419 | 1.73367 | 0.100039 | 0.562323
5 | 0.63547095 | 11.5021 | 1.96978 | 0.650842 | 1.71212 | 0.102517 | 0.570785
6 | 0.657225 | 11.3631 | 1.93951 | 0.647353 | 1.69386 | 0.105336 | 0.581756
7 | 0.672626 | 11.2338 | 1.91872 | 0.644999 | 1.68478 | 0.108673 | 0.594455
8 | 0.6839 | 11.1518 | 1.88739 | 0.642642 | 1.66384 | 0.11028 | 0.603949
9 | 0.692366 | 11.0969 | 1.85035 | 0.64016 | 1.63598 | 0.110749 | 0.611081
10 | 0.69876 | 11.0462 | 1.82517 | 0.638869 | 1.61918 | 0.11189 | 0.61939
Table 1: Hydrodynamization times for $p$-$p$ initial conditions. Here $w=\mathcal{E}_{1}^{1/4}\tau$, where $\mathcal{E}_{1}$ in (49). and $w_{hd}=w(\tau_{hd})$ and similarly for the second subsystem. $\varepsilon(\tau_{0})/\tilde{\varepsilon}(\tau_{0})$ | $\varepsilon(\tau_{0}=10^{-3})$ | $\tau_{hd}$ | $\tilde{\tau}_{hd}$ |  $w_{hd}$/10 | $\tilde{w}_{hd}$ | $\mathcal{E}(\tau_{hd})$ | $\mathcal{E}(\tilde{\tau}_{hd})$
---|---|---|---|---|---|---|---
1 | 0.355793523 | 10.9847 | 4.67132 | 0.74703 | 3.18484 | 0.40913 | 0.98169
2 | 0.4830497244 | 10.7939 | 4.84945 | 0.77733 | 3.11824 | 0.46223 | 1.03793
3 | 0.555141806 | 10.6434 | 4.86938 | 0.78833 | 3.02682 | 0.49198 | 1.08058
4 | 0.6002567 | 10.5158 | 4.8427 | 0.79347 | 2.93962 | 0.51212 | 1.11510
5 | 0.62995522 | 10.4025 | 4.79626 | 0.79540 | 2.85756 | 0.52760 | 1.145555
8 | 0.674906 | 10.1236 | 4.61724 | 0.79488 | 2.64735 | 0.55847 | 1.22146
9 | 0.68260147 | 10.0498 | 4.55592 | 0.79355 | 2.58616 | 0.56612 | 1.24423
10 | 0.688444 | 9.9797 | 4.4934 | 0.79234 | 2.5295 | 0.57243 | 1.26543
Table 2: Hydrodynamization times for Pb-Pb initial conditions.

Figure 15: The figures show a doubly logarithmic plot of the scaled
hydrodynamization time for the perturbative sector(in the left) and
holographic sector (in the right) against the total energy density at the
respective $\tau_{hd}$ and $\tilde{\tau}_{hd}$. The black and the green dots
are the hydrodynamization times for initial conditions which are physically
realizable in the $p$-$p$ and Pb-Pb collision. The total energy density at
hydrodynamization time for Pb-Pb collisions is scaled by a factor of $1/5$ and
$1/2$ to fit in the same plot with $p$-$p$ collision.

From these tables one finds that in our $p$-$p$ collisons
$\mathcal{E}(\tau_{hd})$ and $\mathcal{E}(\tilde{\tau}_{hd})$ are minimal and
the hydrodynamization time of the perturbative sector is large and that of the
holographic sector is small compare to the Pb-Pb collision. This behaviour is
reflected in Fig 15, where the $\mathcal{E}(\tau_{hd})$ and
$\mathcal{E}(\tilde{\tau}_{hd})$ of the Pb-Pb collision is scaled by the
factor of $1/5$ and $1/2$ to fit in the plot with $p$-$p$ collision. Also one
observes that this phenomenologically realizable hydrodynamization time of the
perturbative sector lies in the conformal region while that of the non-
perturbative sector lies in the non-conformal windows for both $p$-$p$ and Pb-
Pb collisions. These features of the hydrodynamization time of perturbative
and holographic sector for the initial conditions of $p$-$p$ and Pb-Pb
collisions are captured in the Fig 10 as:

  * •

In case of the perturbative sector for both $p$-$p$ and Pb-Pb collisions, the
total energy densities at the hydrodynamization time,
$\mathcal{E}(\tau_{hd})$, lies approximately within the conformal region
denoted by the black and the green dots in the figure. In addition, since the
$\mathcal{E}(\tau_{hd})$ is minimal in $p$-$p$ collisions, the
hydrodynamization time of the perturbative sector $\tau_{hd}$ is sufficiently
larger which is evident from the scaling behaviour of $\tau_{hd}$ in the
conformal window.

  * •

For the non-perturbative sector, $\mathcal{E}(\tilde{\tau}_{hd})>\gamma^{-1}$
and lies outside the conformal region for both $p$-$p$ and Pb-Pb collisions as
shown in the figure. Further due to the small total energy density,
$\mathcal{E}(\tau_{hd})$, in $p$-$p$ collisions, the hydrodynamization time
$\tilde{\tau}_{hd}$ of the soft sector of $p$-$p$ collision is smaller than
that of Pb-Pb collisions, as indicated from the growth of $\tilde{\tau}_{hd}$
with $\mathcal{E}(\tilde{\tau}_{hd})$.

For details on the evolution of energy densities in these scenarios, see
Appendix D.

##  5 Discussion

The main physical insights which we have gained from the Bjorken flow
attractor of our hybrid fluid model is how hydrodynamizaiton works in an
asymptotically free gauge theory with both weakly interacting and strongly
interacting degrees of freedom. Our model has only one dimensionful energy
scale, namely $\gamma^{-1/4}$, which sets the intersystem coupling. We have
shown analytically that at very early time, $\tau\ll\gamma^{1/4}$, the
perturbative sector’s energy density always dominates over that of the
strongly coupled non-perturbative sector with an exponent which can be
determined by the transport coefficients and the relaxation times of the two
sectors.

Furthermore, we have studied hydrodynamization by initializing at time earlier
than $\gamma^{1/4}$. When the total energy densities at the respective
hydrodynamization times ($\mathcal{E}(\tau_{hd})$,
$\mathcal{{E}}(\tilde{\tau}_{hd})$) is less than $\gamma^{-1}$, both sectors
hydrodynamize as in a conformal hydrodynamic attractor
($\tau_{hd}\propto\mathcal{E}(\tau_{hd})^{-1/4}$ and
$\tilde{\tau}_{hd}\propto\mathcal{E}(\tilde{\tau}_{hd})^{-1/4}$) up to a pre-
factor $p$. Remarkably, the dominant sector (which is the perturbative sector
unless we fine-tune the initial conditions) always hydrodynamize in a
universal way, with the pre-factor $p$ same as in a decoupled conformal
hydrodynamic attractor. The inter-system coupling does not affect $p$ for the
dominant sector, however the pre-factor for the sub-dominant (typically
strongly interacting) sector increases with the ratio of energy densities at
the initialization time – more the dominance of the dominant sector, more
delayed is the hydrodynamization of the sub-dominant sector. When the energy
densities at the respective hydrodynamization times increases beyond
$\gamma^{-1}$, the hydrodynamization times of both sectors increase rapidly
after reaching their minimum values.

We also have gained new insights with potential relevance for heavy-ion
collision experiments. Setting initial conditions motivated by saturation
scenarios, we see that the total energy density at the hydrodynamization time
of the strongly coupled sector $\mathcal{{E}}(\tilde{\tau}_{hd})$ lies outside
the conformal window, while that corresponding to the time of
hydrodynamization of the perturbative sector, namely
$\mathcal{{E}}({\tau}_{hd})$, lies within the conformal window. It follows
from our general results that the perturbative sector hydrodynamizes later,
while the non-perturbative strongly interacting sector hydrodynamizes earlier
in $p$-$p$ collisions compared to Pb-Pb collisions, since the energy densities
at the hydrodynamization times are smaller in the case of $p$-$p$ collisions.
Thus we get a new perspective on how collective flow develops in small system
collisions.

Another relevant insight of our model is that although the full system behaves
as a single fluid at late time, its effective equation of state and effective
transport coefficients depend on which curve on the attractor surface that the
system converges to, and is thus process dependent. Although for the equation
of state the process-dependence is sub-leading, the effective shear viscosity
depends on how the full energy is finally shared between the two subsectors.

In the future, we hope to study a fuller semi-holographic set-up with a
kinetic theory coupled to a five-dimensional geometry governed by Einstein’s
equations. This is necessary for confronting our approach with experimental
data, especially for the study of dilepton production Naik:2021yph ;
Coquet:2021lca  and hadronization Giacalone:2019ldn  on the attractor surface.

Another interesting direction of study is to understand how stochastic
statistical and quantum fluctuations, which are suppressed in the large $N$
limit, affect the hydrodynamic attractor. This needs to be understood
separately for both weakly interacting and strongly interacting theories. In
our hybrid set-up, these fluctuations are needed for equi-partition of energy
between the perturbative and holographic sectors and thus for complete
thermalization (we will discuss this more in an upcoming work). In the context
of the attractor, this may lead the Bjorken attractor surface to collapse to a
single curve where the ratio of physical subsystem temperatures goes to one
asymptotically, over a long time-scale. Such fluctuations could be enhanced if
the system passes near a critical point,555See Grossi:2020ezz ; Grossi:2021gqi
; Florio:2021jlx  for studies of effects of stochastic critical fluctuations
on hydrodynamics, Kovtun:2012rj  for a review on stochastic statistical
fluctuations in hydrodynamics and Dore:2020jye  for an investigation of
potential hydrodynamic attractors near the critical point. See also
Mitra:2020hbj ; Rodgers:2022fuv  for studies of superfluid Bjorken flow. and
therefore such a study could be interesting for the search of a critical
endpoint in heavy-ion collisions.

###### Acknowledgements.

The authors would like to thank Michal Heller, Edmond Iancu and Jean-Yves
Ollitrault for useful discussions. AM acknowledges support from the Ramanujan
Fellowship and ECR award of the Department of Science and Technology of India
and also the New Faculty Seed Grant of IIT Madras. AS is supported by the
Austrian Science Fund (FWF), project no. J4406.

##  Appendix A Late time gradient expansion

The coefficients present in the late time expansion of our variables are given
by (note that the subscripts 1 and 2 denote the untilded and tilded sectors,
respectively)

$\displaystyle{\epsilon_{11}}$
$\displaystyle=-\frac{8}{3}C_{\eta}{\epsilon_{10}}^{3/4}\tau_{0}^{-1},\;\;{\epsilon_{12}}=-\frac{2}{9}{\epsilon_{10}}^{1/2}\tau_{0}^{-2}\left(-12C_{\eta}^{2}+4C_{\eta}C_{\tau}+3{\epsilon_{10}}^{1/2}{\epsilon_{20}}\tau_{0}^{2}\gamma\right),$
(118) $\displaystyle{\epsilon_{21}}$
$\displaystyle=-\frac{8}{3}C_{\eta}{\epsilon_{20}}^{3/4}\tau_{0}^{-1},\;\;{\epsilon_{22}}=-\frac{2}{9}{\epsilon_{20}}^{1/2}\tau_{0}^{-2}\left(-12\tilde{C}_{\eta}^{2}+4\tilde{C}_{\eta}\tilde{C}_{\tau}+3{\epsilon_{10}}^{1/2}{\epsilon_{20}}\tau_{0}^{2}\gamma\right),$
(119) $\displaystyle f_{10}$
$\displaystyle=0,\;f_{11}=0,\;f_{12}=0,\;f_{13}=\frac{16}{9}C_{\eta}{\epsilon_{20}}^{3/4}\tau_{0}^{-1},\;f_{14}=-\frac{32}{27}C_{\eta}{\epsilon_{10}}^{1/2}\tau_{0}^{-2}\left(3C_{\eta}-C_{\tau}\right),$
(120) $\displaystyle f_{20}$
$\displaystyle=0,\;f_{21}=0,\;f_{12}=0,\;f_{23}=\frac{16}{9}\tilde{C}_{\eta}{\epsilon_{10}}^{3/4}\tau_{0}^{-2},\;f_{24}=-\frac{32}{27}\tilde{C}_{\eta}{\epsilon_{20}}^{1/2}\tau_{0}^{-2}\left(3\tilde{C}_{\eta}-\tilde{C}_{\tau}\right),$
(121) $\displaystyle A_{11}$
$\displaystyle=-\frac{1}{2}{\epsilon_{20}}\gamma,\;A_{12}=\frac{4}{3}\tilde{C}_{\eta}{\epsilon_{20}}^{3/4}\tau_{0}^{-1}\gamma,$
(122) $\displaystyle A_{21}$
$\displaystyle=-\frac{1}{2}{\epsilon_{10}}\gamma,\;A_{22}=\frac{4}{3}C_{\eta}{\epsilon_{10}}^{3/4}\tau_{0}^{-1}\gamma,$
(123) $\displaystyle B_{11}$
$\displaystyle=\frac{1}{6}{\epsilon_{20}}\gamma,\;B_{12}=0,$ (124)
$\displaystyle B_{21}$
$\displaystyle=\frac{1}{6}{\epsilon_{10}}\gamma,\;B_{22}=0,$ (125)
$\displaystyle C_{11}$
$\displaystyle=\frac{1}{6}{\epsilon_{20}}\gamma\tau_{0},\;C_{12}=-\frac{4}{3}\tilde{C}_{\eta}{\epsilon_{20}}^{3/4}\gamma,$
(126) $\displaystyle C_{21}$
$\displaystyle=\frac{1}{6}{\epsilon_{10}}\gamma\tau_{0},\;C_{22}=-\frac{4}{3}C_{\eta}{\epsilon_{10}}^{3/4}\gamma.$
(127)

The equations of the fluctuations are given by

$\displaystyle\delta\phi^{\prime}(\hat{\tau})$
$\displaystyle=\delta\phi(\hat{\tau})\left(-\frac{{\epsilon_{10}}^{1/4}\tau_{0}}{C_{\tau}}\frac{1}{\hat{\tau}^{1/3}}+\frac{2\left(C_{\eta}-2C_{\tau}\right)}{3C_{\tau}}\frac{1}{\hat{\tau}}+\frac{2\left(C_{\eta}C_{\tau}+3\gamma{\epsilon_{10}}^{1/4}{\epsilon_{20}}\tau_{0}^{2}\right)}{9C_{\tau}{\epsilon_{10}}^{1/4}\tau_{0}}\frac{1}{\hat{\tau}^{5/3}}\right)$
$\displaystyle+\delta\varepsilon(\hat{\tau})\left(\frac{4C_{\eta}}{3C_{\tau}}\frac{1}{\hat{\tau}}-\frac{8C_{\eta}}{27{\epsilon_{10}}^{1/4}\tau_{0}}\frac{1}{\hat{\tau}^{5/3}}\right)$
(128) $\displaystyle\delta\varepsilon^{\prime}(\hat{\tau})$
$\displaystyle=\delta\tilde{\varepsilon}(\hat{\tau})\left(\frac{2\gamma{\epsilon_{10}}}{9}\frac{1}{\hat{\tau}^{7/3}}-\frac{8\gamma
C_{\eta}{\epsilon_{10}}^{3/4}}{9\tau_{0}}\frac{1}{\hat{\tau}^{3}}\right)+\delta\tilde{\phi}(\hat{\tau})\left(-\frac{2\gamma{\epsilon_{10}}}{3}\frac{1}{\hat{\tau}^{7/3}}+\frac{8\gamma
C_{\eta}{\epsilon_{10}}^{3/4}}{3\tau_{0}}\frac{1}{\hat{\tau}^{3}}\right)$
$\displaystyle+\delta\varepsilon(\hat{\tau})\left(-\frac{4}{3\hat{\tau}}+\frac{8\gamma\varepsilon_{20}}{9}\frac{1}{\hat{\tau}^{7/3}}-\frac{32\gamma\tilde{C}_{\eta}{\epsilon_{20}}^{3/4}}{9\tau_{0}}\frac{1}{\hat{\tau}^{3}}\right)+\delta\phi(\hat{\tau})\left(\frac{8\gamma\tilde{C}_{\eta}{\epsilon_{20}}^{3/4}}{3\tau_{0}}\frac{1}{\hat{\tau}^{3}}+\frac{1}{\hat{\tau}}\right)$
(129) $\displaystyle\delta\tilde{\phi}^{\prime}(\hat{\tau})$
$\displaystyle=\delta\tilde{\phi}(\hat{\tau})\left(-\frac{{\epsilon_{20}}^{1/4}\tau_{0}}{\tilde{C}_{\tau}}\frac{1}{\hat{\tau}^{1/3}}+\frac{2(\tilde{C}_{\eta}-2\tilde{C}_{\tau})}{3\tilde{C}_{\tau}}\frac{1}{\hat{\tau}}+\frac{2\left(\tilde{C}_{\eta}\tilde{C}_{\tau}+3\gamma\varepsilon_{20}^{1/4}{\epsilon_{10}}\tau_{0}^{2}\right)}{9\tilde{C}_{\tau}{\epsilon_{20}}^{1/4}\tau_{0}}\hat{\tau}^{1/3}\right)$
$\displaystyle+\delta\tilde{\varepsilon}(\hat{\tau})\left(\frac{4\tilde{C}_{\eta}}{3\tilde{C}_{\tau}}\hat{\tau}-\frac{8\tilde{C}_{\eta}}{27\varepsilon_{20}^{1/4}\tau_{0}}\hat{\tau}^{1/3}\right)$
(130) $\displaystyle\delta\tilde{\varepsilon}^{\prime}(\hat{\tau})$
$\displaystyle=\delta\varepsilon(\hat{\tau})\left(\frac{2\gamma{\epsilon_{20}}}{9}\frac{1}{\hat{\tau}^{7/3}}-\frac{8\gamma\tilde{C}_{\eta}{\epsilon_{20}}^{3/4}}{9\tau_{0}}\frac{1}{\hat{\tau}^{3}}\right)+\delta\phi(\hat{\tau})\left(-\frac{2\gamma{\epsilon_{20}}}{3}\frac{1}{\hat{\tau}^{7/3}}+\frac{8\gamma
C_{\eta}{\epsilon_{20}}^{3/4}}{3\tau_{0}}\frac{1}{\hat{\tau}^{3}}\right)$
$\displaystyle+\delta\tilde{\varepsilon}(\hat{\tau})\left(-\frac{4}{3\hat{\tau}}+\frac{8\gamma{\epsilon_{10}}}{9}\frac{1}{\hat{\tau}^{7/3}}-\frac{32\gamma
C_{\eta}\varepsilon_{10}^{3/4}}{9\tau_{0}}\frac{1}{\hat{\tau}^{3}}\right)+\delta\tilde{\phi}(\hat{\tau})\left(\frac{8\gamma
C_{\eta}{\epsilon_{10}}^{3/4}}{3\tau_{0}}\frac{1}{\hat{\tau}^{3}}+\frac{1}{\hat{\tau}}\right)$
(131)

##  Appendix B Phase transition

As was explored in the inviscid case in Kurkela:2018dku , the coupling
equations as presented exhibit a phase transition depending on the value of
$r$. Thermal equilibrium is best explored in the Minkowski background. For the
purposes of this section

$\displaystyle\varepsilon=3P,\quad
P=nT^{4}_{1},\quad\text{and}\quad\mathcal{T}=aT_{1}=\tilde{a}T_{2},$ (132)

where the last relation is the requirement of thermal equilibrium, which means
that the physical system temperature, $\mathcal{T}$, can be expressed in terms
of the individual subsystem temperature via

$\displaystyle\sqrt{-\eta_{00}}\mathcal{T}=\sqrt{-g_{00}}T_{1}=\sqrt{-\tilde{g}_{00}}T_{2}.$
(133)

Further, we introduce the light cone velocity

$\displaystyle v=a/b,\quad\text{and}\quad\tilde{v}=\tilde{a}/\tilde{b},$ (134)

which represents the effective causal light cone of each individual subsystem.
Interactions between the two systems change the effective light cone velocity.
Note that $c>v,\tilde{v}>0$, where $c=1$ is the speed of light.

We first outline some results in the case of identical subsystems. It is
instructive to understand the phase transition by considering the inviscid
Minkowski equations, for details see Kurkela:2018dku .

For identical inviscid conformal subsystems in thermal equilibrium, we have
that the Minkowski coupling equations can be solved until there is one left,
which is given by Kurkela:2018dku

$\displaystyle\gamma
n_{1}\mathcal{T}^{4}=\frac{v^{5}(1-v^{2})(3+v^{2})}{(3+v^{4}-3r(1-v^{2})^{2})^{2}}.$
(135)

Requiring that $0<\mathcal{T}<\infty$ and remembering that $0<v<1,$ we see
that $r>1$. Analytically, a phase transition was found in Kurkela:2018dku
(see Appendix D for more details) when the coupling takes the value

$\displaystyle
r_{c}=\frac{1}{540}\left(195+43\sqrt{15}+\sqrt{30(4082-557\sqrt{15})}\right)\approx
1.114.$ (136)

We can explore this phase transition in the Bjorken case of identical
subsystems by mapping from Bjorken to Minkowski. The identical subsystem
equations for the Bjorken case are given by (18) and (25)-(27), after
identifying the tilded quantities with the untilded ones, i.e.
$\tilde{a}\rightarrow a$.

First, as a sanity check, we do the matching for the inviscid case. We show an
illustrative example with $r=2$ in the left panel of Fig. 16. In the Minkowski
case, we solve (135) directly and arrive at $v=v(\mathcal{T})$. In the Bjorken
case, we parametrically plot the light cone velocity $v(\tau)$ as a function
of temperature $\mathcal{T}=a(\tau)(\varepsilon(\tau)/3)^{1/4}$. Early time is
represented by the green dot, with late time coinciding with the red dot. The
combined choice of initial energy density and initialization time sets the
value of the green point, which in the case of the plot we take
$\gamma\varepsilon(\tau_{0}=0.0001)=0.348$.

Likewise, in the right panel of Fig. 16 we take a value of the coupling below
the critical value of $r$, namely $r=1.01<r_{c}$. In this case, we see the
presence of the phase transition via the multivaluedness of the light cone
velocity as a function of $\mathcal{T}$. Of course, the light cone velocity
can only take one value at a time, which means there is a discontinuous jump
around the critical temperature. As expected, we see that the there is
complete agreement between the two sets of equations.

Next, we consider adding viscous corrections. We now have to contend with the
evolution of $\phi$ given by (20) with values of the transport coefficients
$\eta=4C_{\eta}\varepsilon(\tau)^{3/4}/3$ and
$\tau_{\pi}=C_{\tau_{\pi}}\varepsilon^{-1/4}$ with $C_{\eta}=1/4\pi$ and
$C_{\tau_{\pi}}=(2-\log{2})/2\pi.$ The matching with the inviscid Minkowski
case will not be perfect since we are no longer in thermal equilibrium, but it
will serve as a guide. We find that there is little difference to the above
parametric plots for late times.

Now we go to the full case considered in this paper: unequal subsystems. To
study this numerically, we consider equal initial conditions for both
subsystems. As we see in Fig. 17, due to the choice of the initial conditions,
the both subsystems have initially the same light cone velocities. Then for
intermediate times, the individual subsystems have a markedly different
evolution at intermediate times before the light cone velocities approach the
vacuum value of 1. Interestingly, we can determine situations where the light
cone velocity in one system is multivalued, while in the other it is not. Of
course, it is worthwhile to repeat that we are not in thermal equilibrium,
which means we are not exactly mapping to the physical temperature
$\mathcal{T}$, but instead to an effective temperature.

Figure 16: Light cone velocity as a function of temperature for the inviscid
Bjorken and Minkowski cases of identical subsystems for two different values
of $r$. The orange dashed curve is the result of working directly in the
Minkowski background (135), as studied in Kurkela:2018dku . The blue line
denotes the Bjorken evolution with increasing $\tau$ going from the green to
red dot (right to left on the plot). Left: in the case $r=2$, the green dot is
for $\gamma\varepsilon(\tau_{0}=0.0001)=0.348$, while the red dot denotes
$\tau=20,$ the time we choose to stop the simulation. Right: below the phase
transition, for $r=1.01$ the Bjorken system is initialized at the green dot
with $\gamma\varepsilon(\tau_{0}=0.01)=0.1$.

Figure 17: Light cone velocity as a function of temperature for the viscous
Bjorken and Minkowski cases for a variety of $r$. The orange dashed curve is
the inviscid case in the Minkowski background as shown in Fig. 16. As before,
the evolution in the $\tau$ coordinate follows from the green to red dot
(right to left). The upper solid purple line denotes the more viscous system,
while the lower solid blue line is the weakly coupled system. The initial
conditions are identical for both systems for each $r$ and with initially no
dissipation $\phi(\tau_{0})=\tilde{\phi}(\tau_{0})=0$ at $\tau_{0}=1$, i.e.
top: $r=2$ and $\gamma\varepsilon=\gamma\tilde{\varepsilon}=0.348$, middle:
$r=r_{c}$ with $\gamma\varepsilon=\gamma\tilde{\varepsilon}=0.248$, and
finally bottom: $r=1.2$ with
$\gamma\varepsilon=\gamma\tilde{\varepsilon}=0.3$, while the transport
coefficients are given by (30).

##  Appendix C The disperser surface

On the disperser surface, we have the following scaling behaviors near
$\tau=0$:

$\displaystyle\gamma\tilde{\varepsilon}$
$\displaystyle\rightarrow\sqrt{\frac{r-1}{r}},\gamma\varepsilon\rightarrow{g_{1}}\tau^{\frac{4}{3}(\sigma_{2}-\sigma_{1})},$
(137) $\displaystyle\chi$
$\displaystyle\rightarrow-\sigma_{1}+k_{1}\tau^{\frac{1}{3}(2-\sigma_{2})},\tilde{\chi}\rightarrow-\sigma_{2}+k_{2}\tau^{-\frac{5}{3}\sigma_{1}+\frac{4}{3}\sigma_{2}+\frac{2}{3}},$
(138) $\displaystyle a$ $\displaystyle\rightarrow
a_{10}\tau^{\sigma_{2}-\frac{4}{3}\sigma_{1}-\frac{1}{3}},b\rightarrow
b_{10}\tau^{-\frac{1}{3}(\sigma_{2}+1)},c\rightarrow
c_{10}\tau^{\frac{1}{3}(2-\sigma_{2})},$ (139) $\displaystyle\tilde{a}$
$\displaystyle\rightarrow
a_{20}\tau^{-\frac{1}{3}(\sigma_{2}+1)},\tilde{b}\rightarrow
b_{20}\tau^{-\frac{1}{3}(\sigma_{2}+1)},\tilde{c}\rightarrow
c_{20}\tau^{\frac{1}{3}(2-\sigma_{2})}.$ (140)

We note that $\chi$ and $\tilde{\chi}$ go to $-\sigma_{1}$ and $-\sigma_{2}$
respectively, $\gamma\varepsilon$ vanishes while $\gamma\tilde{\varepsilon}$
goes to $\sqrt{r-1}/\sqrt{r}$ as $\tau\rightarrow 0$. The above scalings
follow on the disperser surface if instead of (43), we require that
$(4/3)\sigma_{1}<\sigma_{2}<1$ with $0<\sigma_{1}<1$ which is satisfied by our
choices of parameters.666More specifically, we get these scalings if
$a\varepsilon/\tilde{a}\rightarrow 0$ as $\tau\rightarrow 0$.

The five identities analogous to those in Eq. (47) which should hold at
$\tau=0$ are:

$\displaystyle\frac{\tilde{b}}{\tau^{-\frac{1}{3}(\sigma_{2}+1)}}=\frac{\tilde{c}}{\tau^{\frac{1}{3}(2-\sigma_{2})}}=\sqrt{\frac{r}{r-1}}\frac{\tilde{a}}{\tau^{-\frac{1}{3}(\sigma_{2}+1)}}$
(141)
$\displaystyle=\sqrt{\frac{3\sqrt{r(r-1)}}{\gamma\tilde{\varepsilon}(4r-1-2(r-1)\sigma_{2})}}\frac{b}{\tau^{-\frac{1}{3}(\sigma_{2}+1)}}=\sqrt{\frac{3\sqrt{r(r-1)}}{\gamma\tilde{\varepsilon}(4r-1+4(r-1)\sigma_{2})}}\frac{c}{\tau^{\frac{1}{3}(2-\sigma_{2})}}$
$\displaystyle=\frac{b}{\tau^{-\frac{1}{3}(\sigma_{2}+1)}}\sqrt{r\gamma\varepsilon\frac{c}{\tau
a_{1}}}$

We also find that the other six coefficients of the early time expansion can
be determined in terms of $k_{1}$ and $g_{1}$ in consistency with the fact
that the disperser surface is two-dimensional.

Figure 18: A disperser solution is studied with initial conditions
$\varepsilon(\tau_{0})=0.1$, $\tilde{\varepsilon}(\tau_{0})=0.3$,
$\phi_{1}(\tau_{0})=0.02$ and $\phi_{2}(\tau_{0})=0.02$ at $\tau_{0}=1$ with
all parameters same as in the case of our attractor solutions except $r=3$.
Left panel: In the limit, $\tau\rightarrow 0$, $\chi$ and $\tilde{\chi}$ go to
$-\sigma_{1}$ and $-\sigma_{2}$ respectively. Right panel: The opposite of
bottom-up happens: $\gamma\tilde{\varepsilon}$ goes to a constant
$\sqrt{(r-1)/r}$ and $\gamma\epsilon$ goes to zero as
$\tau^{\frac{4}{3}(\sigma_{2}-\sigma_{1})}$. Figure 19: A check of the grand
identity (141) for the disperser solution in Fig. 18.

In Fig. 18, the early time behaviors of $\chi$, $\tilde{\chi}$, $\varepsilon$
and $\tilde{\varepsilon}$ in a disperser solution have been shown, and in Fig.
19, we have shown the numerical verification of (141) to a very high accuracy.

##  Appendix D Different initial ratios of subsector energy densities

While at sufficiently early times all attractor solutions have a dominant
weakly coupled component, different initial conditions lead to different
ratios of energy densities at some initial reference time, chosen here as
$\tau_{0}=0.001\gamma^{1/4}$. In figures 20 and 20 we display the evolution of
the energy densities in the subsectors and the total energy density for
different initial ratios, comparing the two sets of initial conditions we
tentatively associated with $p$-$p$ and Pb-Pb collisions in sect. 4.2.

Figure 20: Energy densities for $p$-$p$ collisions (left column) and Pb-Pb
(right column) for a variety of initial ratios
$\varepsilon/\tilde{\varepsilon}$. For both $p$-$p$ and Pb-Pb collisions, as
the ratio of the initial energy densities of the two sectors increases, the
energy of the hard sector dominates for longer time.

For a ratio $\mathcal{E}_{1}:\mathcal{E}_{2}$ less than 4:1, the energy
density in the holographic sector overtakes the one in the perturbative sector
for a stretch of time, whereas in the opposite case it never gains dominance.
At sufficiently large time, the “perturbative” sector becomes again dominant
in each case, but although this is reminiscent of the transition to a weakly
coupled hadron gas after hadronization, we consider our model as a toy of
heavy-ion collisions only in the deconfined early stages.

## References

  * (1) P. M. Chesler and L. G. Yaffe, _Boost invariant flow, black hole formation, and far-from-equilibrium dynamics in N = 4 supersymmetric Yang-Mills theory_ , _Phys. Rev. D_ 82 (2010) 026006, [0906.4426].
  * (2) P. M. Chesler and L. G. Yaffe, _Holography and colliding gravitational shock waves in asymptotically AdS 5 spacetime_, _Phys. Rev. Lett._ 106 (2011) 021601, [1011.3562].
  * (3) 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].
  * (4) P. M. Chesler, _Colliding shock waves and hydrodynamics in small systems_ , _Phys. Rev. Lett._ 115 (2015) 241602, [1506.02209].
  * (5) M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana et al., _Holographic Collisions in Non-conformal Theories_ , _JHEP_ 01 (2017) 026, [1604.06439].
  * (6) M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana et al., _Paths to equilibrium in non-conformal collisions_ , _JHEP_ 06 (2017) 154, [1703.09681].
  * (7) P. Romatschke, _Relativistic Fluid Dynamics Far From Local Equilibrium_ , _Phys. Rev. Lett._ 120 (2018) 012301, [1704.08699].
  * (8) 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].
  * (9) L. Keegan, A. Kurkela, P. Romatschke, W. van der Schee and Y. Zhu, _Weak and strong coupling equilibration in nonabelian gauge theories_ , _JHEP_ 04 (2016) 031, [1512.05347].
  * (10) 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].
  * (11) M. P. Heller, R. A. Janik and P. Witaszczyk, _Hydrodynamic Gradient Expansion in Gauge Theory Plasmas_ , _Phys. Rev. Lett._ 110 (2013) 211602, [1302.0697].
  * (12) M. P. Heller and M. Spalinski, _Hydrodynamics Beyond the Gradient Expansion: Resurgence and Resummation_ , _Phys. Rev. Lett._ 115 (2015) 072501, [1503.07514].
  * (13) M. P. Heller, _Holography, Hydrodynamization and Heavy-Ion Collisions_ , _Acta Phys. Polon. B_ 47 (2016) 2581, [1610.02023].
  * (14) M. Spaliński, _On the hydrodynamic attractor of Yang–Mills plasma_ , _Phys. Lett. B_ 776 (2018) 468–472, [1708.01921].
  * (15) P. Romatschke, _Relativistic Hydrodynamic Attractors with Broken Symmetries: Non-Conformal and Non-Homogeneous_ , _JHEP_ 12 (2017) 079, [1710.03234].
  * (16) 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].
  * (17) 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].
  * (18) J. Casalderrey-Solana, N. I. Gushterov and B. Meiring, _Resurgence and Hydrodynamic Attractors in Gauss-Bonnet Holography_ , _JHEP_ 04 (2018) 042, [1712.02772].
  * (19) 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].
  * (20) M. Strickland and U. Tantary, _Exact solution for the non-equilibrium attractor in number-conserving relaxation time approximation_ , _JHEP_ 10 (2019) 069, [1903.03145].
  * (21) S. Jaiswal, C. Chattopadhyay, A. Jaiswal, S. Pal and U. Heinz, _Exact solutions and attractors of higher-order viscous fluid dynamics for Bjorken flow_ , _Phys. Rev. C_ 100 (2019) 034901, [1907.07965].
  * (22) 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].
  * (23) 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].
  * (24) A. Dash and V. Roy, _Hydrodynamic attractors for Gubser flow_ , _Phys. Lett. B_ 806 (2020) 135481, [2001.10756].
  * (25) M. P. Heller, R. Jefferson, M. Spaliński and V. Svensson, _Hydrodynamic Attractors in Phase Space_ , _Phys. Rev. Lett._ 125 (2020) 132301, [2003.07368].
  * (26) J.-P. Blaizot and L. Yan, _Analytical attractor for Bjorken flows_ , _Phys. Lett. B_ 820 (2021) 136478, [2006.08815].
  * (27) D. Almaalol, A. Kurkela and M. Strickland, _Nonequilibrium Attractor in High-Temperature QCD Plasmas_ , _Phys. Rev. Lett._ 125 (2020) 122302, [2004.05195].
  * (28) V. E. Ambrus, S. Busuioc, J. A. Fotakis, K. Gallmeister and C. Greiner, _Bjorken flow attractors with transverse dynamics_ , _Phys. Rev. D_ 104 (2021) 094022, [2102.11785].
  * (29) Z. Du, X.-G. Huang and H. Taya, _Hydrodynamic attractor in a Hubble expansion_ , _Phys. Rev. D_ 104 (2021) 056022, [2104.12534].
  * (30) J.-P. Blaizot and L. Yan, _Attractor and fixed points in Bjorken flows_ , _Phys. Rev. C_ 104 (2021) 055201, [2106.10508].
  * (31) C. Chattopadhyay, S. Jaiswal, L. Du, U. Heinz and S. Pal, _Non-conformal attractor in boost-invariant plasmas_ , _Phys. Lett. B_ 824 (2022) 136820, [2107.05500].
  * (32) L. J. Naik, S. Jaiswal, K. Sreelakshmi, A. Jaiswal and V. Sreekanth, _Hydrodynamical attractor and thermal particle production in heavy-ion collision_ , 2107.08791.
  * (33) Z. Chen and L. Yan, _Hydrodynamic attractor in the nonconformal Bjorken flow_ , _Phys. Rev. C_ 105 (2022) 024910, [2109.06658].