Homogeneous isotropization and equilibration of a strongly coupled plasma with a critical point

We use holography to investigate the process of homogeneous isotropization and thermalization in a strongly coupled $\mathcal{N} = 4$ Super Yang-Mills plasma charged under a $U(1)$ subgroup of the global $SU(4)$ R-symmetry which features a critical point in its phase diagram. Isotropization dynamics at late times is affected by the critical point in agreement with the behavior of the characteristic relaxation time extracted from the analysis of the lowest non-hydrodynamic quasinormal mode in the $SO(3)$ quintuplet (external scalar) channel of the theory. In particular, the isotropization time may decrease or increase as the chemical potential increases depending on whether one is far or close enough to the critical point, respectively. On the other hand, the thermalization time associated with the equilibration of the scalar condensate, which happens only after the system has relaxed to a (nearly) isotropic state, is found to always increase with chemical potential in agreement with the characteristic relaxation time associated to the lowest non-hydrodynamic quasinormal mode in the $SO(3)$ singlet (dilaton) channel. These conclusions about the late dynamics of the system are robust in the sense that they hold for different initial conditions seeding the time evolution of the far-from-equilibrium plasma.


Introduction and summary
With the advent of the holographic gauge/gravity duality [1][2][3][4], it has become possible to study physical properties of some strongly coupled quantum systems using classical gravity in higher dimensions (for reviews see, e.g., [5][6][7]). Concerning strongly correlated quantum fluids, this framework has made it possible the investigation of several aspects of such systems, such as their thermodynamics and hydrodynamics [8][9][10], quasinormal modes [11][12][13], and also the far-from-equilibrium dynamics describing the relaxation of holographic fluids toward thermodynamic equilibrium in many different settings (see e.g. [14,15] for recent reviews). One of the main attractive features of holography is its unique ability to deal with the entire evolution of a strongly coupled fluid within a single framework, starting from farfrom-equilibrium anisotropic initial states and dynamically evolving them passing through different stages comprising several kinds of characteristic "relaxation times" until reaching thermodynamic equilibrium. These different relaxation times characterize, for instance, the onset of applicability of hydrodynamics (known as the "hydrodynamization time", which also depends on the specific formulation of hydrodynamics considered), the onset of nearly isotropization of the system (when the longitudinal and transverse pressures in a given flow are approximately equal), the onset of applicability of the equilibrium equation of state in nonconformal plasmas (known as the "EoSization time" [16]), and the onset of true thermalization (when all the physical observables of the theory have approximately relaxed to their equilibrium values). In fact, one of the main outcomes of holographic investigations of far-from-equilibrium strongly coupled quantum fluids was the conclusion that in some cases the system may hydrodynamize when the fluid is still significantly anisotropic and farfrom-equilibrium [17][18][19][20]. This finding, together with recent calculations [13,[21][22][23] that showed that the gradient expansion diverges, have significantly changed our understanding of relativistic hydrodynamics (for a review see [24]).
The fast expanding fireball produced in ultrarelativistic heavy ion collisions at RHIC [25][26][27][28] and LHC [29] is probably the most remarkable example of a dynamical system actually realized in nature featuring a very rich and complex time evolution characterized by both hard (i.e., perturbative) and soft (nonperturbative) physics. Just before the heavy ions collide, the gluon density inside these large nuclei at very high energies is expected to saturate producing the so-called color glass condensate (CGC) [30][31][32][33], which has become the starting point for the initial conditions in heavy ion collisions. Right after the collision the medium has an enormous amount of energy concentrated in a very small volume, which starts to rapidly expand passing through different stages. Before 1 fm/c after the collision the system is expected to be a highly dense coherent medium dominated by the dynamics of classical QCD fields called "glasma" [34]. 1 As the system keeps expanding, the glasma decoheres towards a new state of QCD matter called the quark-gluon plasma (QGP) [35][36][37], whose relevant degrees of freedom correspond to deconfined quarks and gluons.
In practice, the QGP produced in heavy ion collisions is well described by viscous hydrodynamics with an equilibrium equation of state and transport coefficients which are compatible with soft physics expectations [38][39][40][41], indicating that at this stage the system is strongly coupled, contrary to the scenario just after the collision where weak coupling physics plays a prominent role. As the QGP keeps expanding and cooling down it eventually enters in the crossover region [42,43] and hadronizes, giving place to a hadron gas. Later stages of the temporal evolution of heavy ion collisions include the regions of chemical freeze-out 2 and the thermal or kinetic freeze-out 3 . After this stage, the produced hadrons are essentially free and the yields of their decays reach the detectors of the experimental apparatus providing a large amount of information about the evolution of the system.
The full dynamical evolution of high energy heavy ion collisions outlined above cannot be completely described by the gauge/gravity duality since the former encompasses both hard and soft physics while the latter can only deal with strongly coupled systems. In fact, it is well-known that the asymptotically free ultraviolet regime of QCD cannot be described by the gauge/gravity duality given that holographic models are usually characterized by strongly coupled ultraviolet fixed points. On the other hand, even though a rigorous top-down construction of a gravity dual of the hydrodynamized and strongly coupled QGP is missing, there are phenomenological bottom-up Einstein-Maxwell-Dilaton (EMD) gauge/gravity constructions which have been shown to successfully describe in practice, not only on a qualitative level but also quantitatively, a plethora of thermodynamical and hydrodynamical observables characterizing the physics of the strongly coupled QGP under many different situations [44][45][46][47][48][49][50][51][52][53][54][55][56]. Therefore, one striking question which is posed in face of the above considerations is the following: in which cases could one hope, at least in principle, to obtain possible useful insights (or even some quantitatively accurate results) for the far-from-equilibrium dynamics of heavy ion collisions using holographic techniques?
The answer to the question above is still unsettled but one may gauge the applicability of gauge/gravity models to the analysis of heavy ion collisions by looking at some recent results comparing lattice QCD calculations (which are performed in equilibrium) with heavy ion experimental data. This kind of comparison may help to identify under which conditions the expanding fireball would probe, through more stages, a strongly coupled regime of QCD.
For instance, in Ref. [41] it was shown that in central heavy ion collisions at RHIC and LHC the experimentally extracted ratio between the pressure and the internal energy density of the QGP, (p/ϵ) exp (T eff ) = 0.21 ± 0.10, is in good agreement with the corresponding lattice QCD estimate, (p/ϵ) lQCD (T eff ) ≈ 0.23, where T eff is half-way between the effective temperatures associated with the collision energies √ s = 200 GeV and √ s = 2.76 TeV (at these high energies, the baryon chemical potential is negligible compared to the temperature of the medium). Moreover, in Ref. [40] a state-of-the-art Bayesian analysis was simultaneously applied to several physical observables while varying the free parameters 2 When inelastic collisions between the produced hadrons cease and the relative ratio between the different kinds of particles remains fixed. 3 When the average distance between the hadrons is large enough to make the nuclear interaction between them effectively negligible, causing the momentum distribution of the particles to remain fixed.
of the model. The results in the space of parameters of the model which match combined heavy ion data measured at RHIC and LHC within the interval √ s = 200 GeV -2.76 TeV were also found to be consistent with results from lattice QCD simulations. These findings show that the QGP produced in these collisions at high energies is not only adequately described by hydrodynamics but also that the (equilibrium) QCD equation of state obtained in lattice simulations may be trustfully used in hydrodynamic simulations of the spacetime evolution of the system at these high energies (i.e., √ s ≳ 200 GeV). Furthermore, it was also shown in Refs. [57,58] that ratios between higher order susceptibilities of baryon and electric charges calculated in equilibrium on the lattice give a good description of experimentally measured ratios between moments of net-proton and net-electric charge multiplicity distributions for collision energies √ s ≥ 39 GeV; however, for √ s < 39 GeV the compatibility observed at higher energies between the set of chemical freeze-out parameters extracted from the independent analysis performed in the baryon and electric charge sectors on the lattice is not warranted [58]. This suggests that at chemical freeze-out the system produced in heavy ion collisions is less close to equilibrium if the system has a higher chemical potential (corresponding to a lower collision energy).
Recently, using a phenomenological holographic EMD model which quantitatively reproduces lattice QCD thermodynamics at zero and nonzero baryon chemical potential [49,52], it was found that an increase in the baryon chemical potential decreases the shear viscosity times temperature to enthalpy density ratio, ηT /(ϵ + p), which gives a measure of the fluidity of the medium [59] 4 . This indicates that the QGP at finite baryon density remains strongly coupled. From the above considerations, one may expect that the gauge/gravity duality is more likely to produce useful insights for the far-from-equilibrium dynamics of heavy ion collisions when the system is doped with a nonzero chemical potential 5 . This is exactly the scenario we are interested in analyzing in the present work.
The literature regarding out-of-equilibrium holographic dynamics is vast, and even if we restrict ourselves to works that deal with models endowed with a chemical potential there is a substantial amount of papers already available. For non-equilibrium aspects of some strongly coupled plasmas doped with a chemical potential which do not employ numerical relativity and probe the thermalization process using non-local observables (e.g. Wilson's loops), see e.g. Refs. [61][62][63][64]. In the far-from-equilibrium context, there is the study of Ref. [65] regarding homogeneous equilibration in a charged plasma without a critical point, while Ref. [60] considers a shock wave analysis with baryon charge, yet without a critical point. One can also study quantum critical points in the context of holographic quenches, as done in Refs. [66,67]. Recent studies about holographic isotropization in the context of Gauss-Bonnet gravity can be found in Ref. [68].
One of the main goals of the present work is to assess the impact of a critical point in the far-from-equilibrium dynamics of a relativistic strongly coupled plasma at finite density, which to the best of our knowledge is a question that has not been previously studied in the literature. In fact, this is the simplest problem one may consider that can lead to interesting insights into the strongly coupled dynamics of the QGP near a critical point at finite density, one of the focus of RHIC's Beam Energy Scan (BES) program. Even though the holographic model we consider here is very different from QCD, we believe such a study is important given that there are no other approaches that can be used to perform real-time, far-from-equilibrium calculations at strong coupling near a critical point.
More specifically, we shall investigate here the homogeneous equilibration dynamics of a far-from-equilibrium top-down holographic plasma at finite density known as the 1-R charge black hole (1RCBH) [69][70][71][72][73][74]. This model describes a strongly coupled conformal N = 4 Super Yang-Mills (SYM) plasma charged under a U (1) subgroup of the global SU (4) R-symmetry and features a critical point (CP) in its phase diagram. Our gravitational setup is different from the one considered in Ref. [65] where the authors studied the homogeneous equilibration of a strongly coupled SYM plasma with a nonzero charge density without a CP. Indeed, as we are going to see in a moment, the action of the 1RCBH model includes a dilaton field that considerably changes the dynamics of the theory and leads to a CP in its phase diagram.
In order to solve the corresponding numerical relativity problem, in this work we follow the pioneering work of Chesler and Yaffe [75] by employing the so-called characteristic formulation of general relativity [76], which is very convenient for asymptotically AdS spacetimes (see, e.g., Ref. [14] for a review). Alternatively, one may also use the ADM formulation [77,78] of numerical relativity in asymptotically AdS spacetimes, as done in Refs. [79,80].
An important remark must be done at this point. As it is well-known, the thermodynamics and the hydrodynamics of the SYM plasma are very different than what is observed in the QGP produced in heavy ion collisions (see, e.g., [81]). Thus, one should not expect to obtain, in general, quantitative insights for the early time dynamics of heavy ion collisions from the study of the far-from-equilibrium dynamics of the 1RCBH model. However, it may be that some qualitative properties derived in the far-from-equilibrium dynamics of this model are robust or "nearly universal", as it happens with the shear viscosity to entropy density ratio of holographic fluids where η/s = 1/4π for a broad class of strongly coupled systems [8]. In order to look for possible signatures of this kind of "universality", one should also consider the far-from-equilibrium dynamics in other holographic models at finite density endowed with a CP.
For instance, the quasinormal mode (QNM) oscillations of the 1RCBH model in the external scalar and vector diffusion channels were found in Ref. [82] to be damped as one increases the U (1) R-charge chemical potential, as long as one is away from criticality. This finding is in qualitative accordance with the observation done in Ref. [49] that by increasing the baryon chemical potential far from the CP there is a damping in the QNM oscillations of the external scalar channel of a phenomenological bottom-up EMD model at finite baryon density that quantitatively describes QGP thermodynamics. This damping of the QNM oscillations caused by increasing the chemical potential away from criticality in two very different holographic models may be a signature of robustness or an indication of a general behavior for strongly coupled systems at finite density.
However, close to the CP, the QNM oscillations of the 1RCBH model in the external scalar and vector diffusion channels revert this trend observed at lower densities and increase as the U (1) R-charge chemical potential is enhanced [82]. On the other hand, the behavior of the QNMs of the QCD-like phenomenological EMD model at finite baryon density of Refs. [49,53] has not yet been investigated close to criticality, a task we postpone for a future work. It would be interesting to check if the same qualitative modification in the QNM oscillations of the 1RCBH model caused by the proximity of the CP is also featured in the phenomenological EMD model, since this could point out to a possible universal effect of the CP on the equilibration dynamics of strongly coupled fluids.
Therefore, even though the 1RCBH model (and the SYM plasma in general) is not a holographic setup suited for direct applications to heavy ion phenomenology, it may potentially contain some general properties displayed by strongly coupled fluids driven out-of-equilibrium. Moreover, one of the main conclusions of the present work will be the distinction of two characteristic equilibration times of the far-from-equilibrium system. Namely, by looking at the imaginary part of the lowest non-hydrodynamic QNMs of the model one may extract an upper bound for some characteristic "relaxation times" of the theory, according to the general reasoning first devised in Ref. [11]. The aforementioned behavior of the QNMs in the external scalar and vector diffusion channels would, therefore, suggest that the "equilibration time" of the 1RCBH model decreases with increasing chemical potential far from the CP, while close to the CP it would instead increase. As we shall see in this work, this is, in fact, the behavior found for the isotropization time of the system, which is dominated by the lowest non-hydrodynamic QNM of the external scalar channel of the theory.
On the other hand, even after (nearly) isotropization is reached, the scalar condensate dual to the bulk dilaton field may still be significantly far from its equilibrium value. Since in the present work the equilibration of the scalar condensate will always be the last equilibration time of the system, we shall associate it with the true thermalization time of the medium. Contrary to the isotropization time, in this model the thermalization time always increase with increasing chemical potential. As we are going to show, this is in consonance with the behavior of the lowest non-hydrodynamic QNM of the dilaton channel. Also, in the near future we intend to extend the present analysis to consider the case of the QCD-like EMD model of Ref. [53], which provided a prediction for the location of the long-sought critical point of the QCD phase diagram in the plane of temperature and baryon chemical potential.
The outline of this work is as follows. In Sec. 2 we present the equations of motion for the 1RCBH model assuming a time-dependent and spatially homogeneous anisotropic Ansatz for the bulk fields. This gives a set of nonlinear partial differential equations (PDEs) to be solved numerically. We also discuss the relevant observables of the dual quantum field theory (QFT) that we need to compute to analyze the homogeneous isotropization and thermalization processes in this setup, namely, the one-point functions associated with the expectation value of the boundary stress-energy tensor, ⟨T µν ⟩, dual to the bulk metric field g µν , the scalar condensate, ⟨O ϕ ⟩, dual to the bulk dilaton field ϕ, and the expectation value of the conserved U (1) R-current, ⟨J µ ⟩, dual to the bulk Maxwell field A µ . The derivation of the general form of these one-point functions via holographic renormalization is presented in Appendix A. In Sec. 3 we perform the near-boundary asymptotic expansion of the bulk fields in order to fix the boundary conditions for the set of PDEs. From this near-boundary analysis, we will be left with a couple of unknown time-dependent ultraviolet coefficients that shall be related with the one-point functions of the dual QFT. In Sec. 4, we briefly review the analytical equilibrium solutions of the 1RCBH model and its thermodynamics. In Sec. 5, we discuss some relevant technical issues related to the numerical time-dependent far-from-equilibrium solutions of the 1RCBH model, such as the choice of the initial conditions, and explain how we numerically solve the set of coupled nonlinear PDEs for the far-from-equilibrium system. Once the numerics is settled, we proceed in Sec. 6 to provide the main results of the homogeneous equilibration dynamics of the 1RCBH model. Moreover, in Subsection 6.7 we perform the match of the late time behavior of the pressure anisotropy and the scalar condensate with the lowest nonhydrodynamic QNMs of the external scalar and dilaton channels, respectively. While the QNMs of the external scalar channel were obtained in Ref. [82], we derive in Appendix B the QNMs of the dilaton channel. Finally, we close the paper in Sec. 7 with an outlook of our main results and also point out future perspectives and ongoing investigations.
Throughout this paper we work with natural units where ℏ = c = k B = 1 and use a mostly plus metric signature. Capital Latin indices {M, N, . . . } denote the coordinates {t, r, x, y, z} of a five dimensional asymptotically AdS spacetime where the gravity theory is defined, whereas Greek indices {µ, ν, . . . } represent the coordinates {t, x, y, z} of the four dimensional boundary.

The holographic model and its equations of motion
The 1RCBH model [69][70][71][72][73][74] first appeared as a solution of the five dimensional N = 8 gauged supergravity action [70], which was later demonstrated to lie within a class of solutions equivalent to near-extremal spinning D3-branes in AdS 5 ×S 5 [73]. The Kaluza-Klein compactification of the five sphere S 5 on the spinning D3-branes solutions leads to the SU (4) R-symmetry and the three independent Cartan subgroups of the R-Symmetry [70]. The general solution is known as the STU model, whilst the 1RCBH model is obtained by considering only one charge, i.e., Q ≡ Q a and Q b = Q c = 0. A thorough discussion of the matter content of the dual QFT at zero temperature may be found in Ref. [83]. For the thermal plasma at finite density, detailed discussions may be found in Refs. [47,82]. For the sake of completeness, in the present work we briefly review the thermodynamics of the 1RCBH plasma in Sec. 4.1.
The gravitational action of the 1RCBH model is given by [69][70][71][72][73][74], where κ 2 5 = 8πG 5 with G 5 being the five dimensional Newton's constant, F M N = ∂ M A N − ∂ N A M , with A M being the Maxwell gauge field, and is the Gibbons-Hawking-York boundary action [84,85] needed to provide a well-posed Dirichlet problem. In this term, γ µν denotes the induced metric at the boundary and K represents its extrinsic curvature. The last term in Eq. (2.1) is the boundary counterterm action needed to remove the divergences of the on-shell action. In Appendix A we give the details regarding the boundary terms and show how to obtain the desired one-point functions of the dual QFT. The expressions for the dilaton potential V (ϕ) and the Maxwell-Dilaton coupling f (ϕ) in the action (2.1) which define the top-down construction corresponding to the 1RCBH model are given by where L is the AdS radius, which is set to unity henceforth for simplicity. Moreover, by Taylor expanding the dilaton potential in powers of ϕ close to the boundary, we obtain which tells us that the mass of the dilaton is given by m 2 = −4. Recalling that the relation between the mass of the dilaton and the scaling dimension of its dual operator in the boundary QFT is given by m 2 = ∆(∆ − 4), one concludes that ∆ = 2. Note also that the dilaton field vanishes at the boundary such that the bulk geometry is asymptotically AdS 5 . Einstein's equations are obtained from the variation of the EMD action with respect to the metric, where, is the stress-energy tensor of the matter fields A µ and ϕ. It is usually simpler to work with the "trace-reversed" form of Einstein's equations. 6 This may be derived by noting that Eq. (2.5) implies that (2.7) 6 The tensorial algebra using a specific chart in this paper is done with the help of the Mathematica's package "Riemann Geometry and Tensor Calculus" (RGTC) [86].
By substituting Eq. (2.6) into Eq. (2.7), one rewrites Einstein's equations as follows On the other hand, Maxwell's equations are obtained by varying the EMD action with respect to the Maxwell field The last equation that one needs to fully specify the dynamics of the model is the Klein-Gordon equation for the dilaton field (2.10) To explore the far-from-equilibrium solutions of the model, we adopt the well-known characteristic formulation for asymptotically AdS spacetimes [75]. We consider here a time-dependent and spatially homogeneous anisotropic Ansatz for the metric field, which is suited to study homogeneous isotropization dynamics, in which one starts with an anisotropic configuration and an energy density that is conserved as time evolves. We work with generalized infalling Eddington-Finkelstein (EF) coordinates, in terms of which one may write the Ansatz for the line element as follows [75], where v represents the EF time, which near the boundary is interpreted as the time of the dual QFT. To see this, recall that the EF time is defined via dv = dt + − g rr g tt dr, (2.12) where g rr and g tt are the radial and temporal diagonal components of an asymptotically AdS 5 metric. Thus, as one goes to the boundary located at r → ∞, one obtains v → t. We also remark that in these generalized infalling EF coordinates, infalling radial null geodesics satisfy v = constant, while outgoing radial null geodesics satisfy dr/dv = A(v, r) [75] 7 . Furthermore, there is a residual diffeomorphism invariance in Eq. (2.11) corresponding to the radial shift r → r + λ(v), where λ(v) is an arbitrary function of the EF time [14]. With respect to the Maxwell and dilaton fields, the Ansätze are The U (1) R-charge chemical potential is associated with the boundary value of the Maxwell field Φ(v, r) in equilibrium, while the dilaton field vanishes at the boundary, as mentioned above.
The resulting equations of motion for the EMD system form a set of coupled nonlinear PDEs, where the prime denotes ∂ r , and  [65], we also defined a bulk "electric field", It is important to note that the Maxwell equation (2.14b) relates ϕ, Σ, and Φ ′ , i.e.
We will show in the next section how one can relate the above unknown constant to the U (1) R-charge density ρ c using results from holographic renormalization. Indeed, Eq.
(2.17) expresses the existence of a Gauss charge coming from the Gauss law of classical electromagnetism.
9. The constraint (2.14g) is useful for checking the accuracy of the numerical solutions.
Before we delve into the numerics and solve the system of PDEs (2.14a)-(2.14g) following the above algorithm, there are still some technical details that we need to take into account. Some of these details are: the boundary conditions, the equilibrium solutions, and the initial data.
In order to fix the boundary conditions for the EMD fields we first need to perform a near-boundary expansion of the bulk fields. This expansion will reveal which ultraviolet coefficients need to be dynamically fixed. Such analysis is done in Sec. 3. The detailed knowledge of the equilibrium solutions, which are discussed in Sec. 4, is very important to determine the final state given some initial geometry. Moreover, it is only possible to predict the equilibrium state because in this homogeneous setup the energy and charge density are constant. Still in Sec. 4, using the formulas derived from holographic renormalization in Appendix A, we provide holographic formulas to calculate the energy density and the pressure coming from ⟨T µν ⟩, the charge density that comes from the temporal component of ⟨J µ ⟩, and the expectation value of the scalar operator dual to the dilaton field, ⟨O ϕ ⟩. In Sec. 5, we discuss different initial data and specify the numerical techniques that we shall employ to solve the equations of motion of the 1RCBH model. After that, we will be ready to present in Sec. 6 our results for the homogeneous isotropization and thermalization of the 1RCBH plasma, emphasizing in particular the dynamics near the critical point.

Near-boundary expansion of the bulk fields
In a five dimensional model where the dilaton field has scaling dimension ∆ = 2 the nearboundary expansions of the bulk fields are given by integer powers of the radial coordinate plus logarithmic terms [87], i.e.
where λ(v) is the radial shift function associated with the aforementioned residual diffeomorphism invariance of the bulk metric [14]. As we shall see in Sec. A, the logarithmic terms in the above expansions vanish due to three main facts: the conformal flatness of the boundary, the scaling dimension of the QFT scalar operator dual to the bulk dilaton field, and the conformal symmetry of the system. Therefore, we can already set the logarithmic corrections to zero. We adopt then the following notation Substituting the expansions (3.1)-(3.5) into the equations of motion (2.14a)-(2.14g), setting the radial shift function λ(v) to zero, and eliminating all possible coefficients in favor of the others, the ultraviolet asymptotic behavior of the EMD fields reads where the dot represents the time derivative ∂ v . Furthermore, we find that this near-boundary analysis cannot determine five coefficients: With the exception of Φ 0 , these coefficients are dynamical, i.e. we need the bulk solution to determine them. The coefficient Φ 0 (v) is fixed by the Dirichlet boundary condition for the Maxwell field which imposes that its boundary value gives the chemical potential associated with the U (1) R-symmetry 9 (3.14) The coefficient H is defined by By working out the equations of motion up to O(r −3 ) in the near-boundary expansions of the bulk fields one concludes that H is a constant. Moreover, as we shall see in a moment, this coefficient H is, up to a numerical factor, the energy of the system which is conserved in this homogeneous setup. On the other hand, the coefficients B 4 (v) and ϕ 2 (v) are timedependent quantities related to the pressure anisotropy and the scalar condensate dual to the dilaton field, respectively. Additionally, the coefficient Φ 2 (v) is actually a constant. This can be shown by expanding the equations of motion near the boundary up to O(r −5 ). Moreover, by exploring relation (2.17) near the boundary using the expansions (3.7)-(3.13), one finds that, We shall discuss how the charge density can be related to Φ 2 in Sec. 4.

Equilibrium solutions
Here we briefly review the main features of the equilibrium solution of the 1RCBH plasma and its thermodynamics [47,82].

Thermodynamics
In equilibrium this model has the following analytical solution (written in a slightly different chart from Eq. (2.11), which we call the "modified EF" coordinates and denote with a tilde), wherer h is the radius of the black hole's event horizon given bỹ As discussed in Refs. [47,82] this model is characterized by two non-negative parameters (Q, M ) or, alternatively, (Q,r h ). They are related to the Hawking temperature of the black hole and the U (1) R-charge chemical potential In fact, standard algebraic manipulations show that (4.10) Since Q/r h is non-negative, Eq. (4.10) implies that µ/T ∈ 0, π/ √ 2 . It also follows from (4.10) that for every value of µ/T ∈ 0, π/ √ 2 there are two different values of Q/r h , which in turn parameterize two different branches of solutions. By analyzing the thermodynamics in Fig. 1, one can see the two branches merge at µ/T = π/ √ 2, which defines a critical point marking a 2nd order phase transition.  Since 1 [88], the entropy density can be determined from Bekenstein-Hawking's relation [89,90] as follows On the other hand, the U (1) R-charge density, ρ c = limr →∞ δS/δΦ ′ , may be expressed as where the stable/unstable branches correspond to lower/upper signs, as in Eq. (4.10). Furthermore, standard thermodynamic relations give us the pressure which can be shown to be p = ε/3 with the internal energy density given by ε = T s−p+µρ. This is the expected result for a conformal field theory (CFT) in 4 spacetime dimensions.

Mapping between FG (A.1) and the modified EF (4.1) coordinates
Once the equilibrium solution is known, one may compute the one-point functions of the dual QFT by using Eqs. (A.23), (A.26), and (A.38) derived in Appendix A. However, we note that these equations are written in the so-called Fefferman-Graham (FG) coordinates and, thus, we need to find a relation between these coordinates and the modified EF coordinates in order to calculate ⟨T µν ⟩, ⟨J µ ⟩, and ⟨O ϕ ⟩ in equilibrium. First we consider the diagonal form of the line element (4.1), (4.14) The task now is to find a relation betweenr and ρ (the radial coordinate of the FG chart discussed in Appendix A). This can be achieved by evaluating the following integral This integral can be solved perturbatively inr close to the boundary in order to determinẽ r(ρ) as a series expansion. Such expansion will enable us to expand the bulk fields near the boundary in the form of Eqs. (A.12)-(A.14) , which is what we need to fix the ultraviolet coefficients γ (4)µν , A (0)µ , and ϕ (0) by comparing these expansions with the equilibrium solution. With these ultraviolet coefficients fixed, we may proceed to calculate ⟨T µν ⟩, ⟨J µ ⟩ and ⟨O ϕ ⟩ using Eqs. (A.38), (A.26), and (A.23), respectively. Thus, by perturbatively evaluating the integral (4.15) close to the boundary using the equilibrium solution, we find thatr The next step is to expand the bulk fields near the boundary and substituter by ρ using Eq. (4.16). This gives the following asymptotic results, At this stage, it is useful to rewrite the energy density (ε), the equilibrium pressure (p), and the U (1) R-charge density (ρ c ) as The above results hold for the 1RCBH plasma in equilibrium.

Mapping between modified EF (4.1) and the original EF (2.11) coordinates
The purpose of this subsection is to find relations between the coefficients of the nonequilibrium solution in the near-boundary expansions (3.7)-(3.13), and the parameters of the equilibrium solution for the 1RCBH plasma (4.1). Evidently, we can only match the equilibrium solution with the asymptotically late time behavior of the non-equilibrium solution when thermodynamic equilibrium is achieved by the system. For this reason, before we derive an expression for r(r), we define the equilibrium coefficients of the nonequilibrium solution as follows, Note that the coefficients H and Φ 2 are not included in this definition because they are constant, as stated before. In order to find the relation betweenr and r, we need to solve the following equation The above expression renders an analytical relation between r andr, where 2 F 1 is the hypergeometric function and ξ is an integration constant which we choose to be in order to obtain limr →∞ r(r) −r = 0. Just as we did in Sec. 4.2 to relater with ρ, we expand the above relation in powers ofr close to the boundary Next, we substitute the above expression into Eqs. (3.7)-(3.13) and take the asymptotic equilibrium limit defined by Eq. (4.26). The result reads On the other hand, the near-boundary behavior of the analytical static solution of the 1RCBH model obtained from Eqs. (4.2)-(4.6) is Comparing Eqs. (4.31a)-(4.31d) and Eqs. (4.32a)-(4.32d), we find that Now we want to relate the non-equilibrium metric in EF coordinates (2.11) with the expression for the metric in FG coordinates (A.1) in order to obtain expressions for the one-point functions ⟨T µν ⟩, ⟨J µ ⟩, and ⟨O ϕ ⟩ far-from-equilibrium. To do so, we need to solve the integral below dr which, again, may be solved perturbatively near the boundary with the help of Eq. (3.7). The result is, Plugging r(ρ) into Eqs. (3.7)-(3.13), one obtains With the above set of equations at hand, we obtain ⟨T µν ⟩, ⟨J µ ⟩, and ⟨O ϕ ⟩ from Eqs.
Analogously to what was done in Eqs. (4.25), one may recast the physical observables of the boundary QFT as follows where ∆p(v) is the pressure anisotropy. The above expressions hold for the far-fromequilibrium 1RCBH plasma undergoing a spatially homogeneous equilibration process. Moreover, in terms of the charge density ρ c , the bulk electric field defined in Eq. (3.16) is given by Thus, we finish here the discussion of the quantities that we need to follow when solving the full nonlinear PDEs of the 1RCBH model far-from-equilibrium. Next, we discuss how to proceed with the numerics to obtain the time evolution of the 1RCBH plasma.
5 Far-from-equilibrium solutions

Field redefinitions
First, we express the system of PDEs in a form where the numerical analysis becomes as simple as possible, mapping the infinite radial domain to a finite one. This may be accomplished by redefining the radial coordinate as follows which means that the boundary now lies at u = 0. How far we go into the bulk in a numerical simulation is an issue which shall be discussed in Sec. 5.2. Moreover, since all the nontrivial dynamics of the system depends on the subleading terms of the nearboundary expansion of the bulk fields, as indicated in Eqs. (3.1)-(3.5), we define the following subtracted fields 10 The equations of motion (2.14a)-(2.14f) rewritten in terms of these subtracted fields 10 Our redefinitions are similar to the ones given in Ref. [15]. We also remark that d+(Σs) ̸ = (d+Σ)s. read with the prime now representing ∂ u . The price we paid to remove the redundant information from the bulk fields is that the equations of motion for them are now longer. Moreover, to solve these equations we need to specify the boundary conditions which, in light of Eqs.
(3.7)-(3.13) and Eqs. (5.2), are given by The time evolution equations for B s and ϕ s , which are obtained from the definitions of d + B and d + ϕ using the redefined fields (5.2), are given by We discuss in Sec. 5.4 the numerical scheme we use to evolve in time the EMD fields provided ∂ v B s and ∂ v ϕ s are known.

Radial position of the black hole event horizon
A delicate aspect of the numerical solution of the far-from-equilibrium equations of motion is how deep into the bulk one should go when integrating the radial dependence of the PDEs. Since we are working with a finite temperature setup, there is a black hole inside the bulk and one needs to cover the entire portion of the bulk between the black hole's event horizon and the boundary of the asymptotically AdS 5 spacetime, since this is the region of the bulk geometry causally connected to the boundary QFT. In the case of timedependent backgrounds, one main difficulty is that the radial position of the event horizon is unknown a priori since it depends on the time evolution of the system. If one cuts off the radial integration before reaching the horizon, then the obtained numerical solution will be inevitably inaccurate; on the other hand, if the radial integration proceeds too deep inside the horizon, the numerical simulation will probably break down due to possible singularities (caustics) associated with strong curvatures in this deep infrared region. One possible way to deal with this issue is to use the radial shift function λ(v) associated with the residual diffeomorphism invariance of the metric and regard it as an auxiliary field to fix the (apparent) horizon position at a specific value of the radial coordinate [14]. However, we do not adopt this strategy here. Instead, as mentioned above Eq. (3.7) we set λ(v) = 0 and let the horizon position fluctuate in each simulation. We then follow the reasoning of Ref. [91] in which an infrared radial cutoff corresponding to u IR ≈ 1.01 is used in the numerical simulations. Moreover, since we are interested here in analyzing the equilibration dynamics of the system in terms of the dimensionless "time combination" vT for different values of µ/T , we choose to set the equilibrium temperature to 1/π for any value of µ/T . This gives the following relation betweenr h and Q coming from Eq. (4.8), Therefore, given some value of µ/T we use Eqs. (4.10) and (5.13) to determine the corresponding values ofr h and Q. In Fig. 2 we show the plots forr h , Q, and u h as functions of µ/T . 11 In particular, Fig. 2 (c) indicates that it is fine to adopt u IR ≈ 1.01 as the inner radial cutoff for numerical integration since this prescription already works well at zero density and the horizon position in equilibrium gets closer to the boundary as the chemical potential is increased. Moreover, in Fig. 2 (d) we illustrate a typical case with the time evolution of the event horizon 12 , showing that also for transient states our inner radial cutoff is always beyond the horizon. 12 To obtain the event horizon, one just needs to solve the equation for the radial null geodesic (by setting ds 2 = 0 in Eq. (2.11) and discarding the piece which does not depend on dv), subjected to the boundary condition r(v → ∞) = r h .

Initial states
The last step required before we numerically simulate the far-from-equilibrium evolution of the 1RCBH plasma is to choose the initial data. In the context of the homogeneous equilibration of the 1RCBH model, we need to specify three inputs to start the numerical integration of the PDEs (5.3a)-(5.3f): • The initial metric anisotropy function, B s (v 0 , u); • The initial profile for the dilaton, ϕ s (v 0 , u); • The value of the equilibrium chemical potential, µ/T .
The initial metric anisotropy will be assumed to have either the form of a "pulse" with a Gaussian parameterization [65], or simply a constant profile The parameter A controls the initial amplitude of the metric anisotropy whereas the parameters σ and u 0 are related to the width and how deep into the bulk the pulse is centered. For the initial dilaton profile, we shall consider again the constant and Gaussian forms and, additionally, we will also consider the case where the initial dilaton profile is already at its equilibrium solution, i.e. 13 Note that for B(v 0 , u) = 0 (corresponding to a pure thermalization process with no initial anisotropy), the use of the last initial condition for the dilaton field automatically gives the equilibrium solution for the 1RCBH plasma. Finally, one needs to choose a value for the equilibrium chemical potential µ/T ∈ [0, π/ √ 2]. The value µ/T = 0 corresponds to the standard N = 4 SYM plasma whereas µ/T = π/ √ 2 corresponds to the critical point. Moreover, by fixing µ/T we also fix the values of Q and u h as explained above in Eq. (5.13).
Regarding the initial energy density (and the pressure, since ε = 3p), one may find it using the expression obtained in Eq. (4.13) from the thermodynamics of the 1RCBH plasma, i.e.
13 Note that these are initial profiles for the subtracted dilaton field defined in Eqs. (5.2). Therefore, the full initial dilaton profile ϕ(v0, u) = u 2 ϕs(v0, u) vanishes at the boundary u = 0 for all these initial data, as it should be.
where the temperature T is given by Eq. (4.8). Hence, far-from-equilibrium solutions with different values of µ/T will have different energy density and equilibrium pressure. Furthermore, since the 1RCBH model is a conformal theory, one cannot use only v to measure the time evolution of the fields, i.e. one needs to compose this EF time with another quantity to produce a dimensionless time measure. Thus, since the temperature T will be the same for all the numerical simulations, as explained in Sec. 5.2, the dimensionless quantity that we use to measure the time flow is vT . A different approach, which we do not explore in this manuscript, would be to fix the energy density for all the values of µ/T and let the temperature T vary among the solutions; in this case, one would use the dimensionless quantity vε 1/4 to compare the time evolutions between different values of µ/T . This is done, for instance, in Ref. [65].

Numerical techniques
In this work, in order to deal with the radial part of the system of PDEs (5.3a)-(5.3f) we developed a new numerical code that employs the so-called pseudospectral (or collocant) method [92][93][94], which is a widespread technique used to solve the equations of motion in the characteristic formulation of numerical relativity due to its accuracy and rapid convergence.
The main idea behind the pseudospectral method is to expand the numerical solution in terms of functions that form a complete basis, converge, and are easy to compute. One family of functions that accomplish these demands are the Chebyshev polynomials of the first kind, T n (x). Indeed, for well-behaved functions, the convergence of the numerical solution is exponential with respect to the number of added polynomials [92,93]. A simple example where the solution is not so well-behaved occurs when there are logarithmic terms in the near-boundary expansion of the bulk fields; this situation led the authors of Ref. [95] to favor the Runge-Kutta method to solve the radial part of the PDEs.
Assuming that the domain of interest is u ∈ [0, u IR ], where u IR is the infrared radial cutoff, we write the numerical approximation X N (u) of a function X(u) as, where {a k } is the set of spectral coefficients. In the problem that we are interested to solve here, X(u) represents the elements of the set In the numerical calculation one needs to specify the set of grid points {u i } where the discretized equations of motion take place. For this work, we adopt the Chebyshev-Gauss-Lobatto grid, which is given by 14 22) where N is the number of grid points, also known as the collocant points.
14 As discussed in Section 5.2, we used in the numerical simulations an infrared cutoff corresponding to uIR ≈ 1.01.
Another important representation of X N is given in terms of the values it assumes at the grid points, i.e.
with p 0 = p N −1 = 2 and p i = 1 otherwise, satisfying the condition, It is then evident that the set of coefficients {X i } are precisely the values of X N (u i ). Moreover, the description of the problem in terms of {X i } is completely equivalent to the description in terms of the spectral coefficients {a i }, the choice between them being just a matter of convenience. In this work, in particular, we will solve the PDEs (5.3a)-(5.3f) using the cardinal basis representation (5.23) to get {X i }.
The last element that we need to define before numerically solving the equations of motion is the pseudospectral differentiation matrix, D ij , which provides the discrete version of X ′ (u), i.e. X ′ i = D ij X j . To obtain this matrix, one needs to go through the derivation of the cardinal functions. Here we just give the final form of the differentiation matrix [92]: 15 26) A nice property of D ij is that, to obtain higher order derivatives, one just needs to exponentiate this matrix D where n is the order of the desired derivative. The discretized version of the differential equations acquire then the following generic form where Q X is a N ×N matrix and ⃗ f X is a vector with N components. Hence, we have transformed a continuum differential problem into a linear algebra problem, i.e. an inversion matrix problem since the solution of Eq. (5.28) is given by (5.29) 15 One can obtain the differentiation matrix in Mathematica using the command: where ugrid denotes the radial grid and N is the number of collocant points.
We note that when going from Eq. (5.28) to Eq. (5.29) we assumed that Q X is invertible. This is not always true. Indeed, in many cases one needs to add an initial/boundary condition to the original Eq. (5.28) to obtain an invertible matrix. At this stage it is instructive to give a simple example to see what should be the form of Q X and ⃗ f X . For instance, the following second order differential equation assumes a discretized version given by, which means that Q X = D 2 and ⃗ f X = (2, . . . , 2) T . Now the discussion regarding the solution of the radial part of the equations of motion is almost done, we only have to discuss the filtering process.
In our calculations, we have found the common "aliasing" (a.k.a. "spectral blocking") problem where high frequency modes have spurious growth and contaminate the numerical computation until it eventually breaks down. Such issue is typical in problems involving nonlinear equations [92,93]. To circumvent this problem, we need to access the spectral coefficients a i and perform a "damping" on the higher modes. To go from {X i } to {a i } we use the Matrix Multiplication Transform (MMT) as defined in Ref. [92], where denoting the Gaussian quadrature weight function and (T i , T i ) denoting the scalar product, Our damping process, in particular, is very efficient. For instance, if we take N = 40 as the number of radial grid points, setting the last three spectral coefficients to zero at each time step suffices to produce a well behaved numerical evolution. We illustrate this process schematically below 16 (5.36) 16 Another common way to do the filtering process is to use the Fast Fourier Transform (FFT). This is explained in the review [14] and it is also used, e.g., in Ref. [19].
For the time evolution of Eqs. (5.11) and (5.12) we have used the Adams-Bashforth (AB) method, which is a well-known explicit multi-step method employed to solve ordinary differential equations. In particular, the third order AB formula is given by, where the subscript n denotes the position of the variable at the nth time step and ∆v is the time step value. Furthermore, since the AB method (5.37) requires the time derivative of a few previous steps it needs to be initialized with another method. In this work we used the Euler method for the first three steps. 17 6 Equilibration dynamics: results for different initial data

