Hydrodynamization in hybrid Bjorken flow attractors

Hybrid fluid models, consisting of two sectors with more weakly and more strongly self-interacting degrees of freedom coupled consistently as in the semi-holographic framework, have been shown to exhibit an attractor surface for Bjorken flow. Retaining only the simple viscid fluid descriptions of both sectors, we find that, on the attractor surface, the hydrodynamization times of both subsectors decrease with increasing total energy density at the respective point of hydrodynamization following a conformal scaling, reach their minimum values, and subsequently rise rapidly. The minimum values are obtained when the respective energy densities are of the order of the inverse of the dimensionful inter-system coupling. Restricting to attractor curves which can be matched to glasma models at a time set by the saturation scale for both $p$-$p$ and Pb-Pb collisions, we find that the more weakly coupled sector hydrodynamizes much later, and the strongly coupled sector hydrodynamizes earlier in $p$-$p$ collisions, since the total energy densities at the respective hydrodynamization times of these sectors fall inside and outside of the conformal window. This holds true also for phenomenologically relevant solutions that are significantly away from the attractor surface at the time we match to glasma models.

E Different initial ratios of subsector energy densities 44 1 Introduction One of the most remarkable discoveries in non-equilibrium dynamics recently has been that quantum many-body systems can hydrodynamize far away from equilibrium [1][2][3][4][5][6][7][8].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.In fact, hydrodynamization can happen in both weakly coupled and strongly coupled systems even when they are far away from equilibrium and have large pressure anisotropies [9].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 [3,7,[10][11][12][13][14][15][16][17][18][19][20][21], see [22] 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.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.[23,24] 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 [25] and subsequently kinetic theory [26,27] is applicable, to later epochs when the soft gluons form a strongly coupled thermal bath.The semi-holographic framework, briefly reviewed in Sec. 2, allows for simple phenomenological constructions where the evolution of both perturbative (weakly self-interacting) and non-perturbative (strongly self-interacting) degrees of freedom can be described simultaneously in a consistent way [28][29][30][31][32][33][34].In the simplest model, one retains only two effective metric couplings, and reduces both sectors to viscous fluids each following the Müller-Israel-Stewart (MIS) [35,36] 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 γ, γ = rγ) 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. 1he existence of hydrodynamic attractors in such hybrid fluid models has been established in [37].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 feature of 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 given the results obtained in [34].
the attractor surface is an initial evolution that is reminiscent of the so-called bottom-up thermalization scenario of Ref. [38], for at sufficiently early time, the energy density of the perturbative sector always dominates over the strongly self-interacting non-perturbative sector.This "bottom-up-like" feature is a consequence of the effective fluid description of the semi-holographic setup, which however can be expected to be a good approximation to a full semi-holographic treatment only at τ ≳ γ 1/4 .Nevertheless, it is a robust feature of the resulting hybrid attractor itself, where the dominance of the weakly coupled sector at earlier times is determined by an exponential factor with an exponent involving only the viscosities and relaxation times of the two sectors.Apart from the overall energy scale, the ratio of the respective energy densities at a conveniently chosen early reference time τ 0 turns out to be a useful parameter for the curves which span the attractor surface, even though τ 0 is outside the phenomenologically relevant range.
In this paper, we systematically study hydrodynamization in these hybrid hydrodynamic attractors following preliminary observations in [37] 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 [39,40]) by a detailed exploration of the attractor surface.In the application to such high-energy collisions, we consider only times τ ≳ γ 1/4 , since a fluid description cannot be expected to be applicable at arbitrarily small times.In particular, our systematic study of hydrodynamization on the attractor surface suggests a novel explanation for rapid hydrodynamization of the softer degrees of freedom in small system collisions as summarized below.
The general characterization of the hydrodynamization times of the two components on the hybrid attractor surface in our simplified model is as follows.Let us denote the hydrodynamization times of the perturbative and non-perturbative (holographic) sectors as τ hd and τ hd , respectively.When the energy density at the hydrodynamization times is in the regime E(τ hd ), E( τ hd )≪γ −1 , we find that the hydrodynamization times of both sectors follow conformal behaviors, i.e.
with p and p being constants.Generically, due to the bottom-up-like feature discussed above, the energy density of the perturbative sector is dominant at some reference time τ 0 ≪ γ 1/4 .For such cases, we find that p ≈ 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, 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 reference time (say) τ 0 = 0.001γ 1/4 is increased from 1 : 1 to 1000 : 1, the constant 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 reference time, the roles in the conformal window are reversed.The p of the holographic sector is then universally ≈ 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 γ −1 , i.e.
For larger E(τ hd ), E( τ 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 matching time τ i ∼ γ 1/4 as corresponding to small systems (high-multiplicity p-p or A-p collisions) [39,40] as opposed to A-A heavy-ion collisions.To quantify our results, we identify the energy scale of the intersystem coupling γ −1/4 with the saturation scale Q p-p s of p-p collisions, and set various initial conditions such that the total energy density ∝ Q 4 s at τ ≈ Q −1 s , with Q s in A-A systems chosen from typical values in the IPsat and bCGC models [41].Generically, we find: • The total energy densities at hydrodynamization time of the perturbative sector (E(τ hd )) are within the conformal window approximately for both p-p and Pb-Pb collisions.Since E(τ hd ) is smaller in p-p collisions, naturally the hydrodynamization time of the perturbative sector is significantly larger in p-p collisions as readily follows from the scaling behavior of τ hd in the conformal window.
• The total energy densities at hydrodynamization time of the holographic sector (E( τ hd )) are outside the conformal window for both p-p and Pb-Pb collisions, i.e.E( τ hd ) > γ −1 .Since E(τ hd ) is smaller in p-p collisions, the hydrodynamization time of the softer sector is smaller in p-p collisions than in the case of Pb-Pb collisions.This follows from the growth of τ hd with E( τ hd ).
Furthermore, we have established that these results remain valid for generic phenomenologically relevant initial conditions which can be matched to the glasma models at time scales of the order the saturation scale for both p-p and Pb-Pb type collisions.For a large class of such initial conditions, our model exhibits an attractor of a system with both weakly and strongly coupled components that is driving the systems towards hydrodynamic behavior well before the onset of hydrodynamics.Our simple hybrid model therefore provides a new perspective on how collective flow can emerge from the softer components even in small system collisions.The organization of the paper is as follows.In Sec. 2, we review our set-up of two fluids undergoing Bjorken flow and coupled mutually via the democratic effective metric coupling [32,33].We further discuss the method to discern the attractor surface in this hybrid system.In Sec. 3, we provide more analytic details of early time behavior on the attractor surface, before turning our attention to the late-time hydrodynamic behavior.In both cases, we expand the discussion beyond our earlier results in [37].Finally in Sec. 4, we systematically study the hydrodynamization on the attractor surface and connect to the phenomenology of p-p and Pb-Pb collisions.

The hybrid action
A minimalist set-up required to obtain our hybrid fluid model involves only two leading democratic effective metric couplings [32,33].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 [33] Above S and S denote the Wilsonian effective actions for the perturbative and nonperturbative sectors respectively, with respective degrees of freedom ψ and ψ, and living in respective effective metrics g µν and gµν , while the physical metric of the full theory is G µν (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 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 g.In what follows, we will not need any explicit description of both S and S. Since both S and S are individually diffeomorphism-invariant in the respective background metrics g and g, we obtain the two Ward identities where the energy momentum tensors are and ∇ and ∇ are built out of g and g, respectively.
The interaction term S int in (2.1) relates these two different background metrics in terms of the actual physical background metric G ultralocally, as follows.We readily see that the two effective metrics appear in the full action as auxiliary fields, and extremizing the action with respect to them leads to imposing The above yield the following coupling equations where γ and γ ′ are the two dimensionful couplings with mass dimension −4.
, where Λ 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, Λ int can be identified with a saturation scale Q s determining the initial conditions.It is evident from (2.5) that higher order effective metric couplings will be suppressed by powers of (T /Q s ) 4 , where T is an effective temperature determining the scales of the stress tensors.
We can compute the energy momentum tensor of the full system by varying the action (2.1) with respect to the physical background metric G µν , which yields In the first equation above, all indices for subsystem variables are lowered using the respective background metrics, while that of the full system (on the left hand side) is lowered using the physical background metric.It can be explicitly checked that the coupling equations (2.5) and the two Ward identities (2.2) indeed imply the local conservation of the full energy-momentum tensor in the physical background metric, i.e.
The mutual effective metric coupling thus leads to a full energy-momentum tensor (2.6) which can be obtained ultra-locally from the subsystem energy-momentum tensors alone, and does not require any further microscopic input.This ensures that a coarse-grained hydrodynamic description of the full energy-momentum tensor in terms of those of the subsystems.It implements the hypothetical Wilsonian perspective that the coarse-grained perturbative description should determine the complete non-perturbative description at any scale.Therefore, assuming that 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,2 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 [31,34,42].At the first step of iteration, the effective metrics are identified with the background metric G µν , and the two systems are solved independently with given initial data.In the next step, the effective metrics of the two subsystems are updated via the coupling equations (2.5) using the respective subsystem stress tensors obtained from the previous iteration.The subsystems are solved once again with the same initial conditions.The iterations are thus continued with fixed initial conditions until they converge implying that the full energy-momentum tensor is conserved in the physical background metric.
This iterative procedure was first proposed in the context of heavy-ion collisions with glasma-like initial conditions for the perturbative sector and an empty AdS-like initial conditions for the holographic sector [30].Subsequent works [31,34,42] 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 [34,42].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.

Hybrid fluid model
Although the above discussion included an action, it is not necessary to describe the evolution of the system, which is governed by the effective conservation laws (2.2) and the algebraic coupling equations (2.5) in the hydrodynamic limit.The only required inputs are the energy-momentum tensors of each sector, t µν and tµν , determined via the constitutive relations, and the physical background metric G µν .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.
u µ g µν u ν = −1 and ũµ gµν ũν = −1. (2.10) The dissipation tensors, Π µν and Πµν , are taken to be orthogonal to their respective four velocities Π µν u µ = 0 and Πµν ũµ = 0. (2.11) The two fluids are assumed to have dissipative dynamics dictated by Müller-Israel-Stewart (MIS) theory 3 [35,36,43], 4 namely where η and η 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 (2.6) will in general not be conformal, T µν g (B) µν ̸ = 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) [33].

Explicit equations for boost-invariant Bjorken flow
We now provide the explicit form of the complete set of equations, (2.2), (2.5) and (2.12).We will work in the Milne coordinates: τ = √ t 2 − z 2 , x, y, and ζ = arctanh(z/t), where x 3 MIS theory with the Weyl covariant version of the convective derivative, as in Eqs.(2.12), is just a consistent truncation of BRSSS theory where certain non-linear transport coefficients are considered to be vanishing. 4The 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 [44,45] 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.) 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 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 τ .Motivated by the form of the background metric, we make the following boost-invariant ansatz for the effective metrics: (2.14) The four velocities are then , 0, 0, 0 and , 0, 0, 0 . (2.15) Here we consider the dissipation tensors, Π µν and Πµν , to be symmetric, transverse and traceless i.e., g µν Π µν = 0 and gµν Πµν .These conditions imply that there is only one independent component of Π µν for each sector, which we call Π η η = −ϕ (the longitudinal non-equilibrium pressure) following the convention in [12], which means that ) (2.17) Then explicitly, the only non-vanishing components of the conservation equations (2.2) are the MIS equations are (2.12) and finally, the algebraic coupling equations (2.14) are b2 − 1 = γ ab 2 c τ where we introduced r ≡ −γ ′ /γ.We see that in the decoupling limit, namely γ = 0, the right hand side vanishes and the effective metrics reduce to the usual background Milne metric (2.13) in which the two fluids evolve independently.It is clear that the above system actually has four dynamical variables, namely the energy densities and the non-equilibrium pressures (ε, ε, ϕ, φ) which follow the first-order equations (2.18), (2.19), (2.20) and (2.21).The six effective metric variables, a, ã, b, b, c, c, can be eliminated by determining them as functions of the four dynamical variables by solving the six algebraic equations (2.22)-(2.27).The evolution of the full system can thus be readily determined numerically for given initial values of (ε, ε, ϕ, φ).
In what follows, we parametrize the transport coefficients in the two conformal subsystems, and the respective MIS relaxation times via dimensionless quantities C η , Cη , C τ and Cτ defined via In the simulations that follow, we will usually choose 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 unless noted otherwise.In the inviscid Minkowski case [33], 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 ⪅ 1.11, which we discuss in Appendix B. Note that the coupling γ is the only dimensionful scale in our framework.Further, we will work with the dimensionless anisotropy parameter, χ = ϕ/(ϵ + P ).Thus the conservation and the MIS equations read where the prime denotes a derivative with respect to the proper time, e.g.ε ′ = ∂ τ ε.The full energy momentum tensor (2.6) is given by where the total energy density, total transverse pressure and longitudinal pressure are (2.39)

