\section{Introduction}
\label{intro}

One of the most remarkable discoveries in non-equilibrium dynamics recently has been that quantum many-body systems can \textit{hydrodynamize} far away from equilibrium \cite{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 \cite{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 \cite{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 \cite{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.\ \cite{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 \cite{Gelis:2010nm} and subsequently kinetic theory \cite{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. \ref{sec:set-up}, 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 \cite{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{\"u}ller-Israel-Stewart (MIS) \cite{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.\footnote{At 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 \cite{Ecker:2018ucc}.}

The existence of hydrodynamic attractors in such hybrid fluid models has been established in \cite{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.~\cite{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 \cite{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
\cite{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) \cite{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 \cite{Kowalski:2007rw}. Generically, we find:
\begin{itemize}
    \item The total energy densities at hydrodynamization time of the perturbative sector ($\mathcal{E}({\tau}_{hd})$) are \textit{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 \textit{larger} in $p$-$p$ collisions as readily follows from the scaling behavior of $\tau_{hd}$ in the conformal window.
    \item The total energy densities at hydrodynamization time of the holographic sector ($\mathcal{E}(\widetilde{\tau}_{hd})$) are \textit{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 \textit{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})$.
\end{itemize}
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.~\ref{sec:set-up}, we review our set-up of two fluids undergoing Bjorken flow and coupled mutually via the democratic effective metric coupling \cite{Banerjee:2017ozx,Kurkela:2018dku}. We further discuss the method to discern the attractor surface in this hybrid system. In Sec.~\ref{sec:early}, 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 \cite{Mitra:2020mei}. Finally in Sec.~\ref{sec:hydrodynamization}, we systematically study the hydrodynamization on the attractor surface and connect to the phenomenology of $p$-$p$ and Pb-Pb collisions.

\section{Set-up}\label{sec:set-up}

\subsection{The hybrid action}

A minimalist set-up required to obtain our hybrid fluid model involves only two leading \textit{democratic effective metric couplings} \cite{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 \cite{Kurkela:2018dku}
\begin{align}\label{Eq:FullAction}
S[\psi,\tilde{\psi},G_{\mu\nu}]&=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}],\nonumber\\
S_{\rm int}&=\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}\nonumber\\
&+\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).
\end{align}
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
\begin{align}\label{ward}
    \nabla_\mu t^{\mu\nu}=0\quad \text{and}\quad \tilde{\nabla}_\mu \tilde{t}^{\mu\nu}=0,
\end{align}
where the energy momentum tensors are
\begin{align}
    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}},
\end{align}
and $\nabla$ and $\tilde{\nabla}$ are built out of ${g}$ and ${\tilde{g}}$, respectively.

The interaction term $S_{int}$ in \eqref{Eq:FullAction} 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
\begin{equation}\label{Eq:coup-var}
    \frac{\delta S}{\delta {g}_{\mu\nu}} =0, \quad \frac{\delta S}{\delta {g}_{\mu\nu}} =0.
\end{equation}
The above yield the following \textit{coupling equations}
\begin{align}\label{eq-coup}
{g}_{\mu\nu}&=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],\nonumber\\
{\tilde{g}}_{\mu\nu}&=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],
\end{align}
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 \eqref{eq-coup} 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 \eqref{Eq:FullAction} with respect to the physical background metric $G_{\mu\nu}$, which yields
\begin{align}
    T^\mu_{\phantom{\mu}\nu}
    &=\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,
    \nonumber\\
\Delta K&=-\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].
\label{fullemt}\end{align}
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 \eqref{eq-coup} and the two Ward identities \eqref{ward} indeed imply the local conservation of the full energy-momentum tensor in the physical background metric, i.e.
\begin{align}
    \nabla_\mu^{(B)}  T^\mu_{\phantom{\mu}\nu}=0.
\end{align}

The mutual effective metric coupling thus leads to a full energy-momentum tensor \eqref{fullemt} 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,\footnote{In \cite{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 \cite{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 \eqref{eq-coup} using the respective subsystem stress tensors obtained from the previous iteration. The subsystems are solved once again with the \textit{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 \cite{Iancu:2014ava}. Subsequent works \cite{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 \cite{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.

\subsection{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 \eqref{ward} and the algebraic coupling equations \eqref{eq-coup} 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
\begin{align}
t^{\mu\nu} &=(\varepsilon+P)u^\mu u^\nu+P {g}^{\mu\nu}+\Pi^{\mu\nu},\\
\tilde{t}^{\mu\nu} &=(\tilde{\varepsilon}+\tilde{P})\tilde{u}^\mu \tilde{u}^\nu+\tilde{P} {\tilde{g}}^{\mu\nu}+\tilde{\Pi}^{\mu\nu},
\end{align}
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 \eqref{ward}.
The four-velocities are properly normalized with respect to their metrics, i.e.
\begin{align}
u^\mu {g}_{\mu\nu} u^\nu =-1 \quad \text{and} \quad \tilde{u}^\mu {\tilde{g}}_{\mu\nu} \tilde{u}^\nu =-1.
\end{align}
The dissipation tensors, $\Pi^{\mu\nu}$ and $\tilde{\Pi}^{\mu\nu}$, are taken to be orthogonal to their respective four velocities
\begin{align}
\Pi^{\mu\nu} u_\mu=0 \quad \text{and} \quad \tilde{\Pi}^{\mu\nu} \tilde{u}_\mu=0.
\end{align}
The two fluids are assumed to have dissipative dynamics dictated by M{{\"u}}ller-Israel-Stewart (MIS) theory \cite{Muller:1967zza,Israel:1979wp},\footnote{The 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 \cite{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
\begin{align}
(\tpu^\mu \nabla_\mu+1)\Pi^{\mu\nu}&=-\eta \sigma^{\mu\nu},\nonumber\\
(\tilde{\tau}_\pi\tilde{u}^\mu \tilde{\nabla}_\mu+1)\tilde{\Pi}^{\mu\nu}&=-\tilde{\eta} \tilde{\sigma}^{\mu\nu},
\label{mis}
\end{align}
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 \eqref{fullemt} 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) \cite{Kurkela:2018dku}.

\subsection{Explicit equations for boost-invariant Bjorken flow}\label{sec:explicit}
We now provide the explicit form of the complete set of equations, \eqref{ward}, \eqref{eq-coup} and \eqref{mis}. 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
\begin{align}\label{bjorken-metric}
G_{\mu\nu}&=\text{diag}(-1,1,1,\tau^2).
\end{align}
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:
\begin{align}
{g}_{\mu\nu}(\tau)=\text{diag}(-a(\tau)^2,b(\tau)^2,b(\tau)^2,c(\tau)^2),\nonumber\\
{\tilde{g}}_{\mu\nu}(\tau)=\text{diag}(-\tilde{a}(\tau)^2,\tilde{b}(\tau)^2,\tilde{b}(\tau)^2,\tilde{c}(\tau)^2).\label{coup}
\end{align}
The four velocities are then
\begin{align}
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).
\end{align}
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 \cite{Heller:2015dha}, which means that
\begin{align}
 \Pi^{\mu\nu} &= \text{diag}\left(0,\frac{\phi(\tau)}{2},\frac{\phi(\tau)}{2},\frac{-\phi(\tau)}{\tau^2}\right) ,\\
\tilde{  \Pi}^{\mu\nu} &= \text{diag}\left(0,\frac{\tilde{\phi}(\tau)}{2},\frac{\tilde{\phi}(\tau)}{2},\frac{-\tilde{\phi}(\tau)}{\tau^2}\right).
 \end{align}
Then explicitly, the only non-vanishing components of the conservation equations \eqref{ward} are
\begin{align}\label{exp-ward}
0&=\partial_\tau \varepsilon+\frac{4}{3}\varepsilon \partial_\tau \log(b^{2}c)+\phi \partial_\tau\log(\frac{b}{c}),\\ \label{exp-wardt}
0&=\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}}),
\end{align}
the MIS equations are \eqref{mis}
\begin{align}\label{exp-mis}
0&=\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,\\
0&=\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}, \label{exp-mist}
\end{align}
and finally, the algebraic coupling equations \eqref{coup} are
\begin{align}\label{coup1}
1-a^2&=\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{]},\\ \label{coup2}
b^2-1&=\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{]},\\ \label{coup3}
c^2-\tau ^2&=\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{]},\\ \label{coup4}
1-\tilde{a}^2&=\gamma\frac{a b^2 c}{\tau} \Big{[} \frac{\varepsilon}{a^2}+ r (-\frac{\varepsilon}{a^2}+\frac{\phi}{b^2}+\frac{2 \varepsilon}{3 b^2}+\tau
   ^2 (\frac{\varepsilon}{3 c^2}-\frac{\phi}{c^2}))\Big{]},
   \\\label{coup5}
\tilde{b}^2-1&=\gamma\frac{a b^2 c}{\tau} \Big{[}\frac{  (\phi+\frac{2 \varepsilon}{3})}{2 b^2}-  r (-\frac{\varepsilon}{a^2}+\frac{\phi}{b^2}+\frac{2 \varepsilon}{3 b^2}+\tau ^2 (\frac{\varepsilon}{3 c^2}-\frac{\phi}{c^2}))\Big{]},
\\ \label{coup6}
\tilde{c}^2-\tau ^2&=\gamma\frac{a b^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}{3 b^2}+\tau ^2 (\frac{\varepsilon}{3 c^2}-\frac{\phi}{c^2}))\Big{]},
   \end{align}
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 \eqref{bjorken-metric} 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 \eqref{exp-ward}, \eqref{exp-wardt}, \eqref{exp-mis} and \eqref{exp-mist}. 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  \eqref{coup1}-\eqref{coup6}. 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
\begin{align}
\eta=C_\eta \frac{\varepsilon+P}{\varepsilon^{1/4}}, \quad &\text{and} \quad \tau_\pi=C_\tau \varepsilon^{-1/4},\\
\tilde\eta=\tilde{C}_\eta \frac{\tilde{\varepsilon}+\tilde{P}}{\tilde{\varepsilon}^{1/4}}, \quad &\text{and} \quad \tilde{\tau}_\pi=\tilde{C}_\tau \tilde{\varepsilon}^{-1/4}.
\end{align}
In the simulations that follow, we will usually choose
\begin{align}\label{parameters}
\tilde{C}_\eta=\frac{1}{4\pi}, \quad \tilde{C}_\tau=\frac{2-\log2}{2\pi},\quad
C_\eta=10 \tilde{C}_\eta,  \quad \text{and} \quad C_\tau=5 C_\eta.
\end{align}
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
\begin{align}\label{parameters2}
 r\equiv-\gamma'/\gamma=2
\end{align}
unless noted otherwise.
In the inviscid Minkowski case \cite{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~\ref{sec:phase}.
{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
\begin{align}
\frac{4}{3} \chi \left(\frac{b'}{b}-\frac{c'}{c}\right)+\frac{4}{3} \left(\frac{2
   b'}{b}+\frac{c'}{c}\right)+\frac{\varepsilon'}{\varepsilon} &= 0,\\
   \frac{4}{3} \tilde{\chi} \left(\frac{\tilde{b}'}{\tilde{b}}-\frac{\tilde{c}'}{\tilde{c}}\right)+\frac{4}{3} \left(\frac{2
   \tilde{b}'}{\tilde{b}}+\frac{\tilde{c}'}{\tilde{c}}\right)+\frac{\tilde{\varepsilon}'}{\tilde{\varepsilon}} &= 0,\\
-\frac{a \chi {\varepsilon}^{1/4}}{C_{\text{$\tau
   $}}}-\Big(\frac{4 C_{\text{$\eta $}}}{3  C_{\text{$\tau
   $}}}-\frac{4
   \chi{}^2}{3}\Big)\Big(\frac{b'}{b}-\frac{c'}{c} \Big)-\chi'&= 0,\\
-\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}'}{\tilde{b}}-\frac{\tilde{c}'}{\tilde{c}} \Big)-\tilde{\chi}'&= 0,
\end{align}
where the prime denotes a derivative with respect to the proper time, e.g. $\varepsilon^\prime=\partial_\tau \varepsilon.$
The full energy momentum tensor \eqref{fullemt} is given by
\begin{align}
    \mathcal{T}^{\mu\nu}=\text{diag}(\mathcal{E},P_T,P_T,P_L),
\end{align}
where the total energy density, total transverse pressure and longitudinal pressure are
\begin{align}\label{total-energy}
   & \mathcal{E}=-\frac{a b^2 c \varepsilon +\tilde{a} \tilde{b}^2 \tilde{c} \tilde{\varepsilon}}{\tau }+\frac{ \gamma}{36 a \tilde{a} c \tilde{c} \tau^2} \Big(a^2 \tilde{a}^2 \Big[2 b^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)\nonumber\\
   &+c^2 (2 \varepsilon +3 \phi ) \left(2
   \tilde{b}^2 r \tau ^2 (\tilde{\varepsilon}-3 \tilde{\phi})+\tilde{c}^2 (2 r-1) (2
   \tilde{\varepsilon}+3 \tilde{\phi})\right)\Big]\nonumber\\
   &-6 a^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)\nonumber
   \\
   &-6 b^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),\\
&P_T=\frac{a b^2 c (2 \varepsilon +3 \phi )+\tilde{a} \tilde{b}^2 \tilde{c} (2 \tilde{\varepsilon}+3 \tilde{\phi})}{6 \tau }\nonumber\\
    &+\frac{\gamma }{36 ac \tilde{a} \tilde{c} \tau ^2} \Big[\Ca \tilde{a}^2 \Big(2 b^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)\nonumber\\
   &+c^3 (2 \varepsilon +3 \phi ) \left(2
   \tilde{b}^2 r \tau ^2 (\tilde{\varepsilon}-3 \tilde{\phi})+\tilde{c}^2 (2 r-1) (2
   \tilde{\varepsilon}+3 \tilde{\phi})\right)\Big)\nonumber\\
   &-6 b^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)\nonumber\\
   &-6 a c \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],\\