Constant metric anisotropy and dilaton profiles
We begin by discussing the results for the initial data specified in Eqs. which are similar to the first condition probed in Ref. [91]. In Fig. 3 we show a 3D plot with the time evolution of B s and ϕ s for µ/T = 2. As discussed before, we choose the dimensionless quantity vT to represent the time flow in the CFT. With B s and ϕ s now fully determined, one may compute the time evolution of the pressure anisotropy ∆p (4.45b) and the scalar condensate ⟨O ϕ ⟩ (4.45d).
The results for the isotropization process describing how the pressure anisotropy ∆p goes to zero as the time evolves, taking into account the initial condition (6.1), are shown in Fig. 4. In Fig. 4 (a), where ∆p is normalized by the (equilibrium) temperature to the fourth, one notes that at the earliest times the pressure anisotropy in insensitive to the value of the dimensionless chemical potential µ/T which, however, becomes relevant as the time evolves. In fact, as we shall discuss in Sec. 6.7, the isotropization time (which cannot be adequately resolved by eyeball in the scale of the plot shown in Fig. 4 (a)) has a non-monotonic dependence on the value of µ/T , which is a direct consequence of the presence of a critical point in the phase diagram of the 1RCBH model. Namely, as it will become clear in the analysis of Sec. 6.7, the isotropization time may decrease or increase for increasing values of µ/T depending on whether µ/T is far or close enough to the critical point, respectively. This is very different from what happens e.g. in the case of the AdS-Reissner-Nordstrom background investigated in Ref. [65], which does not feature a critical point in its phase diagram. It is also interesting to note that, as anticipated around Eq. (5.20), the equilibrium pressure (and energy density) depends on the chosen value of µ/T , which is clear from the curves shown in Fig. 4 (b), where ∆p is normalized by the equilibrium pressure. Moreover, we have verified the effect on the isotropization time obtained by varying the value of the initial metric anisotropy and the effect was negligible. For instance, if one doubles the initial anisotropy, the only effect is to have a steeper downfall of ∆p before it reaches the first minimum. The explanation for this is the surprising low nonlinearity of the system [91,96], which is also quantitatively explored at length in Ref. [65]. The time evolution of the scalar condensate ⟨O ϕ ⟩ for different values of µ/T is presented in Fig. 5. The first remarkable feature of the time evolution of the scalar condensate is that it takes a much longer time to relax to equilibrium than the pressure anisotropy. This is the reason why we adopted the equilibration time associated with the relaxation of the scalar condensate toward equilibrium as the true thermalization time of the 1RCBH plasma. We also observe a qualitatively different behavior for the dependence of the thermalization time with µ/T when compared to the corresponding dependence of the isotropization time (note also that in the case of the AdS-Reisser-Nordstrom background studied in Ref. [65], spatially homogeneous isotropization corresponds already to the true thermalization of the system since there is no scalar condensate in that case). The thermalization time associated with the relaxation of the scalar condensate always increase with increasing µ/T . It is also interesting to note the qualitative similarity between the time evolution of ⟨O ϕ ⟩ and the thermalization process of confining theories studied in Ref. [95] using holographic quenches, even though the 1RCBH model is non-confining.