Determining the attractor surface
As discussed above, the hybrid fluid system has a four-dimensional phase space with the energy densities and non-equilibrium pressures, (ε, ε, ϕ, φ), as the dynamical variables.It was shown in [37] that the full system has a two-dimensional attractor surface.For any initial condition (ε, ε, ϕ, φ)(τ 0 ), the system evolves towards a particular curve on the attractor surface.The curves which compose this attractor surface are characterized by .Anisotropies χ 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 [37], where more plots are given for this particular solution.
dynamical evolutions where the dimensionless anisotropies, χ and χ do not diverge as we trace back in time to the moment of collision τ = 0. Furthermore, where See [15] for similar analytic behavior of the attractor.In Fig. 1 we show the evolution of the anisotropies in the different sectors and the total system for one particular attractor solution and neighboring trajectories (from [37], where more plots are given for this particular solution).
In order to determine the attractor curves on the attractor surface we proceed as follows.We initialize at a time τ 1 and choose our initial energy densities and non-equilibrium pressures, (ε, ε, ϕ, φ)(τ 1 ).We then solve the coupling equations (2.22)-(2.27) at the initial time to determine the metric components at the initial time.Subsequently, we evolve the system backwards in time via (2.18), (2.19), (2.20) and (2.21), and tune the values of (ϕ, φ)(τ 1 ) with fixed (ε, ε)(τ 1 ) such that χ and χ tend to constants for early times following (2.40).Having thus found an attractor curve to a satisfying level of precision (in this paper, 10 −15 ), we then evolve the system forward in time to study the late time approach to hydrodynamics.