&P_L=P_T-\frac{3 \left(a b^2 c \phi +\tilde{a} \tilde{b}^2 \tilde{c} \tilde{\phi}\right)}{2
   \tau }.
\end{align}

\subsection{Determining the attractor surface}\label{sec:detatt}

 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 \cite{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,
 \begin{equation}\label{Eq:ChChtET}
     \lim_{\tau\rightarrow0}\chi = \sigma_1, \quad \lim_{\tau\rightarrow0}\tilde\chi = \sigma_2,
 \end{equation}
where
\begin{align}\label{Eq:s1s2def}
    \sigma_1 = \sqrt{\frac{C_{\eta}}{C_{\tau}}}, \quad
\sigma_2 = \sqrt{\frac{\tilde{C}_\eta}{\tilde{C}_\tau}}.
\end{align}
See \cite{Denicol:2017lxn} for similar analytic behavior of the attractor.
{In Fig.~\ref{fig:attractor} we show the evolution of the anisotropies in the different sectors and the total system for one particular attractor solution
and neighboring trajectories (from \cite{Mitra:2020mei}, where more plots are given for this particular solution).}

\begin{figure}
\center
\includegraphics[width=.75\columnwidth]{chisfullthin}%
\caption{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 \cite{Mitra:2020mei}, where more plots are given for this particular solution.
\label{fig:attractor}}
\end{figure}

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 \eqref{coup1}-\eqref{coup6} at the initial time to determine the metric components at the initial time. Subsequently, we evolve the system backwards in time via \eqref{exp-ward}, \eqref{exp-wardt}, \eqref{exp-mis} and \eqref{exp-mist}, 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 \eqref{Eq:ChChtET}. 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 \eqref{parameters} and \eqref{parameters2}, 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
\begin{align}
   \gamma \varepsilon(\tau_0) &= \gamma\tilde{\varepsilon}(\tau_0) = 0.18, \nonumber \\\gamma\phi_1 (\tau_0) &= 0.10002105 \quad \text{and}\quad \gamma\phi_2 (\tau_0) = 0.0501539602.
\end{align} and we evolve the system backward  in time.
\begin{figure}
\begin{minipage}{0.49\textwidth}
\includegraphics[width=\linewidth]{chiplotsToni1.pdf}
\end{minipage}
\begin{minipage}{0.49\textwidth}
\includegraphics[width=\linewidth]{epsilonplotspaper1.pdf}
\end{minipage}
\caption{{Illustrative behaviour of an attractor curve -- see text in Sec. \ref{sec:detatt} for the choice of parameters. Top: $\chi$ and $\tilde{\chi}$ approach $\sigma_1$ and $\sigma_2$ respectively as $\tau\rightarrow0$.
Bottom: $\varepsilon$ and $\tilde{\varepsilon}$ satisfy the scaling laws \eqref{Eq:scaling1att1} at early time. We also see that $\gamma \varepsilon$ tends towards $\sqrt{(r-1)/r}=1/\sqrt{2}$ (orange dotted line) as claimed in \eqref{Eq:scaling1att1}.} }\label{fig:epsilonplots1}
\end{figure}
We see that in the left panel of Fig.~\ref{fig:epsilonplots1} 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.~\ref{fig:epsilonplots1} further shows that the scalings \eqref{Eq:scaling1att1} for $\varepsilon$ and $\tilde{\varepsilon}$ (see below) come out to be as predicted in the limit $\tau\rightarrow0$.

\section{Mapping the attractor surface at early and late times}\label{sec:early}

\subsection{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 \eqref{Eq:s1s2def}.

Crucially, we need
\begin{equation}\label{Eq:condsn}
    0 <\sigma_1 < 1 \quad {\rm and} \quad \frac{3}{4}\sigma_1 + \frac{1}{4} < \sigma_2 < 1.
\end{equation}
The relations $\sigma_{1,2} <1 $ are imposed by causality \cite{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.~\ref{fig:epsilonplots1}. 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 \eqref{Eq:condsn}, 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:
 \begin{align}\label{Eq:scaling1att1}
\gamma \varepsilon &\rightarrow \sqrt{\frac{r-1}{r}}, \quad \quad
\gamma\tilde{\varepsilon} \rightarrow{{g_2}} \tau ^{\frac{4}{3} (\sigma_2 - \sigma_1)},\nonumber\\
\chi &\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}},
\end{align}
and the following scaling behaviors of the components of the effective metrics:
\begin{align}\label{Eq:scaling1att2}
a  &\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)}, \nonumber\\
\tilde{a}  &\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)}.
\end{align}
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:
\begin{align}
    & 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}},\nonumber\\
    & a_{10} = \sqrt{\frac{r-1}{r}} b_{10},\nonumber\\
    & c_{10} = b_{10},\nonumber\\
    & 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},\nonumber\\
    & b_{20} =   \sqrt{\frac{4r-1 + 2(r-1)\sigma_1}{3 {r(r-1)}}} b_{10},\nonumber\\
    & c_{20} =  \sqrt{\frac{4r-1 - 4(r-1)\sigma_1}{3 {r(r-1)}}}b_{10}, \nonumber\\
    & 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}.