Constant metric anisotropy profile and Gaussian dilaton profile
In this subsection we study the time evolution of the system given the initial data: where we kept the same constant initial metric anisotropy as before but changed the initial condition for the dilaton field and considered now a Gaussian profile. This is a way to look for possible new features depending mostly on the choice for the initial dilaton profile.  In Fig. 6 we display a 3D plot with the time evolution of B s (v, u) and ϕ s (v, u) for µ/T = 2 and the initial conditions set by Eq. (6.2). By comparing with Fig. 3, one notes that the evolution of the metric anisotropy is not very sensitive to the variation of the initial profile for the dilaton field. On the other hand, as expected, the evolution of the dilaton field is significantly affected by the choice of its initial value. This seems to suggest that the backreaction produced on the metric anisotropy by varying the initial profile for the dilaton is small, which will be confirmed in the course of the next subsections.
In Fig. 7 we show our results for the time evolution of the pressure anisotropy ∆p for different values of µ/T and the initial conditions given in Eq. (6.2). By comparing with the results obtained in the previous subsection, displayed in Fig. 4, one notes that the pressure anisotropy does not change much by varying the initial condition for the dilaton, if we keep fixed the constant initial metric anisotropy. This is a direct consequence of the previously mentioned robustness of the metric anisotropy against variations of the initial dilaton profile. On the other hand, the early time dynamics of the scalar condensate ⟨O ϕ ⟩ displayed in Fig. 8 is very different from the result obtained in the previous subsection and shown in Fig. 5, while its late time dynamics is remarkably similar for the two different set of initial conditions given in Eqs. (6.2) and (6.1). Moreover, the same observations done in the previous subsection, that the thermalization time associated with the equilibration of the scalar condensate only happens significantly after the system has already relaxed to an (approximately) isotropic state, and that this thermalization time always increase with increasing µ/T , also hold for the initial conditions given in Eq. (6.2). As we are going to see in the next subsections, these qualitative trends hold for all the different initial conditions considered in the present work, suggesting that they are general features of the equilibration process of the 1RCBH plasma.