Mapping the attractor surface at early and late times
Although in the application to Pb-Pb and p-p collisions, only times τ ≳ γ 1/4 are relevant phenomenologically, we shall begin by mapping out the attractor surface for arbitrary times.Even when the early-time behavior of the attractor surface is an insufficient approximation to a full semi-holographic treatment (let alone actual high-energy collisions), its constituent lines can be usefully parametrized by its parameters at some early reference time τ 0 ≪ γ 1/4 .

Early time behavior: dominance of the weakly coupled sector
The attractor surface has a universal scaling behavior in the limit τ → 0 which is determined in terms of the two dimensionless viscous parameters σ 1 and σ 2 defined in (2.41).Crucially, we need The relations σ 1,2 < 1 are imposed by causality [46].The two relations above imply that σ 2 > σ 1 as should be the case if the second subsystem is more strongly coupled than the first.The stricter lower bound on σ 2 is necessary for reasons to be described below.
We first note that with the above conditions, the anisotropies, χ and χ, go to σ 1 and σ 2 respectively at τ = 0 on the attractor surface, as we see in Fig. 2.Moreover, we find numerically that ε goes to zero while ε goes to a constant, namely γ −1 √ r − 1/ √ r.With these numerical inputs, assuming (3.1), and by expanding the conservation, MIS and coupling equations in τ around τ = 0, we readily determine the following universal early time scaling behavior of the energy densities and the anisotropies: and the following scaling behaviors of the components of the effective metrics: 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: As a result, in the limit τ → 0, the following five identities should hold on any curve on the attractor surface.Fig. 3 demonstrates the numerical check of the above identities on an attractor curve.
The stricter lower bound on σ 2 in Eq. (3.1) is necessary because the above scaling relations follow from solving the coupling equations, etc., with the assumption that ã vanishes as τ → 0. From the explicit scaling of ã in (3.3), it is clear that the lower bound is then self-consistent.Furthermore, it is also clear from (3.2) that χ also does not diverge as as a consequence of the lower bound on σ 2 and σ 1 < 1.
If instead of Eq. (3.1) we impose 0 < σ 2 < 1 and max 0, we obtain that ã diverges as τ (5/3)σ 1 −(4/3)σ 2 −(1/3) as τ → 0 and all other scaling relations discussed above are also altered.Nevertheless, χ and χ still reaches σ 1 and σ 2 as τ → 0 on the attractor surface.However, the above condition involves a lower bound on σ 1 which cannot be motivated otherwise, while (3.1) can be readily satisfied with reasonable weak coupling and strong coupling parameters.Therefore, we do not further investigate the above condition.
10 -6 10 -5 10 -4 0.001 0.010 0.100 1 0.7 0.8 0.9 Figure 3. Check of the identities (3.5) at early times.Similar identities given in (C.5) also holds for the disperser surface, and can be verified to a very high accuracy numerically as evident from Fig.It is easy to see from the above that the physical energy densities of the two subsystems behave as follows at early time and therefore Since σ 2 > σ 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 on the attractor surface which is reminiscent of the so-called bottom-up thermalization scenario in heavy-ion collisions [38]. 5Also note that E 1 always diverges as τ → 0 since σ 1 < 1.However, E 2 may or may not diverge in this limit, without contradicting any of the two conditions in (3.1), depending on the choice of parameters.While all initial conditions evolve towards specific curves on the attractor surface, they also get repelled from a disperser surface in which χ and χ approach −σ 1 and −σ 2 respectively in the limit τ → 0. On the disperser surface, the opposite of the energy distribution in bottom-up thermalization takes place, since both E 2 and E 2 /E 1 diverge as τ → 0. In contrast to (3.2), γ ε goes to the constant √ r − 1/ √ r while γε vanishes as τ → 0 on the disperser surface.For more details, see Appendix C. One can also find two other such half-disperser surfaces, in which χ and χ approach σ 1 and −σ 2 respectively, or χ and χ approach −σ 1 and σ 2 respectively as τ → 0. These can be analyzed in a similar manner.All these three surfaces are not relevant for understanding hydrodynamization with generic initial conditions.
Crucially, γε(τ 0 ) is bounded from above for any choice of ε(τ 0 )/ε(τ 0 ), as otherwise we do not obtain solutions.However, the energy densities and the total energy density can be arbitrarily large.This is similar to the equilibrium solutions studied in [33] where ε and ε were bounded from above, but the total energy density and pressure were unbounded.
We can systematically correct this error by iterative fine-tuning.At each step of the iteration, we fine tune ϕ and φ for each γε(τ 0 ) and ε(τ 0 )/ε(τ 0 ) such that χ = σ 1 and χ = σ 2 at progressively earlier times.Practically, the very first step of iteration gives a very precise estimate of initial values of ϕ and φ on the attractor surface for each γε(τ 0 ) and ε(τ 0 )/ε(τ 0 ), and also the upper bound on γε(τ 0 ) for each ε(τ 0 )/ε(τ 0 ) if we initialize at τ 0 = 10 −3 γ 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.