\end{align}
As a result, in the limit $\tau\rightarrow0$, the following five identities should hold
\begin{align}\label{grand-identity}
&\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}}}\\\nonumber &=\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)}}
\end{align}
on any curve on the attractor surface. Fig.~\ref{fig:abcallplots1} demonstrates the numerical check of the above identities on an attractor curve.

The stricter lower bound on $\sigma_2$ in Eq. \eqref{Eq:condsn} 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 \eqref{Eq:scaling1att2}, it is clear that the lower bound is then self-consistent. Furthermore, it is also clear from \eqref{Eq:scaling1att1} 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. \eqref{Eq:condsn} we impose
\begin{equation}\label{Eq:condsn2}
    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,
\end{equation}
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 \eqref{Eq:condsn} can be readily satisfied with reasonable weak coupling and strong coupling parameters. Therefore, we do not further investigate the above condition.

\begin{figure}[ht]\centering{
\includegraphics[width=0.8\linewidth]{abcplotsToni.pdf}}
\caption{
Check of the identities \eqref{grand-identity} at early times. Similar identities given in \eqref{Eq:grand-identity2} also holds for the disperser surface, and can be verified to a very high accuracy numerically as evident from Fig. \ref{fig:abcplotsdisp}.}\label{fig:abcallplots1}
\end{figure}