Constant metric anisotropy profile and equilibrium dilaton profile
In this subsection we study the time evolution of the system given the initial data: where we kept the same constant initial metric anisotropy as in the two previous subsections but considered now the initial dilaton profile discussed in Eq. (5.19), corresponding to the equilibrium value of the dilaton field. In Fig. 9 we display a 3D plot with the time evolution of B s (v, u) and ϕ s (v, u) for µ/T = 2 and the initial conditions set by Eq. (6.3). By comparing with the results in the two previous subsections, one confirms once again the robustness of the time evolution of the metric anisotropy against different choices for the initial dilaton profile.
In Fig. 10 we display our results for the time evolution of the pressure anisotropy, which are similar to what we have found in the two previous subsections. Regarding the time evolution of the scalar condensate shown in Fig. 11, we see that its early time dynamics is different from the previous cases considered here, although its late time dynamics is very similar to those found before. Furthermore, one also notes that in the present case where the initial condition for the dilaton is already its equilibrium value that the fluctuations in the value of the dilaton as time evolves are generally small. This is expected since in the present case the far-from-equilibrium dynamics of the system is being initially driven solely by the metric anisotropy.  Figure 11: Time evolution of the scalar condensate ⟨O ϕ ⟩ for different µ/T using the initial data (6.3).