Hydrodynamic expansion
At late time, the full system can be mapped to a hydrodynamic derivative expansion with powers of τ −2/3 counting number of derivatives.The hydrodynamic expansion is determined only by two dimensionless parameters, which could be chosen for instance as lim τ →∞ γ 2/3 ετ 4/3 and lim τ →∞ γ 2/3 ετ 4/3 . (3.9) (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 γε(τ 0 ) and ε(τ 0 )/ε(τ 0 ) giving the initialization at time τ 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 τ dependence of ε and ε at late time as follows with • • • denoting viscous corrections and τ 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 τ 0 above is arbitrary.A different choice τ0 simply rescales ϵ 10 and ϵ 20 as follows: For the sake of convenience, we define dimensionless time: Note however that ω 1 and ω 2 are not independent time variables.In fact The hydrodynamic expansion in τ , 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 The above equations show all terms up to second order in derivatives.Defining we can invert (3.17) In terms of w 1 we can readily express the hydrodynamic expansion as a phase space curve With w 1 , w 2 , χ and χ as the phase space variables, the phase space curves representing the hydrodynamic expansion (or the attractor) are then given by the above expressions, χ(w 1 ), w 2 (w 1 ) and χ(w 1 ).As manifest from above, each curve is labelled by two dimensionless, time-reparametrisation invariant parameters, α and β.
The hydrodynamic late-time expansion of the components of the effective metrics are: Remarkably, up to O(w −4 1 ) the volume factors are: We can readily compute the full system energy density which behaves as where 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 α and β.
We can explicitly map the full system to a single expanding fluid at late time and compare to the above discussion.For a single fluid, the evolution of the energy density in the Milne background metric takes the generic form (up to first order) where ζ and η are the bulk and shear viscosity, respectively.Furthermore, the trace of the (first-order) hydrodynamic energy-momentum tensor takes the form We can assume that with k, κ and κ constants.Solving (3.33) with inputs from (3.35) in the large proper time expansion, we obtain and also (3.34) We can determine the constants k, κ and κ simply by matching the above expansions to those obtained by solving our system of equations.We readily see that these constants depend on the specific curve on the attractor surface.Comparing Eq. (3.37) with Eq. (3.31), we explicitly find that the effective viscosity of the full system is which depends manifestly on which curve on the attractor surface the full system evolves to at late time.In fact κ can be obtained from (3.39): Similarly, the hydrodynamic expansion of our hybrid system yields that Tr(T ) = − 8τ We can immediately conclude that • and κ = 0, i.e. the effective bulk viscosity vanishes since the τ −2 term vanishes in (3.41).
Since k and κ both depend on α and β, 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 χ, χ, ε and ε.Perfect hydrodynamics begins when the logarithmic plots of ε(τ ) and ε(τ ) become parallel to each other with slope −4/3, i.e. when they both reach the perfect fluid hydrodynamic tail 1/τ 4/3 , which is evident in Fig. 5 for the first-order hydrodynamic tails of ε and ε for the first and second attractors, respectively.First-order hydrodynamics begins when the logarithmic plots of χ(τ ) and χ(τ ) become parallel to each other with slope −2, i.e. when they both reach the hydrodynamic tail 1/τ 2 given by first-order hydrodynamics (note in the perfect fluid limit χ vanishes).This is shown in Fig. 5 for the first-order hydrodynamic tails of χ and χ for the first and second attractors, respectively.We can also match using leading order hydrodynamic expansion of χ ( χ): to a high accuracy.
In Appendix D we illustrate how to understand the hydrodynamic gradient expansion and see how its Borel resummation can describe the attractor surface in terms of the parameters α and β, defined in (3.12) and (3.13).

Late time approach to attractor: fluctuation analysis
The approach to the attractor at late time can be simply understood by perturbing around the late-time hydrodynamic expansion computed above.
The leading order linear fluctuations of energy densities, anisotropies and effective metric variables about the hydrodynamic expansion, i.e. the attractor solution at late time are defined as follows: and similarly for the tilded sector.Above τ := τ /τ 0 and δ ε(τ ), δϕ(τ ), etc. denote the fluctuations about the hydrodynamic expansion.We have truncated to the first few terms of the hydrodynamic expansions which can be found explicitly in Appendix A, because only these are sufficient to obtain the leading order behavior of the fluctuations.The fluctuations have two independent non-hydrodynamic modes capturing the missing information of the initial conditions in the four-dimensional phase space, and two hydrodynamic modes which span the two-dimensional attractor surface.Since the non-hydrodynamic modes, which do not admit a power series expansion τ −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 (δε(τ ), δ ε(τ ), δϕ(τ ), δ φ(τ )) by using six metric coupling equations.This leaves us with four first-order coupled ordinary differential equations in four variables, which are given in Appendix A. One can further solve these equations until one is left with a pair of decoupled fourth-order differential equations for the anisotropy fluctuations.Depending on the relative values of parameters, the fluctuations, δϕ(τ ) and δ φ(τ ), have different leading order forms, which we detail below.
In the case when C τ ϵ 20 1/4 ̸ = Cτ ϵ 10 1/4 , the two independent non-hydrodynamic solutions (at leading order) are: where m 1 , m 2 , n 1 , n 2 are constants, of which only two are independent: (3.54) will be the dominant exponential and give the relaxation towards the attractor.The powers of τ accompanying this exponential in (3.50) and (3.51) imply that δ φ(τ ) will decay faster towards the attractor surface than δϕ(τ ) in this case.This implies that the strongly coupled sector will hydrodynamize faster than the weakly coupled one.If (3.56) is the dominant mode.By an analogous argument, δϕ(τ ) will decay faster than δ φ(τ ) implying that the weakly coupled sector will hydrodynamize faster than the strongly coupled one.
Full numerical solutions do confirm that indeed when (3.53) is realized, the strongly coupled sector hydrodynamizes faster and the reverse occurs when (3.55) holds.However, realizing (3.55) requires choosing large ratios of ε(τ 0 )/ε(τ 0 ).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 ε goes to a constant while ε vanishes as τ → 0.
the two independent solutions to leading order are: where m 1 , m 2 , n 1 , and n 2 are arbitrary constants only two of which are independent: (3.60)

Matching with late times
Now we are prepared to study how parameters setting initial conditions on the Bjorken flow attractor surface map to α and β, which are defined in (3.12) and (3.13), and which determine the hydrodynamic expansion at late time.In particular, we find that they scale as a power of the total energy density on the attractor surface.For r = 2, the scaling of α and β is found to be α ∼ (γE(τ = γ 1/4 )) 3/2 , (3.61) where E is the total initial energy density.Using the following definitions of the total entropy density and the individual subsystem entropy densities, which are valid when both subsystems reach the perfect fluid limit: we can readily see from the hydrodynamic expansions discussed above that the individual system entropy densities behave at late time as which means that the total entropy scales like We can compare this to the pocket formula in [47] giving the final state particle multiplicity as a power of the total initial energy density for single conformal fluid attractors via (Sτ ) hydro ∼ (eτ ) 2/3 0 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 when the total energy density is evaluated at time γ 1/4 .The pocket formula of [47] states that the total particle multiplicity per unit rapidity6 follows dN dη ≈ (Sτ ) hydro . (3.68) 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) [48].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.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, . 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 γε = 0.8, ε/ε = 40.The green dot denotes the value at τ 0 = 10 −3 γ 1/4 and the red point is at τ = 10 +3 γ 1/4 .as follows from (3.31) and (3.41).An illustrative example for the full time-evolution of the interaction energy is shown in Fig. 6.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.