It is easy to see from the above that the physical energy densities of the two subsystems behave as follows at early time
\begin{equation}\label{eq:E1}
    \mathcal{E}_1 = \varepsilon \frac{a b^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}},
\end{equation}
and therefore
\begin{equation}
  \frac{  \mathcal{E}_2}{\mathcal{E}_1}\sim \tau^{\frac{8}{3}(\sigma_2- \sigma_1)}.
\end{equation}
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 \textit{bottom-up thermalization} scenario \cite{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 \eqref{Eq:condsn}, depending on the choice of parameters.

While all initial conditions evolve towards specific curves on the attractor surface, they also get repelled from a \textit{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 \eqref{Eq:scaling1att1}, $\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~\ref{sec:double-disp}. One can also find two other such \textit{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.

\subsection{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 \cite{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. \ref{fig:upper-bound-plot}, 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. \ref{fig:upper-bound-plot}, 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.

\begin{figure}
    \centering
    \includegraphics[width=0.48\textwidth]{hard_sector_upper_bound_on_attractor}
    \includegraphics[width=0.48\textwidth]{evsE.pdf}
    \caption{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 \eqref{parameters} and \eqref{parameters2}, 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$.
    }
    \label{fig:upper-bound-plot}
\end{figure}

\subsection{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
\begin{equation}
    \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}.
\end{equation}
(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
\begin{equation}
\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},
\end{equation}
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:
\begin{equation}
\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}.
\end{equation}

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
\begin{align}\label{alpha}
\alpha := & \gamma(\epsilon_{20}\tau_0^{4/3}) \sqrt{\epsilon_{10} \tau_0^{4/3}} = \gamma \epsilon_{20}\sqrt{\epsilon_{10}} \tau_0^{2}, \\
\beta := & \gamma(\epsilon_{10}\tau_0^{4/3}) \sqrt{\epsilon_{20} \tau_0^{4/3}} = \gamma \epsilon_{10}\sqrt{\epsilon_{20}} \tau_0^{2}. \label{beta}
\end{align}
These two parameters uniquely label each attractor (hydrodynamic) curve as shown below. It is also useful to define
\begin{equation}
\mathfrak{r}:= (\epsilon_{20}/\epsilon_{10})^{1/4} =\sqrt{\alpha/\beta} = \lim_{\tau\rightarrow\infty}{(\tilde{\varepsilon}/\varepsilon)^{1/4}}.
\end{equation}
For the sake of convenience, we define dimensionless time:
\begin{equation}
\omega_1 := \epsilon_{1,pf}^{1/4}\tau, \quad \omega_2 := \epsilon_{2,pf}^{1/4}\tau.
\end{equation}
Note however that $\omega_1$ and $\omega_2$ are not independent time variables. In fact
\begin{equation}
\omega_2 = \mathfrak{r} \omega_1.
\end{equation}
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
\begin{align}
\varepsilon^{1/4}\tau  &= \omega_1 - \frac{2}{3} C_{\eta}  + \frac{3 \alpha-4 C_{\eta} C_{\tau} }{18} \frac{1}{\omega_1} + \mathcal{O}(\omega_1^{-2}),  \label{invert}\\
\chi &= \frac{4C_{\eta} }{3\omega_1} +  \frac{8C_{\eta} C_{\tau} }{9\omega_1^2} +  \mathcal{O}(\omega_1^{-3}) ,\\
\tilde{\varepsilon}^{1/4}\tau  &= \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}), \\
\tilde{\chi} &= \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}),
\end{align}
The above equations show all terms up to second order in derivatives. Defining
\begin{equation}
w_1:= \varepsilon^{1/4}\tau, \quad w_2:= \tilde{\varepsilon}^{1/4}\tau
\end{equation}
we can invert \eqref{invert}
\begin{equation}
\omega_1 = w_1 + \frac{2}{3} C_{\eta}  - \frac{3 \alpha-4 C_{\eta} C_{\tau} }{18} \frac{1}{w_1} + \mathcal{O}(\omega_1^{-2}).
\end{equation}
In terms of $w_1$ we can readily express the hydrodynamic expansion as a phase space curve
\begin{align}
\chi &= \frac{4C_{\eta} }{3w_1} +  \frac{8C_{\eta} C_{\tau} }{9w_1^2} +  \mathcal{O}(w_1^{-3}) ,\\
w_2
&=\sqrt{ \frac{\alpha}{\beta}}w_1 - \frac{2}{3} \left(\tilde{C}_{\eta} - \sqrt{ \frac{\alpha}{\beta}} C_{\eta} \right)\\&+\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} -24 C_{\eta} ^2}{18} \right)\frac{1}{ w_1} + \mathcal{O}(w_1^{-2}), \nonumber \\
\tilde{\chi}
&=\sqrt{ \frac{\beta}{\alpha}} \frac{4\tilde{C}_{\eta} }{3 w_1} +  \frac{\beta}{\alpha}\frac{8\tilde{C}_{\eta} \left(\tilde{C}_{\eta} +C_{\eta} \sqrt{ \frac{\beta}{\alpha}}+\tilde{C}_{\tau} \right)}{9 w_1^2} +  \mathcal{O}(w_1^{-3}).
\end{align}
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:
 \begin{align}
  a &= 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}), \\
 b &= 1 + \frac{1}{6}\alpha\frac{1}{\omega_1^2} + \mathcal{O}(\omega_1^{-4}),\nonumber\\
 c &= 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}),\\
 \tilde{a} &= 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}), \\
 \tilde{b} &= 1 + \frac{1}{6}\frac{\beta^2}{\alpha}\frac{1}{\omega_1^2} +\mathcal{O}(\omega_1^{-4}),\nonumber\\