Gaussian metric anisotropy profile and constant dilaton profile
Now we consider a different set of initial conditions to scan out possible new features in the equilibration dynamics of the 1RCBH plasma as a function of the chosen initial data. We change the initial metric anisotropy profile B s (v 0 , u) but use the same constant initial profile for the dilaton field ϕ s (v 0 , u) as in subsection 6.1, which is a way to probe if there are some new features depending mostly on the chosen initial anisotropy. We present in this subsection the results for a Gaussian initial metric anisotropy (5.15), which is perhaps the most common initial condition chosen for this field, combined with a constant initial profile for the dilaton field (5.17) where the Gaussian initial anisotropy was chosen to be more or less half-way between the black hole horizon and the boundary. In Fig. 12 we show a 3D plot with the time evolution of B s (v, u) and ϕ s (v, u) using the initial condition given in Eq. (6.4) and µ/T = 2. The shape of both functions are very different than the case considered in subsection 6.1, whose results are displayed in Fig. 3 and where the initial condition for the dilaton was the same as in Eq. (6.4) but the initial metric anisotropy was constant. This signalizes a strong backreaction induced on the dilaton field by changing the initial metric anisotropy, in contrast to the general trend observed in the previous subsections that the backreaction induced on the metric anisotropy by changing the initial dilaton profile is quite small. These general features of the far-from-equilibrium 1RCBH plasma will be further confirmed in subsections 6.5 and 6.6.
The time evolution of the pressure anisotropy ∆p that we obtain from the boundary value of B s (v, u) in the present case is shown in Fig. 13. One can see that the early time  Figure 14: Time evolution of the scalar condensate ⟨O ϕ ⟩ for different µ/T using the initial data (6.4).
dynamics of the pressure anisotropy in the case of the initial condition set in Eq. (6.4) is completely different than the one obtained in Fig. 4 by considering the initial condition given by Eq. (6.1). In particular, the isotropization time has significantly decreased by considering a Gaussian initial metric anisotropy when compared to the isotropization times associated with a constant initial metric anisotropy analyzed in the previous subsections. The time evolution of the scalar condensate ⟨O ϕ ⟩ for the initial condition (6.4) is shown in Fig. 14. Compared with the result presented in Fig. 5 using a constant initial anisotropy (6.1), we note that the early time dynamics of ⟨O ϕ ⟩ is very different depending on the chosen initial data. However, the late time dynamics of the scalar condensate and the associated thermalization time are actually very similar for both sets of initial conditions. These observations are in agreement with the general trend observed in the previous subsections and will be further confirmed in the course of the next subsections.