Hydrodynamization on the attractor surface
Figure 5 shows the existence of hydrodynamic tails of the attractor at a sufficiently late time, giving evidence of hydrodynamization in the hybrid attractor.This preliminary observation motivates a systematic study of the hydrodynamization times of the two sectors on the attractor surface.

Hydrodynamization time
Most commonly, hydrodynamization time is defined as the time when both the longitudinal and transverse pressure can be described by the first-order hydrodynamics up to a certain accuracy [6].For the hybrid hydrodynamic attractor, we define the hydrodynamization time of the hard (more weakly coupled) sector via and similarly for the strongly coupled (tilde) sector, so that τ hd and τhd denote the hydrodynamization times of the hard and soft sectors, respectively.This implies that the ratio of the departure of the anisotropic pressures (ϕ and φ respectively) from that given by first order hydrodynamics (ϕ 1st and φ1st respectively) to the respective isotropic parts of the pressures for both sectors is less than 10 percent after the corresponding hydrodynamization times.In general, the hydrodynamization time of perturbative sector τ hd is larger than that of the holographic sector τhd as shown in Fig 7, where the blue curve corresponds to the soft sector and the red one to the hard sector.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γ −1 ⪅ E(τ hd ), E(τ hd )≪γ −1 (negative slope in the Fig 7), we have found that the hydrodynamization time of both the sectors follows conformal behaviours i.e.
where p and p are constants.We note from Fig. 7 that the range of the energy density at hydrodynamization time in the conformal window is slightly larger in the weakly selfinteracting subsector.The value of p, 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 reference time τ 0 = 10 −3 γ 1/4 , we find p ≈ 0.63.(4.4)This is the same value (within numerical accuracy) that is obtained in the case of conformal attractors.We find that the hydrodynamization time of the perturbative sector in this conformal region is universal as shown in Fig 8 and the same as in conformal attractors studied earlier in the literature.However for the holographic sector as in Fig 9, the value of the p increases from 0.64 to 2.1 as the ratio of the energy densities of the individual sectors at some early time τ 0 = 0.001γ 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.
The doubly logarithmic plot of the scaled hydrodynamization time of the weakly coupled (hard) sector vs the total energy density at the τ hd shows universality and indistinguishability from other conformal attractors, when the hard sector is dominant at a reference time (here τ 0 = 0.001γ 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.
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., p assumes the conformal value 0.63, thus showing universality as in Fig 10 and coinciding with the usual conformal attractors.Meanwhile from Fig 11, 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-dominant 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.
Next, when the energy densities at hydrodynamization times of both the sectors are close to γ −1 i.e., E(τ hd ), E(τ hd ) ∼ γ −1 the hydrodynamization times on the attractor surface takes minimum values as shown in Fig 7. 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 τhd , when the weakly coupled (hard) sector is dominant at an initial reference time (here τ 0 = 0.001γ 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 p for ratios from 1 : 1 to 1000 : 1 of the initial energy densities of the subsectors increases from 0.64, 0.76, 0.93, 1.07, 1.21, 1.76 to 2.1.These different values of p shows the dependency of τhd on the ratio ϵ(τ 0 ) : ε(τ 0 ). 10 -6 10 -5 10 -4 0.001 0.010 0.100 1 0.5 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 τhd when the soft sector is dominant at the reference time τ 0 = 0.001γ 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 p is 0.63, and is the same as in Fig. 8 corresponding to the cases where the hard sector was dominant at initial reference time.
When E(τ hd ) and E(τ hd ) are very large, the hydrodynamization times increase rapidly and do not show any scaling behaviour like that in conformal window.This corresponds to the positive slope region in Fig 7 .-29 - The doubly logarithmic plot of the scaled hydrodynamization time of the weakly coupled (hard) sector vs the total energy density at the τ hd when the soft sector is dominant at the reference time τ 0 = 0.001γ 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. 9, 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.

Rough phenomenological match to p-p and Pb-Pb collisions
In order to contemplate possible phenomenological implications of these features, we interpret attractor curves with lower total energy densities at τ ∼ γ 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 γ −1/4 with the saturation scale Q p-p s of p-p collisions, and we assume that E(τ ) ∝ Q 4 s at τ ∼ Q −1 s with large A-A systems having a larger Q s .In the IPsat and bCGC models studied in [41] 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 For Pb-Pb collisions, where A 1/3 ≈ 6, we take simply as a typical value and get that at τ ∼ Q −1 s,P b the total physical energy density is E Pb-Pb (τ ∼ 0.63) ∼ 6.25 E p-p (τ ∼ 1) in units where γ = 1.
The Tables 1 and 2 show the hydrodynamization times of a selection of attractor solutions of the accordingly defined p-p and Pb-Pb collisions for various ratios of energy densities of the two subsectors at the conveniently chosen early reference time τ 0 = 10 −3 γ 1/4 , which we use to parametrize the attractor surface.1. Hydrodynamization times for p-p initial conditions parametrized by ε(τ 0 )/ε(τ 0 ) at a reference time τ 0 = 10 −3 ≪ τ i ∼ 1 (in units where γ = 1).Here w = E 1/4 1 τ , where E 1 is defined in (3.7), w hd = w(τ hd ), and similarly for the second subsystem.Note that for ε(τ 0 )/ε(τ 0 ) < 1, the hydrodynamization times for the hard (soft) sector are typically larger (smaller).
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 τ hd and τ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 as listed in Tables 1 and 2. The total energy density at hydrodynamization time for Pb-Pb collisions is scaled by a factor of 1/5 and 1/2 to fit in the same plot with p-p collision.
From these tables one finds that in our p-p collisions E(τ hd ) and E(τ hd ) are smaller, while the hydrodynamization time of the perturbative sector is larger and that of the holographic sector is smaller compared to the Pb-Pb collisions.This behaviour is reflected in Fig 12, where the E(τ hd ) and E(τ 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 (maximal value of γE(τ hd ) is about 0.1 for p-p and 0.6 for Pb-Pb collisions) while that of the non-perturbative sector lies in the non-conformal region for both p-p and Pb-Pb collisions (minimal value of γE(τ hd ) is about 0.5 for p-p and 0.9 for 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 7 as mentioned below (note also from this figure that the conformal windows for the two sectors are slightly different from each other): • In case of the perturbative sector for both p-p and Pb-Pb collisions, the total energy densities at the hydrodynamization time, E(τ hd ), lies approximately within the conformal region denoted by the black and the green dots in the figure.In addition, since the E(τ hd ) is minimal in p-p collisions, the hydrodynamization time of the perturbative sector τ hd is sufficiently larger which is evident from the scaling behaviour of τ hd in the conformal window.
• For the non-perturbative sector, E(τ hd ) > γ −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, E(τ hd ), in p-p collisions, the hydrodynamization time τhd of the soft sector of p-p collision is smaller than that of Pb-Pb collisions, as indicated from the growth of τhd with E(τ hd ).
Note from the Tables 1 and 2 that for smaller ratios of the hard to soft energy densities at reference time (ε(τ 0 )/ε(τ 0 ) ≲ 1), the hydrodynamization times for hard/soft sectors are typically larger/smaller, in both p-p and Pb-Pb collisions.For details on the evolution of energy densities in these scenarios, see Appendix E.

Relevance of results for general initial conditions
The above study of hydrodynamization, in the phenomenological context of small and large system collisions, has been restricted to the hydrodynamic attractor surface.However, the attractor surface is not necessarily relevant at τ ∼ Q −1 s , where we use conditions such as (4.6) to match the evolution with a more microscopic description, such as the glasma effective theory, via the total energy density.A natural question then is whether our results for hydrodynamization in large vs small system collisions remain valid for generic evolutions that may be far away from the attractor surface at a glasma matching time τ i ∼ Q −1 s .We have explored this issue by studying evolutions in the context of both p-p and Pb-Pb type collisions where the anisotropies of both hard and soft sectors are far from the attractor values around τ i ∼ Q −1 s , retaining the constraints on the total energy densities used in Sec.4.2, which are E Pb-Pb (τ i ∼ 0.63) ∼ 6.25 and E p-p (τ i ∼ 1) ∼ 1 respectively (in units γ = 1), in order that we could match with glasma models at the time where the latter just cease to be applicable and need to be replaced by kinetic-like descriptions.(Recall that we have identified γ with Q −4 s for p-p collisions.)Numerically, this has been implemented by initializing with values different from known attractor solutions both before and after τ i ∼ Q −1 s , and evolving forwards/backwards in time.While it is hard to explore the full four-dimensional phase space, initializing with different anisotropies (i.e.ϕ and φ) away from the attractor surface typically produce maximal deviations from the attractor surface around τ i ∼ Q −1 s .As shown in Fig. 13, in the Pb-Pb collision scenario, the convergence to the attractor happens much earlier than hydrodynamization. 7The hydrodynamization times for the hard/soft sectors for Pb-Pb collision type scenarios are identical to that of the corresponding attractor curves to which the evolution converges.The same is also true for the hard sector in p-p collision type scenario as shown in Fig. 13.However, for the soft sector in p-p collisions, as shown in Fig. 13, the convergence to the attractor is simultaneous with hydrodynamization, and therefore the hydrodynamization times in off-attractor evolutions do not coincide with that of the corresponding attractor curve, and can be larger than the latter by typically up to 50 percent (gray column in Fig. 13(d)).
Nevertheless, we note that even the largest hydrodynamization time for the off-attractor evolution in p-p type collision type scenario is smaller than the hydrodynamization time 7 Because of the nonlinearity of the dynamics, one cannot give a precise initial rate of approach to the attractor as this depends strongly on the specific solution.In the example shown in Fig. 13, the rate Γinit = |d ln(χ − χattractor)/dτ | at initial times γ −1/4 τi = 1 and 0.63, respectively, for p-p and Pb-Pb, is, however, always much larger than τ −1 hd , for both, the hard and the soft sectors.Close to the attractor, Γinit/τ −1 hd ∼ 60 for the hard sector in both p-p and Pb-Pb, and Γinit/τ −1 hd ∼ 20 (30) for the soft sector in p-p (Pb-Pb).Further away from the attractor, this ratio varies between 40 and 130 (20 and 50) for the hard (soft) sector of Pb-Pb, and between 10 and 70 (2 and 300) for the hard (soft) sector of p-p.   .Left: Generic evolutions of the dimensionless anisotropy which deviate from the attractor surface in Pb-Pb collisions for both hard sector (top) and soft sector (bottom).In all cases, the glasma matching condition E Pb-Pb (τ i ∼ 0.63) ∼ 6.25 (in units where γ = 1) is satisfied (see Sec. 4.2) at the glasma matching time τ i ∼ 0.63 = Q −1 s,P b which is indicated as a black dash-dotted vertical line.All evolutions converge to the corresponding attractor curve which is marked in bold red (for hard) and bold blue (for soft), and this occurs much earlier than hydrodynamization in both cases.The hydrodynamization time for the soft sector is marked as a vertical blue dotted line, while that for the hard sector is beyond the range of the plot.Right: The same for p-p collisions, where the glasma matching condition E p-p (τ i ∼ 1) ∼ 1 is satisfied (see Sec. 4.2) at the glasma matching time τ i ∼ 1 = Q −1 s,p indicated as a black dash-dotted vertical line.All evolutions converge to the corresponding attractor curve which is marked in bold red (for hard) and bold blue (for soft).In the case of the hard sector, the convergence to the attractor occurs also much before hydrodynamization time, but in the soft sector of p-p collisions, the convergence to the attractor is simultaneous with hydrodynamization.The hydrodynamization times for various evolutions are spread over a range marked as a gray column which includes the blue dotted vertical time corresponding to the attractor (the corresponding range in the Pb-Pb soft case is negligibly small).Note that the hydrodynamization time for the soft sector in Pb-Pb collisions is larger than the maximum value for that of the soft sector in p-p collisions.
in Pb-Pb scenario.The latter is determined by the attractor.Therefore, we conclude that the hydrodynamization of the soft sector in p-p type collisions occurs earlier than that in Pb-Pb type collisions in generic evolutions even if they deviate significantly from the attractor surface at τ ∼ γ 1/4 .
Thus the insights gleaned from the study of the attractor surface should apply also for more generic phenomenologically relevant initial conditions.Particularly, the facilitation of hydrodynamization of the soft sector due to smaller energy densities in small system collisions compared to that in large system collisions (since both the hydrodynamization times lie outside of the conformal window) should be valid for generic phenomenologically relevant initial conditions.

Discussion
The main physical insights which we have gained from the Bjorken flow attractor of our hybrid fluid model is how hydrodynamization may work in an asymptotically free gauge theory with both weakly interacting and strongly interacting degrees of freedom.In our model we represent those by two fluid components with viscosities differing by one order of magnitude, coupled in accordance with the semi-holographic framework developed in [28][29][30][31][32][33][34].This model has only one dimensionful energy scale, namely γ −1/4 , which sets the inter-system coupling.We have studied the hydrodynamic attractor of this system and we have shown analytically that at very early time, τ ≪ γ 1/4 , the energy density in the weakly coupled sector always dominates over that of the strongly coupled sector with an exponent which is determined by the transport coefficients and the relaxation times of the two sectors.This feature, which is reminiscent of the bottom-up thermalization scenario of Ref. [38], is, however, only relevant physically when the system is initialized correspondingly early; when in the context of heavy-ion collisions γ −1/4 is identified with the saturation scale Q s , no such hierarchy arises at τ ∼ Q −1 s , where we make contact with glasma descriptions, although the weakly coupled sector is typically the dominant one.
The hydrodynamic attractor of our hybrid fluid model is given by a (phase-space) surface of attractor lines that can be characterized by the ratio of the energy densities in the two sectors at some arbitrary early reference time together with the overall energy, and we have determined the resulting hydrodynamization times in the two sectors for a wide range of these parameters.
When the total energy densities at the respective hydrodynamization times (E(τ hd ), E(τ hd )) is less than γ −1 , both sectors hydrodynamize as in a conformal hydrodynamic attractor (τ hd ∝ E(τ hd ) −1/4 and τhd ∝ E(τ 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 being the 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 -the higher the dominance of the dominant sector, the more delayed is the hydrodynamization of the sub-dominant sector.When the energy densities at the respective hydrodynamization times increases beyond γ −1 , the hydrodynamization times of both sectors increase rapidly after reaching their minimum values.
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 subleading, the effective shear viscosity depends on how the full energy is finally shared between the two subsectors.
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 E(τ hd ) lies outside the conformal window, while that corresponding to the time of hydrodynamization of the perturbative sector, namely E(τ hd ), lies within the conformal window.It follows from our general results that while the perturbative sector hydrodynamizes later, the nonperturbative strongly interacting sector hydrodynamizes earlier in p-p collisions compared to Pb-Pb collisions on the attractor surface.This is because the hydrodynamization times for the soft sector in both cases lie outside of the conformal window, and the total energy density is smaller in the case of p-p collisions than in Pb-Pb collisions at the respective hydrodynamization times.Thus we get a new perspective on how collective flow develops in small system collisions.
Furthermore, we have shown that these results remain valid for generic phenomenologically relevant initial conditions with energy densities matched to the glasma models at time scales of order the saturation scale.In Pb-Pb collisions, the off-attractor evolutions for both sectors converge to the attractor surface much earlier than hydrodynamization, while in p-p collisions this happens only for the hard sector.In these cases, the attractor surface determines the hydrodynamization times.However, the approach to the attractor is simultaneous with hydrodynamization for the soft sector in p-p collsions, and the hydrodynamization times show typically up to 50 percent departure from the corresponding attractor value.Nevertheless, our conclusion that the hydrodynamization time of the soft sector in p-p collisions is smaller than in Pb-Pb collisions remains true even for evolutions which deviate significantly from the attractor surface at the time scales when we match with glasma models.
We therefore conclude that the hybrid attractor of our two component system with a more weakly and a more strongly interacting subsector is relevant for a range of phenomenologically relevant initial conditions well before the onset of hydrodynamics.The hybrid fluid model thus provides a potentially useful model for describing regimes with different importance of weakly and strongly coupled physics that will be interesting to confront with experiment.It would be interesting in particular to explore whether the different off-equilibrium evolution of the soft, more strongly interacting IR sector in small systems could lead to characteristic differences between p-p and Pb-Pb collisions.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, also for the study of dilepton production [48,49] and hadronization [47] 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, 8 and therefore such a study could be interesting for the search of a critical endpoint in heavy-ion collisions.
The equations of the fluctuations are given by δϕ

B Phase transition
As was explored in the inviscid case in [33], the coupling equations as presented exhibit a phase transition depending on the value of r.Thermal equilibrium is best explored in the Minkowski background.For the purposes of this section ε = 3P, P = nT 4  1 , and where the last relation is the requirement of thermal equilibrium, which means that the physical system temperature, T , can be expressed in terms of the individual subsystem temperature via Further, we introduce the light cone velocity v = a/b, and ṽ = ã/ b, (B.3) which represents the effective causal light cone of each individual subsystem.Interactions between the two systems change the effective light cone velocity.Note that c > v, ṽ > 0, where c = 1 is the speed of light.
light cone velocity in one system is multivalued, while in the other it is not.Of course, it is worthwhile to repeat that we are not in thermal equilibrium, which means we are not exactly mapping to the physical temperature T , but instead to an effective temperature.), as studied in [33].The blue line denotes the Bjorken evolution with increasing τ going from the green to red dot (right to left on the plot).Left: in the case r = 2, the green dot is for γε(τ 0 = 0.0001) = 0.348, while the red dot denotes τ = 20, the time we choose to stop the simulation.Right: below the phase transition, for r = 1.01 the Bjorken system is initialized at the green dot with γε(τ 0 = 0.01) = 0.1.14.As before, the evolution in the τ coordinate follows from the green to red dot (right to left).The upper solid purple line denotes the more viscous system, while the lower solid blue line is the weakly coupled system.The initial conditions are identical for both systems for each r and with initially no dissipation ϕ(τ 0 ) = φ(τ 0 ) = 0 at τ 0 = 1, i.e. top: r = 2 and γε = γ ε = 0.348, middle: r = r c with γε = γ ε = 0.248, and finally bottom: r = 1.2 with γε = γ ε = 0.3, while the transport coefficients are given by (2.30).
The five identities analogous to those in Eq.We also find that the other six coefficients of the early time expansion can be determined in terms of k 1 and g 1 in consistency with the fact that the disperser surface is two-dimensional.
In Fig. 16, the early time behaviors of χ, χ, ε and ε in a disperser solution have been shown, and in Fig. 17

D 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, α and β (defined in (3.12) and (3.13)) 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 = Let's use the notation ξ i with i = 1, 2 and χ 1 = χ and χ 2 = χ.At arbitrarily high order the coefficients of these two series expansions r i,n for χ i show factorial growth following indicating the series is asymptotic in nature with zero radius of convergence, which is seen in subsystems in the top panels of Fig. 18.The parameters ξ i,0 and γ i for the weak and strong systems depend on α and β.
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 [54].The Borel transform of the series is given by r i,n are the coefficients from the series (D.1).This Borel transformed function has a finite radius of convergence about the origin in the complex ξ 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 χ B,i (ξ).Since we only have access to a finite number of terms, we will need to approximate the integrand, which we do via Padé approximant χ B,i (ξ) ≈ P N,M,i (ξ) = N i=0 n i ξ i 1 + M j=1 m j ξ j (D.3) where the coefficients of the Padé approximant are fixed to agree with series coefficients (D.2).The Padé series can go up to arbitrary orders with a constraint N + M = k, where k is the order up to which Borel transformed anisotropy is truncated.Most commonly, one considers the symmetric case with N = M = k/2.In our case we have taken the order k to be as high as 500 with precision set to thousand digits.The singularities of the function χ B,i (ξ) appear as concentrated poles of the Padéapproximant in the Borel plane and they depend on the dimensionless parameter α and β.In Fig. 18, we show the pole structure of both sectors in the Borel plane.The poles start from a point nearest to the origin, namely ξ 1,0 and ξ 2,0 for χ B (ξ) and χB (ξ) respectively, and accumulate along the real axis giving an image of branch cut as in the case of decoupled MIS [12].The starting value of the branch cut is proportional to the inverse of the slope of the coefficients shown in Fig. 18 and can be identified with the non-hydrodynamic mode that decays exponentially at late-times.From the analysis in Sec.3.4 we readily see that  [12,16,54] χ i (x) = 1 x C e −ξ/x P N,M,i (ξ)dξ (D.4) where C is the contour from 0 to ∞. However the presence of singularity in χ B,i (ξ) introduces complex ambiguities implying that the integration depends on the choice of contour 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 χ and χ, respectively, as shown in Fig. 18.We have two Stokes coefficients corresponding to the two exponential decay modes discussed in Sec.3.4 and these give the two missing data (other than the parameters α and β appearing in the hydrodynamic expansion) which specify the initial conditions fully in the four-dimensional phase space.

E Different initial ratios of subsector energy densities
While at sufficiently early times all attractor solutions have a dominant weakly coupled component, different initial conditions lead to different ratios of energy densities at some initial reference time, chosen here as τ 0 = 0.001γ 1/4 .In figure 19 we display the evolution of the energy densities in the subsectors and the total energy density for different initial ratios, comparing the two sets of initial conditions we tentatively associated with p-p and Pb-Pb collisions in sect.4.2.For a ratio E 1 : E 2 less than 4:1, the energy density in the holographic sector overtakes the one in the perturbative sector for a stretch of time, whereas in the opposite case it never gains dominance.At sufficiently large time, the "perturbative" sector becomes again dominant in each case, but although this is reminiscent of the transition to a weakly coupled hadron gas after hadronization, we consider our model as a toy of heavy-ion collisions only in the deconfined early stages.

Figure 1
Figure1.Anisotropies χ 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[37], where more plots are given for this particular solution.

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

γ - 1 Figure 7 .
Figure 7.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 τ hd and τhd for ratio ϵ(τ 0 ) : ε(τ 0 ) = 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.

Figure 9 .
Figure 9.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 τhd , when the weakly coupled (hard) sector is dominant at an initial reference time (here τ 0 = 0.001γ 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 p for ratios from 1 : 1 to 1000 : 1 of the initial energy densities of the subsectors increases from 0.64, 0.76, 0.93, 1.07, 1.21, 1.76 to 2.1.These different values of p shows the dependency of τhd on the ratio ϵ(τ 0 ) : ε(τ 0 ).

Figure 10 .
Figure 10.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 τhd when the soft sector is dominant at the reference time τ 0 = 0.001γ 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 p is 0.63, and is the same as in Fig.8corresponding to the cases where the hard sector was dominant at initial reference time. p

Figure 13
Figure 13.Left: Generic evolutions of the dimensionless anisotropy which deviate from the attractor surface in Pb-Pb collisions for both hard sector (top) and soft sector (bottom).In all cases, the glasma matching condition E Pb-Pb (τ i ∼ 0.63) ∼ 6.25 (in units where γ = 1) is satisfied (see Sec. 4.2) at the glasma matching time τ i ∼ 0.63 = Q −1 s,P b which is indicated as a black dash-dotted vertical line.All evolutions converge to the corresponding attractor curve which is marked in bold red (for hard) and bold blue (for soft), and this occurs much earlier than hydrodynamization in both cases.The hydrodynamization time for the soft sector is marked as a vertical blue dotted line, while that for the hard sector is beyond the range of the plot.Right: The same for p-p collisions, where the glasma matching condition E p-p (τ i ∼ 1) ∼ 1 is satisfied (see Sec. 4.2) at the glasma matching time τ i ∼ 1 = Q −1 s,p indicated as a black dash-dotted vertical line.All evolutions converge to the corresponding attractor curve which is marked in bold red (for hard) and bold blue (for soft).In the case of the hard sector, the convergence to the attractor occurs also much before hydrodynamization time, but in the soft sector of p-p collisions, the convergence to the attractor is simultaneous with hydrodynamization.The hydrodynamization times for various evolutions are spread over a range marked as a gray column which includes the blue dotted vertical time corresponding to the attractor (the corresponding range in the Pb-Pb soft case is negligibly small).Note that the hydrodynamization time for the soft sector in Pb-Pb collisions is larger than the maximum value for that of the soft sector in p-p collisions.

Figure 14 .
Figure 14.Light cone velocity as a function of temperature for the inviscid Bjorken and Minkowski cases of identical subsystems for two different values of r.The orange dashed curve is the result of working directly in the Minkowski background (B.4), as studied in[33].The blue line denotes the Bjorken evolution with increasing τ going from the green to red dot (right to left on the plot).Left: in the case r = 2, the green dot is for γε(τ 0 = 0.0001) = 0.348, while the red dot denotes τ = 20, the time we choose to stop the simulation.Right: below the phase transition, for r = 1.01 the Bjorken system is initialized at the green dot with γε(τ 0 = 0.01) = 0.1.

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