\tilde{c} &= 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}).
 \end{align}
Remarkably, up to $\mathcal{O}(w_1^{-4})$ the volume factors are:
\begin{equation}
a b^2 c = 1 + \mathcal{O}(w_1^{-4}), \quad \tilde{a} \tilde{b}^2 \tilde{c}= 1 + \mathcal{O}(w_1^{-4}).
\end{equation}

We can readily compute the full system energy density which behaves as
\begin{equation}\label{Eq:wtotal}
\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}),
\end{equation}
where
\begin{equation}
\omega := \mathcal{E}_{pf}^{1/4}\tau = (\epsilon_{1,pf}+\epsilon_{2,pf} )^{1/4} \tau.
\end{equation}
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 \textit{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)
\begin{equation}\label{Bjorken-gen-eqn}
\mathcal{E}' = -\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),
\end{equation}
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
\begin{equation}\label{Bjorken-trT}
T^\mu_{\,\mu} = 3 \mathcal{P}(\mathcal{E}) -\mathcal{E}  - \frac{3}{\tau} \zeta(\mathcal{E}).
\end{equation}
We can assume that
\begin{align}\label{fullP-form}
\mathcal{P}(\mathcal{E}) &= \frac{\mathcal{E}}{3} + k \mathcal{E}^2 + \cdots\\
\eta(\mathcal{E}) &= \kappa \mathcal{E}^{3/4} + \cdots, \quad \zeta(\mathcal{E}) = \tilde\kappa \mathcal{E}^{3/4} + \cdots.
\end{align}
with $k$,  $ \kappa $ and $\tilde\kappa$ constants. Solving \eqref{Bjorken-gen-eqn} with inputs from \eqref{fullP-form} in the large proper time expansion, we obtain
\begin{align}\label{Eq:EtotalSingle}
\mathcal{E} =& \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}
\nonumber\\ & +\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,
\end{align}
and also
\eqref{Bjorken-trT}
\begin{equation}\label{Tr-T-1}
{\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.
\end{equation}

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. \eqref{Eq:EtotalSingle} with Eq. \eqref{Eq:wtotal}, we explicitly find that the effective viscosity of the full system is
\begin{equation}\label{eff-shear}
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}}
\end{equation}
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 \eqref{eff-shear}:
\begin{equation}\label{kappa}
\kappa = \frac{4}{3}C_{\eta, eff}
\end{equation}
Similarly, the hydrodynamic expansion of our hybrid system yields that
\begin{equation}\label{Tr-T-2}
{\rm Tr}(T) = - \frac{8\tau_0^{8/3}}{3\tau^{8/3}}\gamma \epsilon_{10}\epsilon_{20}+ \cdots.
\end{equation}
We can immediately conclude that
\begin{itemize}
    \item $\mathcal{E}_0 = \epsilon_{10}+ \epsilon_{20}.$
\item $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},$
\item and $\tilde\kappa = 0,$
\end{itemize}
i.e. the effective bulk viscosity vanishes since the $\tau^{-2}$ term vanishes in \eqref{Tr-T-2}.

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. \ref{fig:epsilonplotsfull2} for the first-order hydrodynamic tails of $\varepsilon$ and $\tilde{\varepsilon}$ for the first and second attractors, respectively.

\begin{figure}
\centering
\begin{minipage}{0.48\textwidth}
\includegraphics[width=\linewidth]{epsilonplotsfull11.pdf}
\end{minipage}
\caption{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)$ }\label{fig:epsilonplotsfull2}
\end{figure}
At late time
\begin{equation}
\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}).
\end{equation}
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
\begin{equation}
\epsilon_{10}\tau_0^{4/3} \approx 1.346, \, \,\epsilon_{20}\tau_0^{4/3} \approx 0.370, \,\, \alpha = 0.424,\,\,\beta = 0.801.
\end{equation}
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.~\ref{fig:chiplotsfull2} for the first-order hydrodynamic tails of $\chi$ and $\tilde{\chi}$ for the first and second attractors, respectively.

\begin{figure}[ht]
\centering
\begin{minipage}{0.48 \textwidth}
\includegraphics[width=\linewidth]{chiplotsfull11.pdf}
\end{minipage}
\caption{The hydrodynamic tail showing the onset of first-order hydrodynamics with the same values as in Fig.~\ref{fig:epsilonplotsfull2}.}\label{fig:chiplotsfull2}
\end{figure}