Gaussian metric anisotropy and dilaton profiles
Now we change the initial dilaton profile ϕ s (v 0 , u) but maintain the prior parametrization for the initial metric anisotropy B s (v 0 , u) The 3D plot with the time evolution of B s (v, u) and ϕ s (v, u) is shown in Fig. 15 for the initial condition (6.5) and µ/T = 2. From this plot, it is evident the similarity of the shape of the metric anisotropy B s (v, u) with the previous result shown in 12, i.e. a different parameterization for the initial dilaton profile ϕ s (v 0 , u) has essentially no effect  Figure 17: Time evolution of the scalar condensate ⟨O ϕ ⟩ for different µ/T using the initial data (6.5).
on the time evolution of the metric anisotropy. On the other hand, one notes a very different time evolution for the dilaton ϕ s (v, u), which is expected since we changed the initial parameterization of this field. As a direct consequence of the features above, the time evolution of the pressure anisotropy ∆p depicted in Fig. 16 is essentially the same as in the previous case shown in Fig. 13. The difference between them can be barely seen (what configures an even stronger test of robustness than observed before by considering initial conditions with a fixed constant initial metric anisotropy), indicating that the pressure anisotropy is remarkably robust against the addition of other fields in the gravitational action besides the metric. On the other hand, the early time dynamics of the scalar condensate ⟨O ϕ ⟩ presented in Fig. 17 is very different from the previous case shown in Fig. 14, although its late time dynamics and the associated thermalization time are very similar to the ones obtained in the previous subsections.

Gaussian metric anisotropy profile and equilibrium dilaton profile
Finally, the last case which remains to be analyzed correspond to the initial condition with Gaussian initial metric anisotropy and an initial profile for the dilaton field equal to its equilibrium value Figure 18: (Color online) Results for the time evolution of some fields involved in the 1RCBH setup for the initial condition (6.6) with µ/T = 2: (a) the subtracted metric anisotropy function B s (v, u), and (b) the subtracted dilaton field ϕ s (v, u).
In Fig. 18 we present the 3D plot with the time evolution of B s (v, u) and ϕ s (v, u) using the initial condition (6.6) and µ/T = 2. As before, the metric anisotropy B s (v, u) is effectively unaffected by the change done in the initial condition for the dilaton field, while the dilaton profile ϕ s (v, u) has only very small fluctuations in time, which is expected since we have already started the numerical simulations in the present case with an initial   Figure 20: Time evolution of the scalar condensate ⟨O ϕ ⟩ for different µ/T using the initial data (6.6).
dilaton profile equal to its value in thermodynamic equilibrium. The pressure anisotropy evolves in time according to Fig. 19, which displays essentially the very same behavior found in the last two subsections. On the other hand, the early time dynamics of the scalar condensate depicted in Fig. 20 is different than what we have seen in the previous subsection, displaying only very small fluctuations in time. In any case, the thermalization time associated with the equilibration of the scalar condensate is found to be very similar to what has been observed in all the previous cases considered here. Therefore, after we have exhausted the analysis of the selected sets of initial conditions, we may draw some general conclusions for the homogeneous equilibration dynamics of the 1RCBH plasma analyzed here: 1. The backreaction produced in the time evolution of the dilaton field by changing the initial metric anisotropy is extremely large, in contrast to what happens with the backreaction produced on the time evolution of the metric anisotropy by changing the initial dilaton profile, which is small. This implies that the pressure anisotropy is robust against the addition of other fields in the gravitational action besides the metric; 2. The isotropization (associated with the approximate vanishing of the pressure anisotropy ∆p) always happen (well) before the true thermalization of the system (associated with the equilibration of the scalar condensate ⟨O ϕ ⟩); 3. Regarding the late time dynamics of the scalar condensate and the associated thermalization time, these are robust features of the medium which remain almost unchanged regardless of the initial conditions considered. In particular, the thermalization time always increases with increasing chemical potential. On the other hand, the early time dynamics of the scalar condensate is strongly dependent on the set of initial data chosen to seed the time evolution of the far-from-equilibrium 1RCBH plasma; 4. The time evolution of the pressure anisotropy and the associated isotropization time of the system strongly depend on the chosen initial condition for the metric anisotropy.
However, in the different cases considered above, due to the large scales plotted to cover the large amplitudes of oscillation observed in the time evolution of the pressure anisotropy, it was not possible yet to clearly reveal the crucial role played by the critical point in the isotropization of the system. We analyze this issue in detail in the next subsection, where we confront the late time dynamics of the pressure anisotropy with the lowest non-hydrodynamic QNM of the external scalar channel of the 1RCBH plasma obtained in Ref. [82]. We also compare the late time behavior of the scalar condensate with the lowest non-hydrodynamic QNM of the dilaton channel, to be derived in Appendix B.

Matching the quasinormal modes
In this section we compare the late time behavior of the full nonlinear evolution of the equations of motion with the quasinormal modes of the system, which describe exponen-tially damped collective excitations produced in response to disturbances of a black hole background (see e.g. [97,98] for some recent reviews). The QNM spectra of the theory are responsible for the so-called "quasinormal ringdown" phenomenon describing the linear part of the decaying perturbations of a disturbed black hole at late times.
The purpose of this comparison is threefold: it can give a measure of the degree of nonlinearity of the full out-of-equilibrium solutions of the equations of motion 18 ; it may be also used to provide an independent check of the accuracy of the numerical solutions at late times; finally, it can unveil some of the main differences between the early and the late time equilibration dynamics of the system.
Part of the spectra of QNMs of the 1RCBH plasma were recently obtained in Ref. [82]. By following the classification given in Ref. [47] for the gauge and diffeomorphism invariant perturbations of the EMD fields in the homogeneous zero wavenumber limit, in which case the system has a rotational SO(3) symmetry, one has three different channels to analyze corresponding to the three lowest dimensional representations of SO (3): the tensor/quintuplet channel, the vector/triplet channel, and the scalar/singlet channel. By working with the non-normalizable modes 19 of each channel, one obtains the associated transport coefficients through the use of Kubo formulas, namely, the shear viscosity from the SO(3) quintuplet channel, the charge conductivity and diffusion from the SO(3) triplet channel, and the bulk viscosity from the SO(3) singlet channel [47]. On the other hand, by working with the normalizable modes 20 , one obtains the QNMs of each channel. In Ref. [82] we analyzed in details the QNMs of the quintuplet and triplet channels and in Appendix B we shall obtain the QNMs of the singlet channel. 21 The lowest non-hydrodynamic QNMs of the SO(3) quintuplet and singlet channels will be compared in what follows with the late time dynamics of the pressure anisotropy and the scalar condensate, respectively.
In Fig. 21 we display our numerical results for the time evolution of the pressure anisotropy using the initial conditions set in Eq. (6.1) for different values of µ/T . 22 This time we plot the logarithmic of ∆p which makes it possible to clearly resolve the late time exponential damp of the black hole oscillations. The inset displayed in Fig. 21 gives a zoom of the behavior of the pressure anisotropy in the far-from-equilibrium zone. To check whether the late time dynamics of the pressure anisotropy can be described by the lowest non-hydrodynamic QNM of the SO(3) quintuple (external scalar) channel, we parametrize 18 Note that in the analysis of QNMs, by assuming small perturbations, one just considers quadratic fluctuations of the bulk fields in the disturbed action, which gives linearized equations of motion for these perturbations. 19 These are the solutions of the linearized equations of motion for the perturbations with the conditions that these perturbations are regular at the horizon and normalized to unity at the boundary. 20 These are the solutions of the linearized equations of motions which are regular at the horizon and subjected to the Dirichlet condition that these perturbations vanish at the boundary. 21 The SO(3) quintuplet channel coincides with the perturbation for a massless external scalar field [47,49,82], while the SO(3) singlet channel may be identified more generally with a so-called "dilaton channel", as we are going to explain in Appendix B. 22 We remark that although the early time behavior of this plot changes significantly depending on the initial conditions considered, the qualitative behavior observed at late times is the same for the different initial conditions considered here. ∆p using the following functional form where Re[ω]/T , Im[ω]/T , A, and φ are fixed by fitting the functional form (6.7) to the numerical result for the pressure anisotropy evaluated within the late time interval vT ∈ [1.8, 4.1]. 23 The relevant parameters extracted from this fitting procedure are Re[ω]/T and Im[ω]/T , which may be then compared with the real and imaginary parts of the non-hydrodynamic QNM with lowest imaginary part (in magnitude), corresponding to the dominant, less damped mode in the external scalar channel. Such comparison is done in Fig. 22, from which one notes a good agreement between the late time decay of the pressure anisotropy and the lowest non-hydrodynamic QNM of the external scalar channel obtained in Ref. [82]. Such agreement provides an independent check of the accuracy of our numerical routine at late times, as aforementioned. More importantly, we are now able to understand how the isotropization time depends on the presence of a critical point in the phase diagram of the 1RCBH plasma.
We note from Fig. 21 that the early time dynamics of the pressure anisotropy is qualitatively different from its late time dynamics. At early times, the pressure anisotropy is always damped as one increases µ/T , however, at late times the pressure anisotropy may decrease or increase with increasing chemical potential depending on whether the system is far or close enough to the critical point. Therefore, the isotropization time of the 1RCBH plasma is affected by the critical point in accordance with the behavior of the "relaxation time" extracted from the imaginary part of the lowest non-hydrodynamic QNM of the external scalar channel of the theory [82]. This is explained by the fact that, as shown in Fig. 22, the late time decay of the pressure anisotropy is in fact described by the lowest non-hydrodynamic QNM of the SO(3) quintuplet channel. This is qualitatively different from the thermalization time discussed before, which is associated with the equilibration of the scalar condensate. This process only sets in after a nearly isotropic state has been already achieved by the system and it always increase with increasing µ/T . In Fig. 23 we show our numerical results for the time evolution of the difference between the scalar condensate and its equilibrium value using the initial conditions set in Eq. (6.1) for different values of µ/T . To check whether the late time dynamics of this observable can be described by the lowest non-hydrodynamic QNM of the SO(3) singlet (dilaton) channel, we parametrize it using once more the following functional form the dilaton channel is displayed in Fig. 24. One notes a good agreement between both results, which shows that the late dynamics of the scalar condensate is indeed dominated by the lowest non-hydrodynamic QNM of the dilaton channel. This also gives another independent check of the accuracy of our numerical routine at late times.
In summary, we see that the late time dynamics of the pressure anisotropy and the scalar condensate are dominated by the lowest non-hydrodynamic QNMs of the external scalar and dilaton channels, respectively. Consequently, the behavior of the isotropization and thermalization times of the 1RCBH plasma can be correctly inferred from the analysis of the QNMs of the system, even though the early time dynamics of these observables cannot be described by these linearized perturbations.