We can also match using leading order hydrodynamic expansion of {$\chi$ ($\tilde{\chi}$)}:
\begin{align}
\ln {\chi} \approx - 2 \ln \tau + \ln\frac{16C_{\eta} }{9} +3/4\ln (\epsilon_{10}\tau_0^{4/3}),\nonumber\\ \ln {\tilde{\chi}} \approx - 2 \ln \tau + \ln \frac{16\tilde{C}_{\eta} }{9}  +3/4\ln (\epsilon_{20}\tau_0^{4/3}).
\end{align}
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.

\subsection{Late time approach to attractor: fluctuation analysis}\label{sec:late}\label{sec:fluc}

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:
\begin{align}
\varepsilon(\hat{\tau}) &={\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}),\\
\phi(\hat{\tau}) &= 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}) ,\\
a(\hat{\tau}) &= 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}) ,\\
b(\hat{\tau}) &= 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}),\\
c(\hat{\tau}) &= \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}),
\end{align}
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~\ref{app:grad-exp}, 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~\ref{app:grad-exp}. 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:
\begin{align}\label{sol-phi1}
\delta \phi(\hat{\tau}) &=m_1 \; e^{-\frac{3 {\epsilon_{10}}^{1/4} \; \tau_0}{2 C_{\tau}}  \hat{\tau}^{2/3}}  \hat{\tau}^{-\frac{4}{3} +
\frac{2 C_{\eta}}{3 C_{\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}}},\\
\label{sol-phi2}
\delta \tilde{\phi}(\hat{\tau}) &= n_1 \; e^{-\frac{3 {\epsilon_{10}}^{1/4} \; \tau_0}{2 C_{\tau}}  \hat{\tau}^{2/3}} \hat{\tau}^{
-\frac{10}{3} +
\frac{2 C_{\eta}}{3 C_{\tau}}}   + n_2 \; e^{-\frac{3{\epsilon_{20}}^{1/4} \; \tau_0}{2 C_{\tau_2}}  \hat{\tau}^{2/3}}  \hat{\tau}^{-\frac{4}{3} +
\frac{2 \tilde{C}_{\eta}}{3 \tilde{C}_{\tau}}}
\end{align}
where $m_1, m_2, n_1, n_2$ are constants, of which only two are independent:
\begin{eqnarray}
\label{cp1-c1}
m_1 &=& n_1 \frac{3 C_{\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 ,\nonumber\\
\label{cp2-c2}
n_2 &=& m_2 \frac{3 \tilde{C}_{\eta} \left( - \tilde{C}_{\tau} \; {\epsilon_{10}}^{1/4} + C_{\tau} \; {\epsilon_{20}}^{1/4} \right)}{2 C_{\eta} \tilde{C}_{\tau}^2  {\epsilon_{10}} \gamma}\tau_0.
\end{eqnarray}

If $C_{\tau} \: {\epsilon_{20}}^{1/4} > \tilde{C}_{\tau} \: {\epsilon_{10}}^{1/4}$, i.e. if
\begin{equation}\label{Eq:poss1}
    \mathfrak{r}=\sqrt{\alpha/\beta} > \tilde{C}_{\tau}/C_{\tau}
\end{equation}
then
\begin{equation}
\exp{\left(-\frac{3 {\epsilon_{10}}^{1/4} \; \tau_0}{2 C_{\tau}}  \hat{\tau}^{2/3}\right)} =\exp{\left(-\frac{3 {\epsilon_{10}}^{1/4} \; \tau_0^{1/3}}{2 C_{\tau}}  \tau^{2/3}\right)} = \exp{\left(-\frac{3}{2 C_{\tau}}\left(\frac{\beta}{\sqrt{\alpha\gamma}}\right)^{1/3}\tau^{2/3}\right)}
\end{equation}
will be the dominant exponential and give the relaxation towards the attractor. The powers of $\hat{\tau}$ accompanying this exponential in \eqref{sol-phi1} and \eqref{sol-phi2} 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
\begin{equation}\label{Eq:poss2}
   \mathfrak{r}=\sqrt{\alpha/\beta}< \tilde{C}_{\tau}/C_{\tau},
\end{equation}
then
\begin{equation}
\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)}
\end{equation}
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 \eqref{Eq:poss1} is realized, the strongly coupled sector hydrodynamizes faster and the reverse occurs when \eqref{Eq:poss2} holds. However, realizing \eqref{Eq:poss2} 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.
\begin{equation}
    \mathfrak{r}=\sqrt{\alpha/\beta} = \tilde{C}_{\tau}/C_{\tau} ,
\end{equation}
the two independent solutions to leading order are:
\begin{align}\label{sol-phi1-case2}
\delta \phi(\hat{\tau}) &= e^{-\frac{3 {\epsilon_{10}}^{1/4} \tau_0}{2 C_{\tau}}  \hat{\tau}^{2/3}} \left(   m_1  \; \hat{\tau}^{-\frac{4}{3} + \frac{2 C_{\eta}}{3 C_{\tau}}}   +  m_2 \; \hat{\tau}^{-\frac{8}{3}  + \frac{2 \tilde{C}_{\eta}}{3 \tilde{C}_{\tau}} }  \right),\\
\label{sol-phi2-case2}
\delta \tilde{\phi}(\hat{\tau}) &= e^{-\frac{3 {\epsilon_{10}}^{1/4} \tau_0}{2 C_{\tau}}  \hat{\tau}^{2/3}} \left(   n_1  \;  \hat{\tau}^{-\frac{8}{3} + \frac{2 C_{\eta}}{3 C_{\tau}}}   + n_2 \;  \hat{\tau}^{- \frac{4}{3}  + \frac{2 \tilde{C}_{\eta}}{3 \tilde{C}_{\tau}} }  \right),
\end{align}
where $m_1 $, $ m_2$, $n_1$, and $n_2$ are arbitrary constants only two of which are independent:
\begin{align}
m_1 = - \frac{C_{\eta} C_{\tau}^2 \left( - \tilde{C}_{\eta} C_{\tau} + C_{\eta} \tilde{C}_{\tau} - 2 C_{\tau} \tilde{C}_{\tau}  \right)}{\tilde{C}_{\eta} \tilde{C}_{\tau}^4  {\epsilon_{10}}  \gamma} n_1 ,\nonumber\\
n_2 = - \frac{\tilde{C}_{\eta}  \left( \tilde{C}_{\eta} C_{\tau} - C_{\eta} \tilde{C}_{\tau}  - 2 C_{\tau} \tilde{C}_{\tau}  \right)}{C_{\eta}  \tilde{C}_{\tau} ^4  {\epsilon_{10}}  \gamma}m_2.
\end{align}