Outlook and final remarks
Now we summarize the main findings of the present work concerning the spatially homogeneous equilibration dynamics of the 1RCBH plasma. The relevant one-point functions of the 1RCBH model correspond to the expectation value of the stress-energy tensor ⟨T µν ⟩ dual to the bulk metric g µν , the scalar condensate ⟨O ϕ ⟩ dual to the bulk dilaton field ϕ, and the U (1) R-charge density ρ c = ⟨J t ⟩ associated with the bulk Maxwell field A µ . The charge density is time-independent in this spatially homogeneous setup and it only depends on the value of the chemical potential of the medium. From ⟨T µν ⟩ one extracts the pressure anisotropy ∆p which describes the isotropization of the system, while the approach of the scalar condensate toward its thermodynamic equilibrium value is associated with the true thermalization of the 1RCBH plasma. This is so because the scalar condensate was found here to always equilibrate only after the system has already reached a nearly isotropic state, being always the last equilibration time scale of the system. In other words, generally the isotropization time is shorter than the thermalization time of the 1RCBH plasma.
Moreover, while the thermalization time was found to always increase with increasing values of the U (1) R-charge chemical potential, in agreement with the behavior of the lowest non-hydrodynamic QNM of the dilaton channel, the isotropization time shows a non-monotonic dependence on the chemical potential, namely, it decreases or increases with increasing chemical potential depending on how far or close to the critical point the system is, respectively. This behavior of the isotropization time is in consonance with the behavior of the lowest non-hydrodynamic QNM of the external scalar channel of the 1RCBH plasma [82]. These two observations are explained by the fact that the late time dynamics of the pressure anisotropy and the scalar condensate are dominated by the lowest non-hydrodynamic QNMs of the external scalar and dilaton channels, respectively. These qualitative findings are robust for the 1RCBH model, in the sense that they hold for different initial conditions chosen to seed the time evolution of the far-from-equilibrium plasma.
On the other hand, the early time dynamics of the pressure anisotropy and the scalar condensate are strongly dependent on the chosen initial data. More specifically, the backreaction produced in the time evolution of the dilaton field by changing the initial metric anisotropy was found to be very large, in contrast to what happens with the backreaction produced on the time evolution of the metric anisotropy by changing the initial dilaton profile, which was generally found to be small. This suggests that the pressure anisotropy is robust against the addition of other fields in the gravitational action besides the metric, being only significantly sensitive to the initial profile chosen for the metric anisotropy. In particular, the isotropization time may be significantly modified by changing the initial metric anisotropy. On the other hand, while the early time dynamics of the scalar condensate is strongly dependent on the set of initial data chosen to seed the time evolution of the far-from-equilibrium 1RCBH plasma, its late time dynamics and the associated thermalization time of the medium were found to be remarkably similar for all the different initial conditions considered here.
Some follow-ups of the present work which will appear soon regard the generalization of the analysis of the spatially homogeneous equilibration dynamics of the 1RCBH plasma performed here to consider spatially dependent hydrodynamic flows, such as the Bjorken [99] and the Gubser [100,101] flows. The latter is particulary interesting since the solution for the hydrodynamic Navier-Stokes approximation in this case displays general inconsistencies, such as regions in the fluid with negative temperature [100], which do not appear after resummation [102,103]. We also intend to investigate the collisions of shock waves in the 1RCBH setup in the near future to study the interplay between critical phenomena and the rapid expansion developed by the system under these conditions.
As discussed at length in the introduction, the 1RCBH plasma (as any conformal plasma) is not suited for direct applications in heavy ion phenomenology. Nonetheless, some of the qualitative features disclosed in the present analysis may be general enough and hold for other strongly coupled holographic fluids. Indeed, the fact that the thermalization time of the 1RCBH plasma increases with increasing chemical potential is in consonance with the suggested delayed equilibration of heavy ion collisions at lower energies or higher densities discussed in the introduction. In order to look for possible signatures of far-fromequilibrium universal dynamics of holographic systems, we also intend to generalize the analysis of the present work to the phenomenologically realistic EMD construction of Ref. [53], which provides a quantitative description of QCD thermodynamics at zero and finite baryon chemical potential.
Another perspective for future investigations consists in going beyond the calculation of one-point functions and study also higher order correlation functions [104][105][106][107][108] in far-from-equilibrium holographic settings. In particular, this kind of study may be of relevance for many different ongoing and future low energy heavy ion experiments, such as the BES program [109] being conducted at RHIC, the future fixed target (FXT) experiments [110,111] also at RHIC, the ongoing HADES experiment [112] at GSI, the future Compressed Baryonic Matter (CBM) experiment at FAIR/GSI [113,114], and also experiments in the future NICA facility [115]. 25 This is so because, as recently discussed in Refs. [116,117], the behavior of out-of-equilibrium critical cumulants in the QCD phase diagram may be very different from the equilibrium behavior of these real-time correlation functions [118][119][120][121]. Since ratios between some of these cumulants may be compared with ratios between moments of particle multiplicity distributions measured in heavy ion colli- sions, they are the main smoking gun used in experiments in order to try to identify clear experimental signatures of the QCD critical point in heavy ion collisions. Therefore, an understanding of the behavior of higher order correlation functions in out-of-equilibrium strongly coupled systems may disclose some fundamental insights which may potentially drive the experimental searches for the QCD critical point in the near future.

A Holographic renormalization
In this Appendix we give the details on how one may obtain the one-point functions ⟨T µν ⟩, ⟨J ν ⟩, and ⟨O ϕ ⟩ using the holographic renormalization procedure [87,122,123].
In what follows, and as it is common in the treatment of holographic renormalization, we adopt the Fefferman-Graham (FG) coordinates in which there is an explicit relation between the renormalization group (RG) flow at the boundary QFT and the bulk radial coordinate with the Greek indices running through the coordinates of the dual QFT, x ∈ {t, ⃗ x}, and ρ denoting the radial coordinate in the FG chart, where the boundary lies at ρ = 0. The first step to renormalize the on-shell action is to identify, in a covariant manner, what are the divergent terms. This analysis is done by integrating out the ρ-direction in the on-shell action up to a near-boundary hypersuface ρ = ϵ that acts as a cutoff, defining then a regulated action, S reg = (S bulk + S GHY )| ϵ . Once the divergences of the regulated action are identified, the counterterm action S ct is defined as follows [87,122,123] S ct = −(divergent terms of S reg ). (A. 2) The subtracted action, S sub , which is supposed to be evaluated at the cutoff ρ = ϵ, is given by The renormalized on-shell action is obtained by taking the limit ρ → 0 on the subtracted action, i.e. S ren = lim Once the renormalized on-shell action (A.4) is found, we follow the holographic dictionary and take functional derivatives of the renormalized action with respect to the boundary values of the bulk fields to obtain the corresponding one-point functions in the dual QFT. In particular, for an EMD model, the important one-point functions are 26 δS ren δg µν where g µν = ρ γ µν is the metric of the boundary QFT, which we shall take to be Minkowski at the end of the calculations. The subscript (0) denotes that these fields are computed at the boundary of the asymptotically AdS space; we will give the precise meaning of it below when we expand the fields near the boundary.

A.1 Counterterm action
The counterterm action for the 1RCBH model, which is an EMD model whose bulk scalar field has dimension ∆ = 2, is the same counterterm action for the Coulomb branch flow [87] S ct = 1 where f (0) = f (ϕ = 0), and R[γ], R µν [γ], are the respective Ricci scalar and Ricci tensor of the induced metric at the boundary, γ µν . From now on, though, in order to simplify the notation, we will suppress the explicit metric dependence γ of the curvature tensors evaluated at the boundary of the bulk space, e.g. R ≡ R[γ]. Also, from the term multiplying ln ρ one can already see what is the trace anomaly of the theory, which is zero for the case of the conformal 1RCBH model 27 . Regarding the derivation of the counterterm action, we suggest Ref. [124] for very enlightening and clear discussions about it. Also, one may find insightful discussions about the derivation of counterterm actions in the EMD context using the Hamiltonian approach in Ref. [125]. Moreover, it is important to know that, due to the fact that in the 1RCBH model the scalar field and the Abelian gauge field do not break the original conformal symmetry of the SYM plasma, one may add finite contributions to the counterterm action (A.8), which unveils the scheme dependence of the holographic renormalization procedure. The finite counterterms that one may add are 26 Note that for a bulk scalar field with dimension ∆ = 2, as it is the case of the dilaton in the 1RCBH model, we need to introduce an extra ln ρ term to regulate the expectation value of its dual scalar operator in the boundary QFT. 27 We remark that, although a chemical potential does not induce a trace anomaly, a magnetic field does induce a trace anomaly in the SYM plasma [65].
where {c 1 , c 2 } ∈ R. In short, to see why these terms are finite, one just needs to recall that, with the scaling dimensions of the EMD fields for the 1RCBH model, one obtains Consequently, one may try to simplify the final expressions for the one-point functions of the dual QFT by including some finite counterterms, i.e.
In this work, though, we will not resort to the addition of any extra finite term to the counterterm action (A.8).