\subsection{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 \eqref{alpha} and \eqref{beta}, 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
\begin{align}
    \alpha&\sim (\gamma \mathcal{E})^{3/2},\\
    \beta&\sim (\gamma \mathcal{E})^{3/2},
    \end{align}
where $ \mathcal{E}$ is the total initial energy density. Using the definitions of the total entropy density and the individual subsystem entropy densities as
\begin{align}
    S_{tot} &=s  +\tilde{s},\\
    s&=\frac{4}{3 \tau} \varepsilon^{3/4} b^2 c,\\
    \tilde{s}&=\frac{4}{3 \tau}  \tilde{\varepsilon}^{3/4} \tilde{b}^2 \tilde{c},
\end{align}
we can readily see from the hydrodynamic expansions discussed above that the individual system entropy densities behave at late time as
\begin{align}
    \lim_{\tau\rightarrow\infty}\tau s = \beta/\sqrt{\alpha}, \quad \lim_{\tau\rightarrow\infty}\tau \tilde{s} = \alpha/\sqrt{\beta},
\end{align}
which means that the total entropy scales like
\begin{align} \lim_{\tau\rightarrow\infty}    \tau S_{\rm tot} \equiv (   \tau S_{\rm tot} )_{\rm hydro}&\sim (\gamma \mathcal{E})^{3/4}.
\end{align}
We can compare this to the pocket formula in \cite{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 \cite{Giacalone:2019ldn} states that the total particle multiplicity per unit rapidity\footnote{Here, the rapidity, $\eta$, should not to be confused with the shear viscosity.} follows
\begin{align}
    \frac{dN}{d\eta} \approx  (S\tau)_{\rm hydro}.
\end{align}
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) \cite{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.~\ref{fig:entropy}, 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.~\ref{fig:entropy},  with the soft (strongly self-interacting) sector dominating the entropy production at later time.

\begin{figure}[tbp]
\begin{minipage}{0.5\textwidth}
\includegraphics[width=\textwidth]{entropy_inset.pdf}\end{minipage}
\begin{minipage}{0.5\textwidth}
\includegraphics[width=\textwidth]{ratio_entropy.pdf}
\end{minipage}
\caption{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} = 4 C_{\tau}, \phi (\tau_0) =0.0220931202, \tilde{\phi} (\tau_0) =0.0297056408397,
    \varepsilon (\tau_0) = 0.05, \tilde{\varepsilon} (\tau_0) = 0.1$.
}\label{fig:entropy}
\end{figure}

{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{\vert{\rm Tr} T\vert}{\mathcal{E}}=\frac{\mathcal{E}-P_L - 2P_T}{\mathcal{E}}\sim \tau^{-4/3},$$ as follows from \eqref{Eq:wtotal} and \eqref{Tr-T-2}.  An illustrative example for the full time-evolution of the interaction energy is shown in Fig.~\ref{fig:trace}. 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.}

\begin{figure}
\center
\includegraphics[width=0.75\linewidth]{trace.pdf}
\caption{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$.}\label{fig:trace}
\end{figure}

\subsection{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 \eqref{alpha} and \eqref{beta}) 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},$
\begin{eqnarray} \label{chis}
\chi = \frac{4 C_{\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)\, , \nonumber  \\
\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)\, .
\end{eqnarray}
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
\begin{equation*}
    \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)
\end{equation*}
indicating the series is asymptotic in nature with zero radius of convergence, which is seen in subsystems in the top panels of Fig.~\ref{fig:slope_pade}. 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 \cite{Aniceto:2011nu}.
The Borel transform of the series is given by
\begin{eqnarray} \label{gborel}
 \chi_{B,i} (\xi)= \sum_{n=0}^{\infty} \frac{r_{i,n}}{n!} {\xi}^n
 \end{eqnarray}
 $ r_{i,n}$ are the coefficients from the series (\ref{chis}). 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{\'e} approximant
 \begin{equation}
 \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}
 \end{equation}