A.2 One-point functions
With the counterterm action at hand, we now have the subtracted on-shell action (A.3), which means that we can proceed with the functional derivatives to extract the one-point functions given in Eqs. (A.5)-(A.7). The analysis for the scalar field carried out here is based on Appendix C of Ref. [124], whilst the vector field analysis is based on Ref. [126]. It is clear from Eqs. (A.5)-(A.7) that we need to expand the fields near the boundary. Thus, we perform the FG expansion of the EMD fields 28 Note that in order to obtain the one-point functions, it suffices to expand the fields up to O(ρ) since the remaining terms vanish as ρ → 0. The independent terms of the above expansions are {ϕ (0) , ϕ (0,1) , A (0)µ , A (2)µ , γ (0)µν , γ (4)µν } and, thus, any other coefficient may be recast in terms of these independent ones. Furthermore, we are keeping here the analysis fairly general for any EMD model with ∆ = 2; we shall only specialize to the 1RCBH background at the end of the calculations. Next, one substitutes the above near-boundary expansions for γ µν , A µ , and ϕ into Eqs. (A.5)-(A.7) to obtain the explicit formulas for the one-point functions. However, since the required algebra is not so simple, we give some further details below. First, let us provide the formulas for the variation of the regularized on-shell action with respect to 28 Note that in the subscripts (n, m) of the coefficients of these expansions, n = 0 denotes the leading order term in ρ (as it goes to zero) and m denotes the power of ln ρ. This is reminiscent of the general form of these expansions presented in Eqs. the sources 29 with K µν = ρ∂ ρ γ µν denoting the extrinsic curvature of the boundary and K = γ µν K µν its trace. We also have defined the following objects Moreover, the leading order solution of Maxwell's equations (2.9) under the expansion (A.12)-(A.14) give us a relation between A (2,1)µ and A (0)µ , i.e.
When we integrated the ρ coordinate by parts in the variation of the integrals we considered that the normal vector n M = (−2ρ, 0, 0, 0, 0) at the boundary of the manifold is "outward-pointing", i.e. gMN n M n M = 1.
which leads us to the final form the U (1) R-current, Notice, however, that the last term of the above equation is absent in the 1RCBH background since ∂ µ A (0)ν = 0. Regarding the one-point function ⟨T µν ⟩, the algebra is a little bit more complicated. Thus, in order to simplify the analysis, we will only focus on the finite contributions of Eq. (A.15) 30 . The finite terms coming from T reg µν are On the other hand, T ct µν also contributes to the finite part of the total stress-energy tensor, i.e.
Hence, the counterterms do impact on the finite result for the one-point functions, even though their original purpose was to eliminate the divergences of the on-shell action.
Proceeding with the tensorial algebra to obtain the one-point function of the stressenergy tensor, the next step is to sum Eq. (A.27) with Eq. (A.28), i.e.
Now we are close to give the full expression for ⟨T µν ⟩. The final step consists in expressing the coefficients of the metric expansion, such as γ (2) , in terms of the curvature tensors of the boundary metric γ (0)µν . This step, though, requires more laborious algebra, and more definitions are needed in order to do it in an simple way.
A convenient way to expand the EMD equations of motion using the FG coordinates is to use the ADM decomposition of general relativity [77,78] in which one considers spacelike foliations keeping ρ = constant at each foliation. We suggest at this point Ref. [125] for a more complete discussion about this subject. Using the ADM decomposition, 30 Since all the divergences are mutually canceled out by taking into account the counterterms. from where the Gauss-Codazzi equations are derived, the (µν)-components of Einstein's equations become One can now expand the above equation near the boundary using Eqs. (A.12)-(A.14), obtaining which is valid for any five dimensional EMD theory with a bulk scalar field with dimension ∆ = 2 and for any kind of boundary. Specializing the above results for the 1RCBH background will vastly simplify the expression for ⟨T µν ⟩. First, the conformal boundary of such theory is flat (i.e., Minkowski), which means that all the curvature tensors are identically zero. Second, if one takes the trace of Eq. (A.36), the resulting expression reads κ 2 5 ⟨T µ µ ⟩ = −2ϕ (0) ϕ (0,1) + ϕ 2 (0,1) − Hence, one arrives at a very important result: in a conformal EMD model with ∆ = 2 the logarithmic terms on the near-boundary expansions of the bulk fields are absent. This sentence justifies the prior assumption made in Sec. 3 when we set to zero the logarithmic terms of the near-boundary asymptotic expansions. Furthermore, we remark that the finite counterterm contribution (A.9) does not modify the trace anomaly of the theory. From the previous discussion, the stress-energy tensor for the 1RCBH model reads, which is the one-point function adopted throughout this paper. We remark once more that, due to the fact ∂ µ A (0)ν = 0 in our setup, the Maxwell terms in Eq. (A.38) will vanish. A minimal internal consistency check of the above results may be done by looking at the trace Ward Identity [125,127], i.e. where A denotes the trace anomaly of the theory. For the specific case of the EMD theory with ∆ = 2, the anomaly is given by where, (A.41)

B Quasinormal modes for the SO(3) singlet (dilaton) channel
In this Appendix we obtain the QNMs of the 1RCBH model for the SO(3) singlet (dilaton) channel in the homogeneous, zero wavenumber limit complementing the study carried out in Ref. [82] where the SO(3) quintuplet (external scalar) and triplet (vector diffusion) channels have been analyzed in detail.
As discussed in Ref. [47], at zero wavenumber the EMD system features a rotational SO(3) symmetry under which the gauge and diffeomorphism invariant perturbations of the system are organized into a quintuplet channel corresponding to the five traceless spatial components of the perturbation of the metric field h ij , a triplet channel corresponding to the three spatial components of the perturbation of the Maxwell field a i , and a singlet channel corresponding to a combination involving the dilaton perturbation φ and the trace of the spatial part of the perturbation of the metric field, namely, where one sees that the background dilaton field ϕ couples the dilaton fluctuation with the spatial trace of the graviton. This S-perturbation is analogous to the Z 2 -mode of the socalled "non-conformal channel" discussed in Ref. [128]. 31 This was called a non-conformal channel because this mode is intrinsically associated with the background dilaton field, which in the Einstein-dilaton models analyzed in Ref. [128] was responsible for breaking conformal symmetry. However, the dilaton field may also preserve conformal symmetry when in the presence of other fields, as in the case of the 1RCBH model studied here. Therefore, more generally, one could say that this S-perturbation defines the "dilaton channel". This nomenclature is also adequate due to the fact that this mode shares the same near-boundary asymptotics of the background dilaton field [47,128]. The linearized equation of motion for the S-perturbation derived in Ref. [47] translates as follows to the modified EF coordinates (4.1), 32 where ω is the frequency of the plane-wave Ansatz for the S-perturbation, which shall give the QNMs of the dilaton channel under appropriate boundary conditions to be discussed below.
In order to solve the eigenvalue problem to be derived next for the quasinormal frequencies of the system we make use of the pseudospectral method as done in Ref. [82]. For this, we begin by mapping the radial coordinater ∈ [r h , ∞) into a new radial coordinatẽ r h /r =:ũ ∈ [0, 1]. By doing so and substituting the equilibrium 1RCBH backgrounds (4.1) -(4.6) into the equation of motion (B.2), one is left with a linear differential equation for the perturbation S(ũ) on top of the equilibrium 1RCBH backgrounds. In general, the near-boundary, ultraviolet asymptotic behavior of the S-perturbation is given by, S(ũ) ∼ũ 4−∆ G(ũ) +ũ ∆ Y (ũ) + · · · , asũ → 0. Since in the 1RCBH model ∆ = 2, we have two degenerate exponents equal to two. As discussed in detail in Ref. [82], the correct eigenvalue problem for the QNMs is obtained by working with the subleading, normalizable mode of the relevant perturbation (in the dilaton channel considered here, the normalizable mode is associated with the expectation value of the boundary operator dual to the dilaton field, ⟨O ϕ ⟩), which in the present case corresponds to set the Dirichlet boundary condition G(0) = 0 with Y (0) ̸ = 0. Then, by substituting S(ũ) =:ũ 2 Y (ũ) into the equation of motion for S(ũ) and defining the dimensionless quasinormal frequencyω ≡ ω/T , one obtains 32 In order to go from the domain-wall coordinates used in Ref. [47] to the modified EF coordinates (4.1), one simply uses that d/dr = (∂v/∂r)∂v + ∂r and d/dt = (∂v/∂t)∂v, where ∂v/∂r = −grr/gtt = e B−A /h and ∂v/∂t = 1.
the following differential equation for Y (ũ), The numerical routine used to solve the GEP equation (B.3) for the QNMsω(µ/T ) employed the pseudospectral method and it follows the same general steps discussed in 33 Note that in the EF coordinates the infalling wave condition at the horizon is imposed by just requiring regularity of the solutions there. This is one of the main reasons why the EF coordinates are very convenient to deal with the calculation of QNMs. Another reason is that in the domain-wall coordinates one would obtain instead of the GEP (B.3) a quadratic eigenvalue problem (QEP) forω, which is far more demanding in terms of computational costs. Note also what the QNMsω are only functions of the dimensionless combination µ/T controlled by the background parametersr h and Q (here we restrict our analysis to the QNMs evaluated on the stable branch of black hole solutions). Ref. [82] and we refer the interested reader to consult it for technical details. In Figs. 25 and 26 we show the behavior of the QNMs of the dilaton channel as a function of µ/T . In Fig. 27 we show how the characteristic "equilibration time" associated with the inverse of minus the imaginary part of the lowest non-hydrodynamic QNM of the dilaton channel behaves as a function of µ/T . Remarkably, one notes that this characteristic equilibration time is qualitatively different from the equilibration times obtained in Ref. [82] for the external scalar and vector diffusion channels, since the latter are reduced as one increases µ/T far from the CP while increasing close to the CP. In the dilaton channel, however, the equilibration time always increases with increasing µ/T , in agreement with the late time behavior of the scalar condensate as we have discussed in Subsection 6.7.