where the coefficients of the Pad{\'e} approximant are fixed to agree with series coefficients (\ref{gborel}). The Pad{\'e} 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\'e-approximant in the Borel plane and {they depend} on the dimensionless parameter $ \alpha$ and $ \beta$.  In Fig.~\ref{fig:slope_pade}, 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 \cite{Heller:2015dha}.   The starting value of the branch cut is proportional to the inverse of the slope of the coefficients shown in Fig.~\ref{fig:slope_pade} and can be identified with the non-hydrodynamic mode that decays exponentially at late-times. From the analysis in Sec.~\ref{sec:fluc} we readily see that $ \xi_{1,0} = \xi_{2,0} =\frac{ 3 } {2 C_{\tau}} \Big(\frac{\beta}{ \sqrt{\alpha \gamma}} \Big)^{1/3}= 3\epsilon_{10}^{1/4}\tau_0^{1/3} /(2 C_{\tau})=\xi_0 $ when the condition \eqref{Eq:poss1} is satisfied, and $ \xi_{1,0} = \xi_{2,0} =\frac{ 3 } {2 C_{\tau}} \Big(\frac{\alpha}{ \sqrt{\beta \gamma}} \Big)^{1/3}= 3\epsilon_{20}^{1/4}\tau_0^{1/3} /(2 C_{\tau})=\xi_0 $ when the condition \eqref{Eq:poss2} holds.

With the analytic continuation of the Borel transformed function, we can now perform the inverse integration \cite{Aniceto:2011nu,Heller:2015dha,Casalderrey-Solana:2017zyh}
\begin{eqnarray} \label{Eq:inverse_Borel}
\chi_i (x) = \frac{1}{x}  \int_\mathcal{C} e^{-\xi/x} P_{N,M,i} ({\xi}) d \xi
\end{eqnarray}
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. \ref{fig:slope_pade}. We have two Stokes coefficients corresponding to the two exponential decay modes discussed in Sec. \ref{sec:fluc} 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.

\begin{figure}[ht]
\includegraphics[width=0.5\linewidth]{Strong_slope1.pdf}
\includegraphics[width=0.5\linewidth]{Weak_slope1.pdf} \\
\includegraphics[width=0.5\linewidth]{Strong_pade1.pdf}
\includegraphics[width=0.5\linewidth]{Weak_pade1.pdf}
\caption{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 } {2 C_{\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$.}\label{fig:slope_pade}
\end{figure}

\section{Hydrodynamization on the attractor surface}\label{sec:hydrodynamization}
Figures \ref{fig:epsilonplotsfull2} and \ref{fig:chiplotsfull2} 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.

\subsection{Hydrodynamization time}\label{sec:hydro-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 \cite{Attems:2017zam}. For the hybrid hydrodynamic attractor, we define the hydrodynamization time of each of the sectors through
\begin{eqnarray}
\frac{|\Delta P_L|}{P} := \frac{|\phi-\phi_{1st}|}{P}< 0.1, \hspace{0.4cm}\text{for} \hspace{0.4cm} \tau>\tau_{hd},
\end{eqnarray}
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 \ref{fig:ratio_5:1}, where the blue curve corresponds to the hard sector and the red one to the holographic sector.

\begin{figure}[H]
\centering
\includegraphics[scale=0.8]{ratio5_1tau.pdf}
\caption{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.} \label{fig:ratio_5:1}
\end{figure}

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 \ref{fig:ratio_5:1}).
{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 \textit{universal} as shown in Fig \ref{fig:hardconformality} and the same as in conformal attractors studied earlier in the literature. However for the holographic sector as in Fig \ref{fig:softconformality}, 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.

\begin{figure}[ht]
\centering
\includegraphics[scale=0.9]{hardhydro_universality.pdf}
\caption{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.} }\label{fig:hardconformality}
\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[scale=0.9]{softhydro_universality.pdf}
\caption{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 \textbf{\textcolor{Green}{0.64}}, \textbf{\textcolor{BrickRed}{0.76}}, \textbf{\textcolor{Blue}{0.93}}, \textbf{ \textcolor{RedViolet}{1.07}}, \textbf{\textcolor{SpringGreen}{1.21}}, \textbf{ \textcolor{cyan}{1.76}} to \textbf{\textcolor{Sepia}{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.
}\label{fig:softconformality}
\end{figure}

{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 \textit{universality} as in Fig \ref{fig:ratiol1softconformality} and coinciding with the usual conformal attractors. Meanwhile from Fig \ref{fig:ratiol1hardconformality}, 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.

\begin{figure}[ht]
\centering
\includegraphics[scale=0.9]{hd-soft-conformal.pdf}
\caption{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.~\ref{fig:hardconformality} corresponding to the cases where the hard sector was dominant at initial reference time.} }\label{fig:ratiol1softconformality}
\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[scale=0.9]{hd-hard-conformal.pdf}
\caption{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.~\ref{fig:softconformality}, 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.} }\label{fig:ratiol1hardconformality}
\end{figure}

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 \ref{fig:ratio_5:1}

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 \ref{fig:ratio_5:1}.

\subsection{Rough phenomenological match to $p$-$p$ and Pb-Pb collisions }\label{sec:pheno-pp-pb}

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 \cite{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
\begin{equation}
    Q_{s,A}^2\sim 0.4 (x/0.001)^{-0.3} A^{1/3} Q_{s,p}^2.
\end{equation}
For Pb-Pb collisions, where $A^{1/3}\approx 6$, we take simply
\begin{equation}
Q_{s,Pb}^2 = 2.5\,  Q_{s,p}^2
\end{equation}
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 \ref{pp-table} and \ref{pb-table} show the hydrodynamization time of {accordingly defined} $p$-$p$ and Pb-Pb collisions for various ratios of initial energy densities of the two subsectors.

\begin{table}[H]
\centering
 \begin{tabular}{|c| c| c| c| c|c|c|c|c| }
 \hline
$\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})$ \\ [0.5ex]
 \hline\hline
 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\\[1ex]
 \hline
 \end{tabular}
 \caption{Hydrodynamization times for $p$-$p$ initial conditions.
{Here $w= \mathcal{E}_1^{1/4}\tau$, where $\mathcal{E}_1$ in \eqref{eq:E1}. and $w_{hd}=w(\tau_{hd})$ and similarly for the second subsystem.}}\label{pp-table}
\end{table}

\begin{table}[H]
\centering
 \begin{tabular}{|c| c| c| c| c| c|c|c|c|}
 \hline
$\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})$\\ [0.5ex]
 \hline\hline
 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 \\[1ex]

 \hline
 \end{tabular}
\caption{Hydrodynamization times for Pb-Pb initial conditions.}\label{pb-table}
\end{table}

\begin{figure}[H]
\centering
\includegraphics[scale=0.72]{hardhydropheno.pdf}
\includegraphics[scale=0.72]{softhydropheno.pdf}
\caption{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.} \label{fig:phenopp-pb}
\end{figure}

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 \ref{fig:phenopp-pb}, 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 \ref{fig:ratio_5:1} as:
\begin{itemize}
    \item 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.

    \item 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})$.
\end{itemize}
For details on the evolution of energy densities in these scenarios, see Appendix \ref{sec:thermalization}.

\section{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 \cite{Naik:2021yph,Coquet:2021lca} and hadronization \cite{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,\footnote{See \cite{Grossi:2020ezz,Grossi:2021gqi,Florio:2021jlx} for studies of effects of stochastic critical fluctuations on hydrodynamics, \cite{Kovtun:2012rj} for a review on stochastic statistical fluctuations in hydrodynamics and \cite{Dore:2020jye} for an investigation of potential hydrodynamic attractors near the critical point. See also \cite{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}.

\acknowledgments

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.