Holographic entropy production in a Bjorken expanding
hot and dense strongly coupled quantum fluid

Romulo Rougemont [email protected] Instituto de Física, Universidade Federal de Goiás, Av. Esperança - Campus Samambaia, CEP 74690-900, Goiânia, Goiás, Brazil Willians Barreto [email protected] Centro de Ciências Naturais e Humanas, Universidade Federal do ABC, Av. dos Estados 5001, 09210-580 Santo André, São Paulo, Brazil Centro de Física Fundamental, Universidad de Los Andes, Mérida 5101, Venezuela

Abstract

We analyze the time evolution of several physical observables, namely the pressure anisotropy, the scalar condensate, the charge density, and also, for the first time, the non-equilibrium entropy for a Bjorken expanding strongly coupled $\mathcal{N}=4$ Supersymmetric Yang-Mills plasma charged under an Abelian $U(1)$ subgroup of the global $SU(4)$ R-symmetry. This represents a far-from-equilibrium, hot and dense strongly coupled quantum fluid with a critical point in its phase diagram. For some sets of initial data preserving all the energy conditions, dynamically driven transient violations of the dominant and the weak energy conditions are observed when the plasma is still far from the hydrodynamic regime. The energy conditions violations get stronger at larger values of the chemical potential to temperature ratio, $\mu/T$ , indicating that those violations become more relevant as the strongly coupled quantum fluid approaches its critical regime. For some of those energy conditions violations, it is observed a clear correlation with different plateau structures formed in the far from equilibrium entropy, indicating the presence of transient, early time windows where the Bjorken expanding plasma has zero entropy production even while being far from equilibrium. The hydrodynamization of the pressure anisotropy and also the much later thermalization of the scalar condensate are generally found to be delayed, within small relative tolerances, as $\mu/T$ is increased towards criticality. The value of $\mu/T$ in the medium is enhanced by increasing its initial charge density, and/or also by reducing its initial energy density.

I Introduction

The analysis of the evolution of the medium produced in ultrarelativistic heavy ion collisions [1, 2, 3, 4, 5], which leads to the formation of tiny droplets of a strongly coupled quark-gluon plasma [6, 7, 8, 9] in high energy particle colliders, constitutes one of the main contemporary motivations to study the real time physics of far from equilibrium media, and also how hydrodynamization and thermalization are achieved at later times [10, 11].

Although QCD is the fundamental microscopic theory describing the strong nuclear interaction, a self-consistent, first principles microscopic description of heavy ion collisions is currently unavailable. This is mainly due to the fact that during the evolution of the medium created in such collisions, QCD matter passes through several different regimes which include strongly coupled physics [9], requiring therefore a non-perturbative treatment of QCD. The main of such approaches is lattice QCD, which is greatly successful in the description of QCD spectroscopy [12, 13] and QCD thermodynamics at zero [14, 15, 16, 17] and moderate values of baryon chemical potential [18, 19, 20, 21]. However, mainly due to the so-called sign problem [22, 23, 24, 25], current lattice QCD calculations face severe difficulties regarding the calculation of physical observables at large values of baryon density, which are relevant for the search of a putative critical point in the QCD phase diagram in the plane of temperature and baryon chemical potential, and also concerning the calculation of physical observables evaluated at real time, which are important in the analysis of non-equilibrium phenomena.

In face of the aforementioned limitations, several different alternative and simpler approaches based on phenomenological models, effective theories, and also toy models are usually employed to make either quantitative predictions which may be tested against experimental heavy ion data, or to obtain possible qualitative insights on some general aspects of the quark-gluon plasma.

One such alternative or complementary approach is the holographic gauge/gravity duality [26, 27, 28, 29]. Indeed, different holographic models can be devised, with either a more phenomenological and quantitative approach to QCD [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53], or with a more qualitative focus on possible broad and general properties of strongly coupled media [54, 55] (which may not be necessarily tied to specific quantitative results in QCD). In this very regard, probably the most important and general outcome of the holographic gauge/gravity duality is the almost universal result for the shear viscosity to entropy density ratio in strongly coupled quantum fluids, $\eta/s=1/4\pi$ , which is valid for any isotropic and translationally invariant gauge/gravity dual with two derivatives of the metric field in the gravity action [56, 57, 58].

Realistically emulating all the stages of heavy ion collisions in holography is something currently unfeasible, mainly because the holographic system in the gauge/gravity approach is always strongly coupled, therefore missing the correct time evolution of the effective coupling of the medium in actual heavy ion collisions. However, an interesting and important line of investigation can be established on the grounds of understanding the physical possibilities, and also putative general properties, arising from the real time dynamics of far from equilibrium strongly coupled quantum fluids, in contrast to the scenarios realized in classical weakly coupled approaches like kinetic theory [59, 60, 61]. Contrasting such different approaches may lead e.g. to a better qualitative understanding on the essential physical properties behind the several different stages comprised in actual heavy ion collisions. Particularly, in the context of far from equilibrium holography, different physical observables have been studied under several different kinds of dynamic gauge/gravity models, see e.g. [62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95] for a non-exhaustive list.

In the present work, we proceed with our analyses of far from equilibrium holographic dynamics in the context of a strongly coupled conformal plasma expanding according to the inhomogeneous boost invariant Bjorken flow [96]. We investigate how several physical observables of a top-down holographic Einstein-Maxwell-Dilaton toy model for a hot and dense plasma with a critical point in its phase diagram, called ‘One R-Charged Black Hole’ (1RCBH) model [97, 98, 99, 100, 101, 102] (see also [35]), are affected by different variations of the far from equilibrium initial data (like the initial charge and energy densities), leading to different values of the chemical potential to temperature ratio of the medium, $\mu/T$ . In a previous work [85], the Bjorken flow for this model was analyzed mainly in terms of how the hydrodynamization time associated to the convergence of the pressure anisotropy to the Navier-Stokes result is affected by increasing the value of $\mu/T$ towards criticality. By considering variations of the initial charge density of the medium, while keeping fixed both the initial profile for the metric anisotropy and the initial energy density, it was observed that, on average, the fluid takes more time to converge to the Navier-Stokes hydrodynamic constitutive relation as $\mu/T$ increases, with the hydrodynamization time of the pressure anisotropy displaying a particularly sharp rise close to the critical point of the model.

On the other hand, there are also several important aspects of the 1RCBH model undergoing Bjorken flow which have not been analyzed in [85], and whose study constitutes the specific focus of the present work. Namely: a) in order to see the late time effective thermalization of the scalar condensate, which is associated to its convergence to the corresponding thermodynamic equilibrium result, in the present work we carry out our numerical simulations for much longer times (what clearly demands a numerical code implemented using efficient programming languages for numerical purposes, as we will discuss in appendix C); b) in order to analyze several qualitatively different physical possibilities for the time evolution of the system, we deploy a more general and qualitatively varied set of initial data, now including also the analysis for variations of the initial energy density of the medium and variations of the initial profile for the bulk dilaton field (besides variations of the initial charge density and of the initial profile for the metric anisotropy); c) more importantly, for the first time, we also analyze the holographic entropy production [63, 103] in the 1RCBH model expanding according to the Bjorken flow dynamics.

Regarding the holographic entropy of a fluid in thermodynamic equilibrium, it is well known that it is dual to the Bekenstein-Hawking black hole entropy [104, 105] measured by the area of the event horizon of a static black hole in equilibrium within the higher dimensional bulk.¹¹1For other notions of the concept of entropy also employed in the holographic gauge/gravity duality in different contexts, see e.g. [106, 107, 108, 109, 110, 111, 112, 113]. For a recent non-holographic calculation of entropy production in Bjorken flow within the context of the QCD phase diagram, see [114]. However, when the fluid is far from equilibrium, there are clear indications that the area of a dynamic event horizon is no longer an adequate measure of the entropy of the medium. In fact, as argued e.g. in [115], one expects that entropy production is local in time, therefore, associating the far from equilibrium entropy to the area of a dynamic event horizon seems unnatural, since in such a context the event horizon can only be determined by knowing the entire future evolution of the black hole geometry and, consequently, it is a global rather than a local observable. Moreover, in the case of the holographic conformal soliton flow [116], which corresponds to an ideal fluid and, thus, has zero entropy production at all times, the area of the event horizon diverges [115], providing an explicit example that it is not an adequate measure of holographic entropy production. On the other hand, the area of the apparent horizon is calculated locally in time, as expected for an entropy function. Furthermore, for the conformal soliton flow, the area of the apparent horizon is constant in time [115], as expected for the entropy of an ideal fluid.²²2As discussed in [115], for the conformal soliton flow the system does not settle down to a stationary state at late times, and the apparent horizon does not converge to the event horizon in this case, contrary to what happen in cases where dissipation drives the evolution of the system. Indeed, in many other works in the holographic literature the non-equilibrium entropy has been associated to the area of a dynamic apparent horizon [63, 65, 66, 71, 70, 79, 117, 118, 119], and we also follow the same approach. Since here the apparent horizon later converges to the event horizon close to equilibrium, the holographic non-equilibrium entropy so defined is assured to converge to the usual Bekenstein-Hawking thermodynamic entropy in equilibrium, what should be observed at late times in the evolution of the dynamic system.

In two previous works [90, 91], by analyzing the dynamical evolution of some initial data, we pointed out some notable correlations between far from equilibrium plateau structures observed in the entropy, indicating the presence of transient, early time windows with zero entropy production in the Bjorken expanding $\mathcal{N}=4$ Supersymmetric Yang-Mills (SYM) plasma at zero density, and some later transient violations of classical energy conditions [120, 121] in the medium. Although violations of classical energy conditions by quantum effects are well known (see e.g. the discussions in [122, 123]), the physical possibility that such violations maybe can also happen during the dynamical evolution of the QCD medium produced in heavy ion collisions is something not currently contemplated in the literature, since phenomenological models of pre-hydrodynamic stage of heavy ion collisions are commonly modeled using classical kinetic theory approaches [124], where such violations do not take place due to the positiveness of the single particle distribution function [125].³³3It would be important to check, within transport approaches, whether such violations can be observed outside the classical regime due to quantum effects, what requires moving from classical kinetic theory to the quantum Kadanoff-Baym equations [126]. More generally, as shown in [90, 91] in the context of the holographic Bjorken flow, and as also known from some other far from equilibrium holographic dynamics [127], it is important to recognize that strongly coupled quantum fluids may generally display dynamically driven violations of classical energy conditions, even for some initial data satisfying all such conditions.

In the present work, we shall see that such violations are generally enhanced by increasing the value of $\mu/T$ in the hot and dense medium described by the holographic 1RCBH model. We will also discuss the correlations observed when plateau structures are produced in the non-equilibrium entropy and later violations of energy conditions, and how the value of $\mu/T$ affects such correlations.

This work is organized as follows. In section II, we review the holographic formulation of the Bjorken flow for the 1RCBH model, providing more details than in the previous work [85], what will be useful for the reader interested in reproducing our results, and we also briefly present some specific results for the thermodynamics of the model which will be used in the analysis of the effective thermalization of the scalar condensate and in some analytical consistency checks for the late time numerics of our calculations, involving also the non-equilibrium entropy. In section III, we discuss the form of the set of initial data which will be analyzed in the present work and how we shall organize their different variations, besides also defining the set of normalized observables we will consider to study different physical aspects of the dynamical evolution of the selected initial data; we also briefly review the formulas for the holographic non-equilibrium entropy and the dominant and weak energy conditions. In sections IV, V, and VI we analyze, respectively, the effects associated to variations in the initial charge density, the initial energy density, and the initial dilaton profile, while keeping the remaining initial data fixed. We shall work with a subset of initial metric anisotropies originally considered for the SYM plasma in [90, 91], which will imply several qualitatively different physical possibilities for the dynamic evolution of the system, allowing us to extract possible general conclusions from the analysis presented. In section VII we present a summary with our main conclusions and future perspectives, while the appendices are devoted to discussions on numerical error analysis, further physical consistency checks of our numerical solutions, and also an analysis of our numerical code’s performance, clearly showing the need for using efficient programming languages for numerical purposes in the present work.

We use a mostly plus metric signature and natural units where $c=\hbar=k_{B}=1$ .

II Holographic Bjorken flow for the 1RCBH plasma

The 1RCBH model [97, 98, 99, 100, 101, 102] (see also [35, 128, 83]) is a solution of the five dimensional $\mathcal{N}=8$ gauged supergravity action [98] which, in turn, corresponds to a consistent truncation of type IIB superstrings on AdS ${}_{5}\times S^{5}$ . This was shown to lie within a class of solutions equivalent to near-extremal spinning D3-branes on AdS ${}_{5}\times S^{5}$ [101]. The Kaluza-Klein compactification of $S^{5}$ on the spinning D3-branes solutions leads to a global $SU(4)$ R-charge symmetry, which has three independent Cartan subgroups $U(1)_{a}\times U(1)_{b}\times U(1)_{c}$ associated with three distinct conserved charges $(Q_{a},Q_{b},Q_{c})$ of the black hole background [98]. The general solution is known as the STU model, while the particular case of the 1RCBH model is obtained by considering only one charge by setting $Q\equiv Q_{a}$ and $Q_{b}=Q_{c}=0$ . The effective five dimensional theory so obtained is holographycally dual to a $\mathcal{N}=4$ SYM plasma at finite temperature with a chemical potential associated to the conserved R-charge transforming under the Abelian $U(1)_{a}$ subgroup of the global $SU(4)$ R-charge symmetry. The effective five dimensional bulk action comprises the metric field, a scalar dilaton field, and an Abelian Maxwell field,

\displaystyle S=\frac{1}{2\kappa_{5}^{2}}\int d^{5}x\sqrt{-g}\left[R-\frac{f(\phi)}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}(\partial_{\mu}\phi)^{2}-V(\phi)\right],

(1)

where $\kappa_{5}^{2}\equiv 8\pi G_{5}$ , with $G_{5}$ being the five dimensional Newton’s constant, and the dilaton potential $V(\phi)$ and the coupling function $f(\phi)$ between the dilaton and Maxwell fields are given by,

\displaystyle V(\phi)=-\frac{1}{L^{2}}\left(8e^{\phi/\sqrt{6}}+4e^{-\sqrt{2/3}\,\phi}\right),\qquad f(\phi)=e^{-2\sqrt{2/3}\,\phi},

(2)

with $L$ being the asymptotic AdS₅ radius (which from now on we set to unity for simplicity). The bulk action (1) is supplemented by two boundary terms: the Gibbons-Hawking-York action [129, 130], needed for the well-posedness of the boundary Dirichlet problem; and the counterterm action [83], required in the holographic renormalization procedure [131, 132, 133] in order to consistently remove the boundary divergences of the on-shell action.

Several different aspects of the 1RCBH model have been already studied in the literature. For instance, the thermodynamics of the model was detailed discussed in [35, 128], where it was shown that this model has a critical point at $\mu/T=\pi/\sqrt{2}$ , which is identified by the divergences of the second (and higher order) derivatives of the pressure, like the heat capacity and the charge susceptibility. Since the model has vanishing trace anomaly, it is conformal and all the dimensionless combinations of physical observables in equilibrium are functions of the dimensionless ratio $\mu/T$ (instead of functions of independent $T$ and $\mu$ , since in the conformal case there is no dimensionful scale in the vacuum of the theory, contrary to what happens e.g. in QCD, where $\Lambda_{\textrm{QCD}}$ sets the scale with respect to which other dimensionful quantities, like $T$ and $\mu$ , can be measured). For each possible value of $\mu/T$ there are two competing branches of solutions, one corresponding to thermodynamically stable black hole backgrounds, and another one corresponding to unstable black holes. The stable branch of black hole backgrounds constitutes the set of physically relevant solutions of the 1RCBH model at finite $\mu/T$ . Due to its conformal nature, the phase diagram of the 1RCBH model is a line in the $\mu/T$ axis (instead of a plane with independent $T$ and $\mu$ directions), which starts at $\mu/T=0$ , corresponding to the AdS₅-Schwarzschild solution, dual to the SYM plasma at finite temperature and zero chemical potential, and ends at the critical point $\mu/T=\pi/\sqrt{2}$ — due to this peculiar feature, although there is a critical point in the phase diagram of the 1RCBH model (where several observables diverge with characteristic critical exponents [35, 128]), there is no phase transition, since in this model $\mu/T$ can not be increased past its critical value.

Also the conductivity and charge diffusion in the 1RCBH model were analyzed in [35]. The spectra of quasinormal modes and their critical behavior were discussed in [128, 83]. The holographic renormalization and the far from equilibrium homogeneous isotropization dynamics were presented in [83]. The evolution of holographic complexity [134], the holographic entanglement entropy and mutual information [135], the entanglement of purification [136], and the critical behavior of the hydrodynamic derivative series [137] were also analyzed for the 1RCBH model. The far from equilibrium Bjorken flow dynamics and the analysis of the hydrodynamization times of the medium, including the vicinity of the critical point, were presented in [85].

Having discussed the above introductory remarks, the main purposes of the present section of the paper are the following: first, we briefly present below some specific results for the thermodynamics of the model which will be used in the analysis of the effective thermalization of the scalar condensate and in some analytical consistency checks for the late time numerics of our calculations, involving also the non-equilibrium entropy. Next, as the main focus of the present section, we will review the holographic formulation of the Bjorken flow dynamics of the 1RCBH model in more detail than in the previous short paper of Ref. [85].

Our new results will be mainly presented in sections IV, V, and VI, after a discussion on the organization of the initial data and the normalization of the relevant physical observables in section III.

II.1 Some important thermodynamic results

Refer to caption — Figure 1: Analytical equilibrium results for dimensionless ratios involving (a) the scalar condensate, (b) the entropy density, and (c) the charge density, all of them normalized by the energy density. We display the thermodynamically stable and unstable branches of equilibrium 1RCBH solutions. In the stable branch, at $\mu/T=0$ , one recovers the pure thermal SYM plasma results.

Let $X$ be any physical observable and let us define $\hat{X}\equiv\kappa_{5}^{2}X=4\pi^{2}X/N_{c}^{2}$ . By considering Eqs. (4.21) and (4.24) of [83], one obtains for the dual quantum field theory in thermal equilibrium the following result for the dimensionless ratio of the scalar condensate $\langle\hat{O}_{\phi}\rangle$ normalized by the square root of the energy density $\hat{\epsilon}$ of the medium,

\displaystyle\frac{\langle\hat{O}_{\phi}\rangle}{\hat{\epsilon}^{1/2}}\biggr{|}_{\textrm{(eq.)}}=\kappa_{5}\,\frac{\langle O_{\phi}\rangle}{\epsilon^{1/2}}\biggr{|}_{\textrm{(eq.)}}=\frac{2Q^{2}}{3M},

(3)

where $Q$ and $M$ are the charge and the mass of the equilibrium black hole solution. On the other hand, from Eq. (2.9) of [128], one gets,

\displaystyle M=\tilde{r}_{\textrm{EH}}\sqrt{\tilde{r}_{\textrm{EH}}^{2}+Q^{2}},

(4)

where $\tilde{r}_{\textrm{EH}}$ is the radial location of the event horizon in equilibrium. By substituting Eq. (4) into Eq. (3), it follows that,

\displaystyle\frac{\langle\hat{O}_{\phi}\rangle}{\hat{\epsilon}^{1/2}}\biggr{|}_{\textrm{(eq.)}}=\frac{2\left(Q/\tilde{r}_{\textrm{EH}}\right)^{2}}{3\,\sqrt{1+\left(Q/\tilde{r}_{\textrm{EH}}\right)^{2}}}.

(5)

By using Eq. (2.12) of [128], the result above can be written as below,

\displaystyle\frac{\langle\hat{O}_{\phi}\rangle}{\hat{\epsilon}^{1/2}}\biggr{|}_{\textrm{(eq.)}}=\frac{4\left[\displaystyle{\frac{1\pm\sqrt{1-\left(x/x_{c}\right)^{2}}}{x/x_{c}}}\right]^{2}}{3\,\sqrt{1+2\left[\displaystyle{\frac{1\pm\sqrt{1-\left(x/x_{c}\right)^{2}}}{x/x_{c}}}\right]^{2}}},

(6)

where we defined the notation $x/x_{c}\equiv(\mu/T)/(\pi/\sqrt{2})$ as in [85],⁴⁴4Note that the notation $x\equiv\mu/T$ has nothing to do with the spatial coordinate $x$ (which is not actually relevant in the present work). In particular, $x_{c}\equiv\pi/\sqrt{2}$ denotes the critical point. and the lower (upper) signs specify the thermodynamically stable (unstable) branch of equilibrium 1RCBH solutions. The corresponding results are displayed in Fig. 1 (a).

Also, by using Eqs. (2.15) and (2.17) of [128], and by taking into account the fact that in a four dimensional conformal theory, as in the case of the 1RCBH plasma, $\epsilon=3p$ , one obtains that,

	$\displaystyle\frac{\hat{s}^{4/3}}{\hat{\epsilon}}\biggr{\|}_{\textrm{(eq.)}}=\kappa_{5}^{2/3}\,\,\frac{s^{4/3}}{\epsilon}\biggr{\|}_{\textrm{(eq.)}}$	$\displaystyle=(4\pi^{2})^{1/3}\,\,\frac{(s/N_{c}^{2}T^{3})^{4/3}}{3p/N_{c}^{2}T^{4}}\biggr{\|}_{\textrm{(eq.)}}$
		$\displaystyle=(4\pi^{2})^{1/3}\,\,\frac{\left[\displaystyle{\frac{\pi^{2}}{16}}\left(3\pm\sqrt{1-\left(x/x_{c}\right)^{2}}\right)^{2}\left(1\mp\sqrt{1-\left(x/x_{c}\right)^{2}}\right)\right]^{4/3}}{\displaystyle{\frac{3\pi^{2}}{128}}\left(3\pm\sqrt{1-\left(x/x_{c}\right)^{2}}\right)^{3}\left(1\mp\sqrt{1-\left(x/x_{c}\right)^{2}}\right)},$		(7)

where, as before, the lower (upper) signs specify the thermodynamically stable (unstable) branch of equilibrium 1RCBH solutions. The corresponding results are shown in Fig. 1 (b). Moreover, in a completely analogous way, one can also easily obtain from Eqs. (2.16) and (2.17) of [128] the following results for the dimensionless ratio,

	$\displaystyle\frac{\hat{\rho}^{4/3}}{\hat{\epsilon}}\biggr{\|}_{\textrm{(eq.)}}=\kappa_{5}^{2/3}\,\,\frac{\rho^{4/3}}{\epsilon}\biggr{\|}_{\textrm{(eq.)}}$	$\displaystyle=(4\pi^{2})^{1/3}\,\,\frac{(\rho/N_{c}^{2}T^{3})^{4/3}}{3p/N_{c}^{2}T^{4}}\biggr{\|}_{\textrm{(eq.)}}$
		$\displaystyle=(4\pi^{2})^{1/3}\,\,\frac{\left[\displaystyle{\frac{x}{16}}\left(3\pm\sqrt{1-\left(x/x_{c}\right)^{2}}\right)^{2}\right]^{4/3}}{\displaystyle{\frac{3\pi^{2}}{128}}\left(3\pm\sqrt{1-\left(x/x_{c}\right)^{2}}\right)^{3}\left(1\mp\sqrt{1-\left(x/x_{c}\right)^{2}}\right)},$		(8)

which are displayed in Fig. 1 (c).

The thermodynamically stable equilibrium results derived above will be used to obtain the value of $\mu/T$ in the medium from the late time evolution of each far from equilibrium initial data and also to analyze the late time effective thermalization of the scalar condensate in the Bjorken flow of the 1RCBH model. Since these are analytical results, they will be also employed to make important physical consistency checks of our numerical results in the late time evolution of the system.

II.2 Bjorken flow dynamics

Now we review in more details the holographic formulation of the Bjorken flow for the 1RCBH model [85].

Let us begin by writing down the general Einstein-Maxwell-Dilaton equations of motion coming from the extremization of the bulk action (1),


$\displaystyle R_{\mu\nu}-\frac{g_{\mu\nu}}{3}\left(V(\phi)-\frac{f(\phi)}{4}F_{\alpha\beta}^{2}\right)-\frac{1}{2}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{f(\phi)}{2}F_{\mu\rho}F_{\nu}^{\rho}$	$\displaystyle=0,$	(9a)
$\displaystyle\nabla_{\mu}(f(\phi)F^{\mu\nu})=\frac{1}{\sqrt{-g}}\partial_{\mu}(f(\phi)F^{\mu\nu})$	$\displaystyle=0,$	(9b)
$\displaystyle\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\phi)-\partial_{\phi}V(\phi)-\frac{\partial_{\phi}f(\phi)}{4}F_{\mu\nu}^{2}$	$\displaystyle=0.$	(9c)

Now we particularize the above results to the Bjorken flow symmetry [96], which is a geometry commonly used in the context of heavy ion collisions to provide a simplified modeling of the expansion of the medium near the beam axis. The Bjorken symmetry corresponds to boost invariance along the beam line, which we take as the $z$ axis, where the fluid expands at the speed of light, plus translation and $O(2)$ rotation invariance in the transverse $xy$ plane. In practice, this means that in the Bjorken flow one is completely neglecting the dynamics of the medium in the transverse plane, which is the reason why this flow is usually employed only near mid-rapidity in the context of heavy ion collisions. The Bjorken symmetry is more naturally described in terms of the Milne coordinates $(\tau,\xi,x,y)$ , where the proper time $\tau$ and the spacetime rapidity $\xi$ are defined by,

\displaystyle\tau\equiv\sqrt{t^{2}-z^{2}},\qquad\xi\equiv\textrm{arctanh}(v)=\textrm{arctanh}\left(\frac{z}{t}\right)=\frac{1}{2}\ln\left(\frac{t+z}{t-z}\right),

(10)

in terms of which the boundary four dimensional Minkowski metric is written as,

\displaystyle ds^{2}_{\textrm{Mink${}_{4}$}}=-d\tau^{2}+\tau^{2}d\xi^{2}+dx^{2}+dy^{2}.

(11)

The expectation value of the energy-momentum tensor is written as $\langle T^{\mu}_{\nu}\rangle=\mathrm{diag}(-\varepsilon,p_{L},p_{T},p_{T})$ , where for a conformal system, as in the case of the 1RCBH model, the longitudinal and transverse pressures are written as functions of the energy density and its time derivative as follows,

\displaystyle p_{L}(\tau)=-\epsilon(\tau)-\tau\partial_{\tau}\epsilon(\tau),\qquad p_{T}(\tau)=\epsilon(\tau)+\frac{\tau}{2}\partial_{\tau}\epsilon(\tau),

(12)

so that the pressure anisotropy of a conformal fluid undergoing Bjorken flow is given by,

\displaystyle\frac{\Delta p}{\epsilon}\equiv\frac{p_{T}-p_{L}}{\epsilon}=2+\frac{3}{2}\,\tau\frac{\partial_{\tau}\epsilon}{\epsilon}.

(13)

Using five dimensional generalized infalling Eddington-Finkelstein coordinates suited to the holographic implementation of the characteristic formulation of general relativity [63], the ansatz for the bulk Einstein-Maxwell-Dilaton fields compatible with the Bjorken symmetry can be written as follows [85],

\displaystyle ds^{2}=2d\tau\left[dr-A(\tau,r)d\tau\right]+\Sigma(\tau,r)^{2}\left[e^{-2B(\tau,r)}d\xi^{2}+e^{B(\tau,r)}(dx^{2}+dy^{2})\right],\ \ A_{\mu}dx^{\mu}=\Phi(\tau,r)d\tau,\ \ \phi=\phi(\tau,r),

(14)

where $r$ is the radial coordinate of the asymptotically AdS₅ spacetime, whose boundary is at $r\rightarrow\infty$ , where $\tau$ becomes the usual four dimensional proper time of the Bjorken flow. In these coordinates, infalling radial null geodesics satisfy $\tau=\textrm{constant}$ , while outgoing radial null geodesics satisfy $dr/d\tau=A(\tau,r)$ . The metric in Eq. (14) is also invariant under radial diffeomorphism shifts of the form $r\to r+\lambda(\tau)$ , with arbitrary $\lambda(\tau)$ .

By substituting the ansatz (14) into the general Einstein-Maxwell-Dilaton equations of motion (9a) — (9c), one obtains the following set of coupled $1+1$ partial differential equations [85],⁵⁵5In writing down the constraint Eq. (15a) coming from one component of Maxwell’s equations, we made use of the other nontrivial component of Maxwell’s equations, given by Eq. (15b), taking also into account that, from the definitions of $d_{+}$ and $\mathcal{E}$ , one may rewrite $(d_{+}\Phi)^{\prime}=-\partial_{\tau}\mathcal{E}-A^{\prime}\mathcal{E}-A\mathcal{E}^{\prime}$ . Moreover, in writing down the constraint Eq. (15h) coming from Einstein’s equations, we made use of the Hamiltonian constraint (15g) also following from Einstein’s equations.


$\displaystyle\partial_{\tau}\mathcal{E}+A\mathcal{E}^{\prime}+\left(3\frac{d_{+}\Sigma}{\Sigma}+\frac{\partial_{\phi}f}{f}d_{+}\phi\right)\mathcal{E}$	$\displaystyle=0,$	(15a)
$\displaystyle\frac{3\Sigma^{\prime}}{\Sigma}+\frac{\partial_{\phi}f}{f}\phi^{\prime}+\frac{\mathcal{E}^{\prime}}{\mathcal{E}}$	$\displaystyle=0,$	(15b)
$\displaystyle 4\Sigma(d_{+}\phi)^{\prime}+6\phi^{\prime}d_{+}\Sigma+6\Sigma^{\prime}d_{+}\phi+\Sigma\mathcal{E}^{2}\partial_{\phi}f-2\Sigma\partial_{\phi}V$	$\displaystyle=0,$	(15c)
$\displaystyle(d_{+}\Sigma)^{\prime}+\frac{2\Sigma^{\prime}}{\Sigma}d_{+}\Sigma+\frac{\Sigma}{12}\left(2V+f\mathcal{E}^{2}\right)$	$\displaystyle=0,$	(15d)
$\displaystyle\Sigma\,(d_{+}B)^{\prime}+\frac{3}{2}(B^{\prime}d_{+}\Sigma+\Sigma^{\prime}d_{+}B)$	$\displaystyle=0,$	(15e)
$\displaystyle A^{\prime\prime}+\frac{1}{12}\left(18B^{\prime}d_{+}B-\frac{72\Sigma^{\prime}d_{+}\Sigma}{\Sigma^{2}}+6\phi^{\prime}d_{+}\phi-7f\mathcal{E}^{2}-2V\right)$	$\displaystyle=0,$	(15f)
$\displaystyle\Sigma^{\prime\prime}+\frac{\Sigma}{6}\left(3\left(B^{\prime}\right)^{2}+\left(\phi^{\prime}\right)^{2}\right)$	$\displaystyle=0,$	(15g)
$\displaystyle d_{+}(d_{+}\Sigma)+\frac{\Sigma}{2}(d_{+}B)^{2}-A^{\prime}d_{+}\Sigma+\frac{\Sigma}{6}(d_{+}\phi)^{2}$	$\displaystyle=0,$	(15h)

where $X^{\prime}\equiv\partial_{r}X$ is the directional derivative along infalling radial null geodesics, while $d_{+}X\equiv[\partial_{\tau}+A(\tau,r)\partial_{r}]X$ is the directional derivative along outgoing radial null geodesics, and we also defined $\mathcal{E}\equiv-\Phi^{\prime}$ . Notice there are five background functions to be determined, $\{A(\tau,r),\Sigma(\tau,r),B(\tau,r),\phi(\tau,r),\mathcal{E}(\tau,r)\}$ , and there are also five dynamical equations (15b) — (15f), besides three constraint equations corresponding to Eqs. (15a), (15g), and (15h). As we are going to discuss in a moment, Eqs. (15b) — (15g) constitute a nested set of equations of motion which can be systematically integrated using numerical techniques, while the constraint Eqs. (15a) and (15h) can be used to monitor the accuracy of the numerical solutions so obtained, as it will be discussed in the numerical error analysis to be presented in appendix B.

Let us now work out the near-boundary ultraviolet expansions of the bulk fields taking into account the boundary conditions associated to the holographic Bjorken flow. The bulk metric given in Eq. (14) must flow to the AdS₅ geometry in the ultraviolet, which is conformally equivalent to the four dimensional Minkowski metric (11) at the boundary. In order to satisfy this condition, one must impose the following boundary conditions for the metric coefficients,

\displaystyle A(\tau,r\to\infty)\sim\frac{r^{2}}{2},\qquad B(\tau,r\to\infty)\sim-\frac{2}{3}\ln(\tau),\qquad\Sigma(\tau,r\to\infty)\sim\tau^{1/3}r.

(16)

Indeed, by substituting (16) into the five dimensional metric in Eq. (14) one obtains in the generalized infalling Eddington-Finkelstein coordinates,

\displaystyle\lim_{r\to\infty}ds^{2}=ds^{2}_{\textrm{AdS${}_{5}$}}=2d\tau dr+r^{2}ds^{2}_{\textrm{Mink${}_{4}$}},

(17)

where $r^{2}$ is the global conformal factor of AdS₅. Moreover, regarding the bulk dilaton field in the 1RCBH model, one notices from the dilaton potential in Eq. (2) that the conformal dimension of the associated scalar condensate in the dual quantum field theory at the boundary is $\Delta_{\phi}=2$ , consequently, the dilaton field must flow to zero in the ultraviolet as implemented by the following boundary condition,

\displaystyle\phi(\tau,r\to\infty)\sim\frac{\phi_{4-\Delta_{\phi}}(\tau)}{r^{4-\Delta_{\phi}}}=\frac{\phi_{2}(\tau)}{r^{2}}.

(18)

Finally, concerning the nontrivial component of the Maxwell field, it is related to the R-charge chemical potential of the quantum fluid at the boundary according to the following boundary condition,

\displaystyle\Phi(\tau\to\infty,r\to\infty)=\mu.

(19)

Taking into account the above boundary conditions, the near-boundary ultraviolet expansions of the bulk fields are given by [131, 83, 85],


$\displaystyle A(\tau,r)$	$\displaystyle=\frac{1}{2}[r+\lambda(\tau)]^{2}-\partial_{\tau}\lambda(\tau)+\sum_{n=1}^{\infty}\frac{a_{n}(\tau)}{r^{n}},$	(20a)
$\displaystyle B(\tau,r)$	$\displaystyle=-\frac{2}{3}\ln(\tau)+\sum_{n=1}^{\infty}\frac{b_{n}(\tau)}{r^{n}},$	(20b)
$\displaystyle\Sigma(\tau,r)$	$\displaystyle=\tau^{1/3}[r+\lambda(\tau)]+\sum_{n=0}^{\infty}\frac{s_{n}(\tau)}{r^{n}},$	(20c)
$\displaystyle\phi(\tau,r)$	$\displaystyle=\sum_{n=2}^{\infty}\frac{\phi_{n}(\tau)}{r^{n}},$	(20d)
$\displaystyle\Phi(\tau,r)$	$\displaystyle=\Phi_{0}(\tau)+\sum_{n=2}^{\infty}\frac{\Phi_{n}(\tau)}{r^{n}},$	(20e)

where the vast majority of the ultraviolet expansion coefficients $\{a_{n}(\tau),b_{n}(\tau),s_{n}(\tau),\phi_{n}(\tau),\Phi_{n}(\tau)\}$ may be fixed as functions of just a few undetermined coefficients and their time derivatives. This can be accomplished by substituting the above ultraviolet expansions into the Einstein-Maxwell-Dilaton equations of motion and them solving the resulting algebraic equations order by order in powers of $r$ . By working with the ultraviolet expansions up to order $n=8$ one identifies four undetermined coefficients in the aforementioned procedure: $\{a_{2}(\tau),\phi_{2}(\tau),\Phi_{0}(\tau),\Phi_{2}(\tau)\}$ . The values of the ultraviolet coefficients $\{a_{2}(\tau),\phi_{2}(\tau)\}$ at the initial time slice $\tau_{0}$ can be freely chosen since they correspond to two of the four initial data of the 1RCBH model undergoing Bjorken flow, as we are going to discuss in a moment;⁶⁶6More precisely, as we are going to discuss, $\phi_{2}(\tau)$ is the boundary value of the subtracted dilaton field, $\phi_{2}(\tau)=\phi_{s}(\tau,u=0)$ , with the initial profile of the subtracted dilaton field being one of the four initial data of the system. Once this profile is chosen, one automatically has the specification of the initial value $\phi_{2}(\tau_{0})=\phi_{s}(\tau_{0},u=0)$ . once their initial values are chosen, their subsequent time evolutions are determined by numerically solving the nested set of partial differential equations previously derived, as it will be discussed in more details afterwards in the text. Moreover, at order $n=6$ in the aforementioned algebraic procedure, one fixes that $\Phi_{2}(\tau)=c/\tau$ , where $c$ is an undetermined constant; however, from the Bjorken flow hydrodynamic analysis of the local charge conservation equation, $\nabla_{\mu}\langle\hat{J}^{\mu}\rangle=0$ , one concludes that for a conformal fluid, like the 1RCBH model, the time evolution of the R-charge density is given by $\hat{\rho}(\tau)\equiv\langle\hat{J}^{\tau}\rangle=\rho_{0}/\tau$ , where $\rho_{0}$ is a parameter setting the initial charge density of the medium [85]. Since the holographic renormalization procedure fixes that $\hat{\rho}(\tau)\equiv\langle\hat{J}^{\tau}\rangle=-\langle\hat{J}_{\tau}\rangle=-\Phi_{2}(\tau)$ [83], one then concludes that $\Phi_{2}(\tau)=c/\tau=-\rho_{0}/\tau$ , where $\rho_{0}$ is also one of the freely chosen initial data of the system [85]. Finally, none of the other ultraviolet expansion coefficients of the bulk fields directly depend on $\Phi_{0}(\tau)$ , which may be formally identified with the value of the R-charge chemical potential through the boundary condition (19) applied to the series expansion in (20e). However, as it is clear from Eq. (19), the value of $\mu/T$ in the Bjorken flow of the 1RCBH model is not an initial data, and it can only be estimated in the late time evolution of the system [85], when the fluid approaches the equilibrium regime, as it will be explained afterwards in the text.

Before discussing the general schematics to integrate the nested set of $1+1$ partial differential equations of motion of the system, let us consider the radial integration of the Maxwell equation (15b), where we make use of the specific form of Maxwell-dilaton coupling function given in (2),

\displaystyle 3\int\frac{d\Sigma}{\Sigma}-2\sqrt{\frac{2}{3}}\int d\phi+\int\frac{d\mathcal{E}}{\mathcal{E}}=0\,\,\,\Rightarrow\,\,\,\mathcal{E}=-\Phi^{\prime}=e^{\bar{c}}\,\Sigma^{-3}\,e^{2\sqrt{2/3}\,\phi},

(21)

where $\bar{c}$ is a radial integration constant, which can be fixed by evaluating (21) close to the boundary at $r\to\infty$ . In order to do it, we consider the near-boundary expansions for $\mathcal{E}(\tau,r\to\infty)\sim-\partial_{r}(\Phi_{0}+\Phi_{2}r^{-2})=2\Phi_{2}r^{-3}$ and for $e^{\bar{c}}\,\Sigma(\tau,r\to\infty)^{-3}\,e^{2\sqrt{2/3}\,\phi(\tau,r\to\infty)}\sim e^{\bar{c}}\,\tau^{-1}r^{-3}$ , in which case one concludes that $e^{\bar{c}}=2\Phi_{2}(\tau)\tau=-2\rho_{0}$ , where we used the aforementioned relation $\Phi_{2}(\tau)=-\rho_{0}/\tau$ . Consequently,

\displaystyle\mathcal{E}(\tau,r)=-2\rho_{0}\,\Sigma(\tau,r)^{-3}\,e^{2\sqrt{2/3}\,\phi(\tau,r)},

(22)

so that the bulk ‘radial electric field’ $\mathcal{E}=-\Phi^{\prime}$ is completely determined in terms of the metric coefficient $\Sigma$ and the dilaton field $\phi$ .

At this point we may lay down the general reasoning used to numerically integrate the nested set of $1+1$ partial differential equations of motion describing the Bjorken flow dynamics of the 1RCBH model:

i.

On the hypersurface defined at the initial time slice $\tau_{0}$ , one chooses the initial profiles for the metric anisotropy coefficient $B(\tau_{0},r)$ and for the dilaton field $\phi(\tau_{0},r)$ , besides also the initial values for the charge density $\rho_{0}$ and for the dynamical ultraviolet coefficient $a_{2}(\tau_{0})$ ;⁷⁷7As we will discuss in details in the next subsections II.3 and II.4, one also has to choose the initial value for the radial shift function $\lambda(\tau_{0})$ , while its time evolution may be conveniently obtained by requiring that the radial position of the apparent horizon remains fixed during the time evolution of the system [69].
ii.

Next one radially solves the Hamiltonian constraint (15g) to obtain $\Sigma(\tau_{0},r)$ , which at this step automatically fixes the value of $\mathcal{E}(\tau_{0},r)$ through Eq. (22);
iii.

Next one radially solves Eq. (15d) to obtain $d_{+}\Sigma(\tau_{0},r)$ ;
iv.

Next one radially solves Eq. (15e) to obtain $d_{+}B(\tau_{0},r)$ ;
v.

Next one radially solves the dilaton Eq. (15c) to obtain $d_{+}\phi(\tau_{0},r)$ ;
vi.

Next one radially solves Eq. (15f) to obtain $A(\tau_{0},r)$ ;
vii.

At this step, from the definition of the directional derivative along outgoing radial null geodesics, $d_{+}\equiv\partial_{\tau}+A(\tau,r)\partial_{r}$ , one already has $\{\partial_{\tau}B(\tau_{0},r),\partial_{\tau}\phi(\tau_{0},r)\}$ , besides also $\partial_{\tau}a_{2}(\tau_{0})$ (as we will discuss in details in the next subsections II.3 and II.4), which together with the initial profiles chosen for the metric anisotropy and the dilaton field, and the initial value chosen for $a_{2}(\tau)$ , comprise the set of initial conditions required to evolve $\{B(\tau_{0},r),\phi(\tau_{0},r),a_{2}(\tau_{0})\}$ to the next time slice $\tau_{0}+\Delta\tau$ using some discrete integration technique (here we employ the fourth order Adams-Bashforth method to integrate in time, while the radial integrations are performed using the pseudospectral method [138]);
viii.

Repeat the previous steps to obtain all the fields in the current time slice and iteratively evolve the system to the next time slices until reaching any desired end time $\tau_{\textrm{end}}$ for the numerical simulations.

Furthermore, as aforementioned, the constraint Eqs. (15a) and (15h) are used to check the numerical accuracy of the solutions obtained with the above general algorithm, as we shall discuss in appendix B.

II.3 Renormalized 1-point functions, field redefinitions and the apparent horizon

The first few terms in the ultraviolet near-boundary expansions of the bulk fields (20a) — (20e) explicitly read as below,


$\displaystyle A(\tau,r)$	$\displaystyle=\frac{r^{2}}{2}+\lambda r+\frac{\lambda^{2}-2\partial_{\tau}\lambda}{2}+\frac{a_{2}}{r^{2}}+\frac{-4\lambda a_{2}+\partial_{\tau}a_{2}}{2r^{3}}+\mathcal{O}\!\left(r^{-4}\right),$	(23a)
$\displaystyle\Sigma(\tau,r)$	$\displaystyle=\tau^{1/3}r+\frac{1+3\tau\lambda}{3\tau^{2/3}}-\frac{1}{9\tau^{5/3}r}+\frac{5+9\tau\lambda}{81\tau^{8/3}r^{2}}-\frac{20+60\tau\lambda+54\tau^{2}\lambda^{2}+27\tau^{4}\phi_{2}^{2}}{486\tau^{11/3}r^{3}}+\mathcal{O}\!\left(r^{-4}\right),$	(23b)
$\displaystyle B(\tau,r)$	$\displaystyle=-\frac{2}{3}\ln(\tau)-\frac{2}{3\tau r}+\frac{1+2\tau\lambda}{3\tau^{2}r^{2}}-\frac{2+6\tau\lambda+6\tau^{2}\lambda^{2}}{9\tau^{3}r^{3}}$
	$\displaystyle+\frac{6+24\tau\lambda+36\tau^{2}\lambda^{2}+24\tau^{3}\lambda^{3}-36\tau^{4}a_{2}-2\tau^{4}\phi_{2}^{2}-27\tau^{5}\partial_{\tau}a_{2}-3\tau^{5}\phi_{2}\partial_{\tau}\phi_{2}}{36\tau^{4}r^{4}}$
	$\displaystyle-\frac{1}{180\tau^{5}r^{5}}\left[24+120\tau^{2}\lambda^{2}(2+2\tau\lambda+\tau^{2}\lambda^{2})-48\tau^{4}a_{2}+315\tau^{5}\partial_{\tau}a_{2}+35\tau^{5}\phi_{2}\partial_{\tau}\phi_{2}+15\tau^{6}(\partial_{\tau}\phi_{2})^{2}\right.$
	$\displaystyle\left.-\,20\tau\lambda(-6+36\tau^{4}a_{2}+2\tau^{4}\phi_{2}^{2}+27\tau^{5}\partial_{\tau}a_{2}+3\tau^{5}\phi_{2}\partial_{\tau}\phi_{2})+135\tau^{6}\partial^{2}_{\tau}a_{2}+15\tau^{6}\phi_{2}\partial^{2}_{\tau}\phi_{2}\right]+\mathcal{O}\!\left(r^{-6}\right),$	(23c)
$\displaystyle\phi(\tau,r)$	$\displaystyle=\frac{\phi_{2}}{r^{2}}+\frac{-2\lambda\phi_{2}+\partial_{\tau}\phi_{2}}{r^{3}}+\frac{36\tau\lambda^{2}\phi_{2}+\sqrt{6}\tau\phi_{2}^{2}+(3-36\tau\lambda)\partial_{\tau}\phi_{2}+9\tau\partial_{\tau}^{2}\phi_{2}}{12\tau r^{4}}+\mathcal{O}\!\left(r^{-5}\right),$	(23d)
$\displaystyle[d_{+}B](\tau,r)$	$\displaystyle=-\frac{1}{3\tau}+\frac{1}{3\tau^{2}r}-\frac{1+\tau\lambda}{3\tau^{3}r^{2}}+\frac{6+12\tau\lambda+6\tau^{2}\lambda^{2}+36\tau^{4}a_{2}+2\tau^{4}\phi_{2}^{2}+27\tau^{5}\partial_{\tau}a_{2}+3\tau^{5}\phi_{2}\partial_{\tau}\phi_{2}}{18\tau^{4}r^{3}}+\mathcal{O}\!\left(r^{-4}\right),$	(23e)
$\displaystyle[d_{+}\Sigma](\tau,r)$	$\displaystyle=\frac{\tau^{1/3}r^{2}}{2}+\frac{(1+3\tau\lambda)r}{3\tau^{2/3}}+\frac{-1+2\tau\lambda+3\tau^{2}\lambda^{2}}{6\tau^{5/3}}+\frac{10}{81\tau^{8/3}r}+\frac{-100-120\tau\lambda+972\tau^{4}a_{2}+81\tau^{4}\phi_{2}^{2}}{972\tau^{11/3}r^{2}}+\mathcal{O}\!\left(r^{-3}\right),$	(23f)
$\displaystyle[d_{+}\phi](\tau,r)$	$\displaystyle=-\frac{\phi_{2}}{r}+\mathcal{O}\!\left(r^{-2}\right).$	(23g)

The holographic renormalization procedure for the 1RCBH model [83] allows one to write down the renormalized 1-point functions of the dual quantum field theory, namely, the expectation values of the boundary energy-momentum tensor $\langle T_{\mu\nu}\rangle$ , the $U(1)$ vector current $\langle J_{\mu}\rangle$ , and the scalar condensate $\langle O_{\phi}\rangle$ , in terms of the bulk gravity ultraviolet coefficients $\{a_{2}(\tau),\phi_{2}(\tau),\Phi_{2}(\tau)=-\rho_{0}/\tau\}$ and their time derivatives. In order to obtain those relations, one first needs to relate the holographic radial coordinates $r$ and $\bar{\rho}$ written,⁸⁸8We use here the bar to denote the radial Fefferman-Graham coordinate $\bar{\rho}$ in order to avoid confusing its notation with the charge density $\rho(\tau)$ . respectively, in the generalized infalling Eddington-Finkelstein coordinates, where we developed our previously discussed analysis of the holographic Bjorken flow dynamics, and in the Fefferman-Graham coordinates used in the holographic renormalization procedure [83]. For the holographic Bjorken flow this relation is explicitly given by [91],

\displaystyle r(\bar{\rho})=\frac{1}{\sqrt{\bar{\rho}}}-\frac{a_{2}(\tau)\bar{\rho}^{3/2}}{4}-\frac{\partial_{\tau}a_{2}(\tau)\bar{\rho}^{2}}{10}+\mathcal{O}(\bar{\rho}^{5/2}).

(24)

By substituting Eq. (24) into the near-boundary expansions (23a) — (23g), one identifies the relevant ultraviolet coefficients in Fefferman-Graham coordinates entering the holographic renormalization formulas for the 1-point functions of the the dual quantum field theory at the boundary [83]. The final results for the renormalized 1-point functions of the 1RCBH plasma undergoing Bjorken flow read as follows [85],


$\displaystyle\hat{\epsilon}(\tau)$	$\displaystyle\equiv\kappa_{5}^{2}\langle T_{\tau\tau}\rangle(\tau)=-3a_{2}(\tau)-\frac{1}{6}\phi_{2}(\tau)^{2},$	(25a)
$\displaystyle\hat{p}_{T}(\tau)$	$\displaystyle\equiv\kappa_{5}^{2}\langle T_{x}^{x}\rangle(\tau)=\kappa_{5}^{2}\langle T_{y}^{y}\rangle(\tau)=-3a_{2}(\tau)-\frac{1}{6}\phi_{2}(\tau)^{2}-\frac{3}{2}\tau\partial_{\tau}a_{2}(\tau)-\frac{1}{6}\tau\phi_{2}(\tau)\partial_{\tau}\phi_{2}(\tau),$	(25b)
$\displaystyle\hat{p}_{L}(\tau)$	$\displaystyle\equiv\kappa_{5}^{2}\langle T_{\xi}^{\xi}\rangle(\tau)=3a_{2}(\tau)+\frac{1}{6}\phi_{2}(\tau)^{2}+3\tau\partial_{\tau}a_{2}(\tau)+\frac{1}{3}\tau\phi_{2}(\tau)\partial_{\tau}\phi_{2}(\tau),$	(25c)
$\displaystyle\hat{\rho}(\tau)$	$\displaystyle\equiv\kappa_{5}^{2}\langle J^{\tau}\rangle(\tau)=-\Phi_{2}(\tau)=\frac{\rho_{0}}{\tau},$	(25d)
$\displaystyle\langle\hat{O}_{\phi}\rangle(\tau)$	$\displaystyle\equiv\kappa_{5}^{2}\langle\mathcal{O}_{\phi}\rangle(\tau)=-\phi_{2}(\tau).$	(25e)

From the above holographic formulas for the physical observables of the 1RCBH plasma, one can immediately confirm that the boundary system is conformal, since it has vanishing trace anomaly: $\langle\hat{T}_{\mu}^{\mu}\rangle=-\hat{\epsilon}+\hat{p}_{L}+2\hat{p}_{T}=0$ .

In order to facilitate the numerical implementation of the general steps previously discussed to integrate the equations of motions of the system, we make use of field redefinitions corresponding to subtractions of the leading order near-boundary terms in the ultraviolet expansions of the bulk fields (23a) — (23g), such that the boundary values of the subtracted fields go to radial constants at the boundary,


$\displaystyle u^{2}A_{s}(\tau,u)$	$\displaystyle\equiv A(\tau,u)-\frac{1}{2u^{2}}-\frac{\lambda}{u}-\frac{\lambda^{2}}{2}+\partial_{\tau}\lambda,$	(26a)
$\displaystyle u^{3}\Sigma_{s}(\tau,u)$	$\displaystyle\equiv\Sigma(\tau,u)-\frac{\tau^{1/3}}{u}-\frac{1+3\tau\lambda}{3\tau^{2/3}}+\frac{u}{9\tau^{5/3}}-\frac{(5+9\tau\lambda)u^{2}}{81\tau^{8/3}},$	(26b)
$\displaystyle u^{4}B_{s}(\tau,u)$	$\displaystyle\equiv B(\tau,u)+\frac{2\ln(\tau)}{3}+\frac{2u}{3\tau}-\frac{(1+2\tau\lambda)u^{2}}{3\tau^{2}}+\frac{(2+6\tau\lambda+6\tau^{2}\lambda^{2})u^{3}}{9\tau^{3}},$	(26c)
$\displaystyle u^{2}\phi_{s}(\tau,u)$	$\displaystyle\equiv\phi(\tau,u),$	(26d)
$\displaystyle u^{3}[d_{+}B]_{s}(\tau,u)$	$\displaystyle\equiv[d_{+}B](\tau,u)+\frac{1}{3\tau}-\frac{u}{3\tau^{2}}+\frac{(1+\tau\lambda)u^{2}}{3\tau^{3}},$	(26e)
$\displaystyle u^{2}[d_{+}\Sigma]_{s}(\tau,u)$	$\displaystyle\equiv[d_{+}\Sigma](\tau,u)-\frac{\tau^{1/3}}{2u^{2}}-\frac{1+3\tau\lambda}{3\tau^{2/3}u}-\frac{-1+2\tau\lambda+3\tau^{2}\lambda^{2}}{6\tau^{5/3}}-\frac{10u}{81\tau^{8/3}},$	(26f)
$\displaystyle u[d_{+}\phi]_{s}(\tau,u)$	$\displaystyle\equiv[d_{+}\phi](\tau,u),$	(26g)
$\displaystyle\mathcal{E}_{s}(\tau,u)$	$\displaystyle=\mathcal{E}(\tau,u),$	(26h)

where we defined the new radial coordinate $u\equiv 1/r$ , in terms of which the boundary lies at $u=0$ . In terms of the above subtracted fields one can easily see that the corresponding boundary values are given by the following radial constants [85],


$\displaystyle A_{s}(\tau,u=0)$	$\displaystyle=a_{2},$	(27a)
$\displaystyle\Sigma_{s}(\tau,u=0)$	$\displaystyle=-\frac{20+60\tau\lambda+54\tau^{2}\lambda^{2}+27\tau^{4}\phi_{2}^{2}}{486\tau^{11/3}},$	(27b)
$\displaystyle B_{s}(\tau,u=0)$	$\displaystyle=\frac{6+24\tau\lambda+36\tau^{2}\lambda^{2}+24\tau^{3}\lambda^{3}-36\tau^{4}a_{2}-2\tau^{4}\phi_{2}^{2}-27\tau^{5}\partial_{\tau}a_{2}-3\tau^{5}\phi_{2}\partial_{\tau}\phi_{2}}{36\tau^{4}},$	(27c)
$\displaystyle\phi_{s}(\tau,u=0)$	$\displaystyle=\phi_{2},$	(27d)
$\displaystyle[d_{+}B]_{s}(\tau,u=0)$	$\displaystyle=\frac{6+12\tau\lambda+6\tau^{2}\lambda^{2}+36\tau^{4}a_{2}+2\tau^{4}\phi_{2}^{2}+27\tau^{5}\partial_{\tau}a_{2}+3\tau^{5}\phi_{2}\partial_{\tau}\phi_{2}}{18\tau^{4}},$	(27e)
$\displaystyle[d_{+}\Sigma]_{s}(\tau,u=0)$	$\displaystyle=\frac{-100-120\tau\lambda+972\tau^{4}a_{2}+81\tau^{4}\phi_{2}^{2}}{972\tau^{11/3}},$	(27f)
$\displaystyle[d_{+}\phi]_{s}(\tau,u=0)$	$\displaystyle=-\phi_{2}.$	(27g)

The equations of motion to be numerically solved as functions of $(\tau,u)$ are then given by the original equations of motion rewritten in terms of the subtracted fields, whose boundary values were obtained above.

Before proceeding to a discussion on the numerics in the next subsection II.4, we briefly discuss the calculation of the apparent horizon of the far from equilibrium black hole solutions.

The apparent horizon is the outermost trapped null surface inside the event horizon, separating a spacetime region where the geodesics are directed outward with light rays moving outward and a region where the light rays along the same geodesics move inward. Consequently, inside the apparent horizon all light rays move inward — see e.g. Fig. 2 of Ref. [63] for a clear illustration. The apparent horizon converges to the event horizon at late times, when the black hole geometry approaches equilibrium. Since the apparent horizon lies inside the event horizon, by cutting off the radial integration of the bulk equations of motion at some position inside the apparent horizon, one assures that the radial domain of the bulk geometry causally connected to observers at the boundary is adequately taken into account and no physical information is lost in such integration procedure.

In the holographic Bjorken flow, the radial position of the apparent horizon within the bulk generally presents wide fluctuations during the time evolution of the system. This is clearly inconvenient if one wants to cutoff the radial integration at some value of the holographic radial coordinate and hold this cutoff fixed for any value of time, since these wide fluctuations may eventually lead the radial cutoff to lie beyond, instead of behind the apparent horizon for some values of time, what could lead to loss of information and inaccurate physical results. Fortunately, it is possible to fix this issue by using the radial diffeomorphism shift invariance mentioned below Eq. (14), which involves the function $\lambda(\tau)$ . Indeed, since it is an arbitrary function of time, $\lambda(\tau)$ may be chosen in a way such as to keep the radial position of the apparent horizon fixed during the time evolution of the system [69].

For any metric field of the form shown in Eq. (14), by assuming that the radial position of the apparent horizon $r_{\textrm{AH}}$ remains constant in time, one calculates it as the value of the radial coordinate which satisfies the following condition [69],

\displaystyle[d_{+}\Sigma](\tau,r_{\textrm{AH}})=0.

(28)

By requiring that $\partial_{\tau}r_{\textrm{AH}}(\tau)=0$ and that Eq. (28) remain valid at all times, it follows that $\partial_{\tau}[d_{+}\Sigma](\tau,r_{\textrm{AH}})=0$ , implying that $d_{+}[d_{+}\Sigma](\tau,r_{\textrm{AH}})=A(\tau,r_{\textrm{AH}})\partial_{r}[d_{+}\Sigma](\tau,r_{\textrm{AH}})$ . Substituting this result into the constraint Eq. (15h), and then combining it with other components of Einstein’s equations, one obtains that,

\displaystyle A(\tau,u_{\textrm{AH}})=\frac{6([d_{+}B](\tau,u_{\textrm{AH}}))^{2}+2([d_{+}\phi](\tau,u_{\textrm{AH}}))^{2}}{2V+f\mathcal{E}^{2}},

(29)

where we used Eq. (28). By using now Eq. (26a) evaluated at the apparent horizon position $u=u_{\textrm{AH}}$ , it follows that,

\displaystyle\partial_{\tau}\lambda(\tau)

\displaystyle=u_{\textrm{AH}}^{2}A_{s}(\tau,u_{\textrm{AH}})+\frac{1}{2u_{\textrm{AH}}^{2}}+\frac{\lambda(\tau)}{u_{\textrm{AH}}}+\frac{\lambda^{2}(\tau)}{2}-A(\tau,u_{\textrm{AH}}).

(30)

Once one chooses the initial value $\lambda(\tau_{0})$ , the function $\lambda(\tau)$ is evolved in time by using Eq. (30), which keeps the radial position of the apparent horizon fixed during the time evolution. In practice, as done in [90, 91], we choose $\lambda(\tau_{0})=0$ and solve Eq. (28) using Eq. (26b) with the Newton-Raphson method. In this way, the apparent horizon remains fixed, within some numerical tolerance, in its initial position calculated with the initial condition $\lambda(\tau_{0})=0$ .

II.4 Numerical solver

In order to numerically integrate the nested set of $1+1$ partial differential equations of motion of the system, we discretize the radial and time directions. We employ the pseudospectral method [138] to integrate in the radial direction and make use of the fourth order Adams-Bashforth method to integrate in the time direction. In both cases, the general algorithmic steps are the same discussed e.g. in [83, 91], and we refer the interested reader to those works for the details.

Here we supplement the information concerning the extra steps required to implement these algorithms in the case of the 1RCBH model undergoing Bjorken flow. The complete set of freely chosen initial data needed to be specified on the initial time slice $\tau_{0}$ is given by $\{B_{s}(\tau_{0},u),\phi_{s}(\tau_{0},u),a_{2}(\tau_{0}),\rho_{0};\lambda(\tau_{0})\}$ . In order to evolve the system to the next time slices, one also needs to obtain the expressions for the time derivatives of these data. As discussed at the end of the previous subsection II.3, we set here $\lambda(\tau_{0})=0$ and $\partial_{\tau}\lambda(\tau)$ is given by Eq. (30). On the other hand, the physically different initial data will be determined by different choices for the subset $\{B_{s}(\tau_{0},u),\phi_{s}(\tau_{0},u),a_{2}(\tau_{0}),\rho_{0}\}$ , which will be discussed in the next section III. It remains, therefore, the task of obtaining the expressions for the time derivatives $\{\partial_{\tau}B_{s}(\tau,u),\partial_{\tau}\phi_{s}(\tau,u),\partial_{\tau}a_{2}(\tau)\}$ required to evolve the system in time.

The expressions for $\partial_{\tau}B_{s}$ and $\partial_{\tau}\phi_{s}$ can be derived from the expressions for $d_{+}B=\partial_{\tau}B+A\partial_{r}B$ and $d_{+}\phi=\partial_{\tau}\phi+A\partial_{r}\phi$ rewritten in terms of the subtracted fields and the radial coordinate $u=1/r$ . After some tedious algebraic manipulations, one arrives at the following results,

$\displaystyle\partial_{\tau}B_{s}(\tau,u)$	$\displaystyle=\frac{[d_{+}B]_{s}}{u}-\frac{2}{3\tau^{4}u}-\frac{2A_{s}}{3\tau}+\frac{2uA_{s}}{3\tau^{2}}-\frac{2u^{2}A_{s}}{3\tau^{3}}+4u^{3}A_{s}B_{s}+\frac{2B_{s}}{u}+\frac{B_{s}^{\prime}}{2}+u^{4}A_{s}B_{s}^{\prime}$
	$\displaystyle+\,\left(4B_{s}-\frac{2}{u\tau^{3}}+\frac{4uA_{s}}{3\tau}-\frac{2u^{2}A_{s}}{\tau^{2}}+uB_{s}^{\prime}\right)\lambda+\left(\!-\frac{1}{3\tau^{3}}\!-\!\frac{7}{3\tau^{2}u}\!+\!2uB_{s}\!-\!\frac{2u^{2}A_{s}}{\tau}+\frac{u^{2}B_{s}^{\prime}}{2}\right)\lambda^{2}$
	$\displaystyle-\,\left(\frac{1}{\tau^{2}}+\frac{4}{3\tau u}\right)\lambda^{3}-\frac{\lambda^{4}}{\tau}+\left(\frac{2}{3\tau^{3}}-4uB_{s}-u^{2}B_{s}^{\prime}+\frac{2\lambda}{\tau^{2}}+\frac{2\lambda^{2}}{\tau}\right)\partial_{\tau}\lambda,$	(31)
$\displaystyle\partial_{\tau}\phi_{s}(\tau,u)$	$\displaystyle=\frac{[d_{+}\phi]_{s}}{u}+\frac{\phi_{s}^{\prime}}{2}+u^{4}A_{s}\phi_{s}^{\prime}+\frac{\phi_{s}}{u}+2u^{3}A_{s}\phi_{s}+\left(u\phi_{s}^{\prime}+2\phi_{s}\right)\lambda+\left(\frac{u^{2}\phi_{s}^{\prime}}{2}+u\phi_{s}\right)\lambda^{2}-\left(2u\phi_{s}+u^{2}\phi_{s}^{\prime}\right)\partial_{\tau}\lambda,$	(32)

where $X_{s}^{\prime}(\tau,u)\equiv\partial_{u}X_{s}(\tau,u)$ is evaluated at any constant time slice by applying the pseudospectral finite differentiation matrix [138] to the numerically known solution $X_{s}(\tau,u)$ , which is expressed as a vector with $N$ components corresponding to the values of $X_{s}(\tau,u)$ on top of the $N$ collocation points of the Chebyshev-Gauss-Lobatto radial grid [83, 91].

Regarding the expression for $\partial_{\tau}a_{2}(\tau)$ , it can be determined directly from Eq. (27c) once $\partial_{\tau}\phi_{2}(\tau)$ is known, which in turn follows from Eqs. (23d) and (26d),

\displaystyle\partial_{\tau}\phi_{2}(\tau)=\phi_{s}^{\prime}(\tau,u=0)+2\lambda(\tau)\phi_{2}(\tau).

(33)

In the algorithms employed in [83, 85, 90, 91] to integrate in the radial direction, the boundary values of the fields are calculated separately, which can be done here by using Eqs. (27a) — (27g) and the following results also coming from the near-boundary expansions (23a) — (23g),

$\displaystyle\partial_{\tau}B_{s}(\tau,u=0)$	$\displaystyle=-\frac{2}{3\tau^{5}}-\frac{2\lambda}{\tau^{4}}-\frac{2\lambda^{2}}{\tau^{3}}-\frac{2\lambda^{3}}{3\tau^{2}}+\frac{2\partial_{\tau}\lambda}{3\tau^{3}}+\frac{2\lambda\partial_{\tau}\lambda}{\tau^{2}}+\frac{2\lambda^{2}\partial_{\tau}\lambda}{\tau}-\frac{7\partial_{\tau}a_{2}}{4}-\frac{7\phi_{2}\partial_{\tau}\phi_{2}}{36}-\frac{\tau(\partial_{\tau}\phi_{2})^{2}}{12}$
	$\displaystyle-\,\frac{3\tau\partial^{2}_{\tau}a_{2}}{4}-\frac{\tau\phi_{2}\partial^{2}_{\tau}\phi_{2}}{12},$	(34)
$\displaystyle\partial_{\tau}\phi_{s}(\tau,u=0)$	$\displaystyle=\partial_{\tau}\phi_{2},$	(35)

where, in particular, the following second derivatives are obtained from Eqs. (23c) and (23d), respectively,

$\displaystyle\partial^{2}_{\tau}a_{2}(\tau)$	$\displaystyle=-\frac{4B^{\prime}_{s}(\tau,u=0)}{3\tau}-\frac{8}{45\tau^{6}}-\frac{8\lambda}{9\tau^{5}}-\frac{16\lambda^{2}}{9\tau^{4}}-\frac{16\lambda^{3}}{9\tau^{3}}-\frac{8\lambda^{4}}{9\tau^{2}}+\frac{16a_{2}}{45\tau^{2}}+\frac{16\lambda a_{2}}{3\tau}+\frac{8\lambda\phi_{2}^{2}}{27\tau}$
	$\displaystyle-\,\frac{7\partial_{\tau}a_{2}}{3\tau}+4\lambda\partial_{\tau}a_{2}-\frac{7\phi_{2}\partial_{\tau}\phi_{2}}{27\tau}+\frac{4\lambda\phi_{2}\partial_{\tau}\phi_{2}}{9}-\frac{(\partial_{\tau}\phi_{2})^{2}}{9}-\frac{\phi_{2}\partial^{2}_{\tau}\phi_{2}}{9},$	(36)
$\displaystyle\partial^{2}_{\tau}\phi_{2}(\tau)$	$\displaystyle=\frac{2\phi^{\prime\prime}_{s}(\tau,u=0)}{3}-4\lambda^{2}\phi_{2}-\sqrt{\frac{2}{27}}\,\phi_{2}^{2}-\frac{\partial_{\tau}\phi_{2}}{3\tau}+4\lambda\partial_{\tau}\phi_{2}.$	(37)

We close this section by remarking that for the numerical calculations carried out in the present work we used $N=21$ collocation points in the radial grid and a time step size of $\Delta\tau=12\times 10^{-5}$ . However, for three specific initial conditions to be discussed next we needed to increase the number of collocation points to $N\sim 30$ in order to eliminate numerically spurious oscillations in the normalized scalar condensate at early times.

III Non-equilibrium entropy, initial data, physical observables and energy conditions

In this section, we discuss the holographic formula for the non-equilibrium entropy density, the form of the initial data $\{B_{s}(\tau_{0},u),\phi_{s}(\tau_{0},u),a_{2}(\tau_{0}),\rho_{0}\}$ which will be analyzed in the present work, besides presenting the conventions we will use to plot the dimensionless ratios for the physical observables, and also the form of the dominant and weak energy conditions for a conformal fluid undergoing Bjorken flow.

Concerning the holographic calculation of the non-equilibrium entropy density, it proceeds as follows [63, 90, 91]. The area of the apparent horizon reads,

\displaystyle A_{\textrm{AH}}(\tau)=\int d^{3}x\sqrt{-g}\biggr{|}_{u=u_{\textrm{AH}}}=\int dxdyd\xi\sqrt{-g}\biggr{|}_{u=u_{\textrm{AH}}}=\sqrt{-g}\biggr{|}_{u=u_{\textrm{AH}}}\mathcal{A}=|\Sigma(\tau,u_{\textrm{AH}})|^{3}\mathcal{A},

(38)

where $\mathcal{V}(\tau)=\tau\mathcal{A}=\tau\int dxdyd\xi$ is the expanding volume of the medium in Bjorken flow. Analogously to the Bekenstein-Hawking relation, one obtains for the non-equilibrium holographic entropy the following relation, written in terms of the area of the apparent horizon,

\displaystyle S_{\textrm{AH}}(\tau)=\frac{A_{\textrm{AH}}(\tau)}{4G_{5}}=\frac{2\pi|\Sigma(\tau,u_{\textrm{AH}})|^{3}\mathcal{A}}{\kappa_{5}^{2}},

(39)

and the entropy density is, therefore,

\displaystyle\hat{s}_{\textrm{AH}}(\tau)\equiv\kappa_{5}^{2}\,s_{\textrm{AH}}(\tau)=\kappa_{5}^{2}\,\frac{S_{\textrm{AH}}(\tau)}{\mathcal{V}(\tau)}=\frac{2\pi|\Sigma(\tau,u_{\textrm{AH}})|^{3}}{\tau},

(40)

where $\Sigma(\tau,u_{\textrm{AH}})$ is calculated in terms of the numerical result for $\Sigma_{s}(\tau,u_{\textrm{AH}})$ through Eq. (26b) evaluated at the radial location of the apparent horizon.

Regarding the initial profile for the subtracted metric anisotropy coefficient, we are going to work with the following general form [90, 91] (which is considerably broader than the forms considered e.g. in [85, 139]),

B_{s}(\tau_{0},u)=\Omega_{1}\cos(\gamma_{1}u)+\Omega_{2}\tan(\gamma_{2}u)+\Omega_{3}\sin(\gamma_{3}u)+\sum_{i=0}^{5}\beta_{i}u^{i}+\,\frac{\alpha}{u^{4}}\left[-\frac{2}{3}\ln\left(1+\frac{u}{\tau_{0}}\right)+\frac{2u^{3}}{9\tau_{0}^{3}}-\frac{u^{2}}{3\tau_{0}^{2}}+\frac{2u}{3\tau_{0}}\right],

(41)

such that one needs to choose the values of the parameters $\{\Omega_{i},\gamma_{i},\beta_{i},\alpha\}$ in (41), and also the values of the initial data $\{\phi_{s}(\tau_{0},u),a_{2}(\tau_{0}),\rho_{0}\}$ in order to fully specify a given initial condition. We take $\tau_{0}=0.2$ as the initial time of our numerical simulations and evolve each initial condition up to $\tau_{\textrm{end}}=35$ . This end time used in our present numerical simulations is about five times larger than considered in previous works [85, 90, 91], what is needed in order to be able to see the effective thermalization of the scalar condensate at late times in the evolution of the system. The set of parameters $\{\Omega_{i},\gamma_{i},\beta_{i},\alpha\}$ in (41) which we are going to explore in this paper is provided in Table 1. With these different profiles for $B_{s}(\tau_{0},u)$ it is possible to generate many physically different possibilities for the time evolution of the 1RCBH plasma expanding according to the Bjorken flow dynamics.

$B_{s}\#$	$\Omega_{1}$	$\gamma_{1}$	$\Omega_{2}$	$\gamma_{2}$	$\Omega_{3}$	$\gamma_{3}$	$\beta_{0}$	$\beta_{1}$	$\beta_{2}$	$\beta_{3}$	$\beta_{4}$	$\beta_{5}$	$\alpha$
1	0	0	0	0	0	0	0.5	-0.5	0.4	0.2	-0.3	0.1	1
2	0	0	0	0	0	0	-0.2	-0.5	0.3	0.1	-0.2	0.4	1
3	0	0	0	0	0	0	0.1	-0.4	0.3	0	-0.1	0	1
4	0	0	1	1	0	0	0	0	0	0	0	0	1
5	0	0	0	0	0	0	-0.2	-0.5	0	0	0	0	1
6	0	0	0	0	0	0	-0.2	-0.4	0	0	0	0	1
7	0	0	0	0	0	0	-0.2	-0.6	0	0	0	0	1
8	0	0	0	0	0	0	-0.3	-0.5	0	0	0	0	1
9	0	0	0	0	1	8	0	0	0	0	0	0	1
10	1	8	0	0	0	0	-0.2	-0.5	0	0	0	0	1
11	0.5	8	0	0	0	0	-0.2	-0.5	0	0	0	0	1

Table 1: Set of parameters for the initial profile of the subtracted metric anisotropy (41) analyzed in this work.

In what concerns the initial profile for the subtracted dilaton field, we are going to consider the following general form,

\phi_{s}(\tau_{0},u)=\sum_{i=0}^{3}a_{i}u^{i}+\alpha_{0}\,e^{-u^{2}/\sigma}+\alpha_{1}\cos(\omega_{1}u)+\alpha_{2}\sin(\omega_{2}u),

(42)

The set of parameters $\{a_{i},\alpha_{i},\sigma,\omega_{i}\}$ in (42) which we are going to explore in this paper is provided in Table 2. ⁹⁹9Notice that $\phi_{s}1$ in Table 2 is the trivial dilaton profile, $\phi_{s}(\tau_{0},u)=0$ .

$\phi_{s}\#$	$a_{0}$	$a_{1}$	$a_{2}$	$a_{3}$	$\alpha_{0}$	$\alpha_{1}$	$\alpha_{2}$	$\sigma$	$\omega_{1}$	$\omega_{2}$
1	0	0	0	0	0	0	0	2	0	0
2	0	1	1	0	0	0	0	2	0	0
3	-0.8	1.1	0.6	0	0	0	0	2	0	0
4	0.7	-0.9	0	0.4	0	0	0	2	0	0
5	0	0	0	0	0.5	0	0	2	0	0
6	0	0	0	0	0	1	0	2	1	0
7	0	0	0	0	0	0	1	2	0	5

Table 2: Set of parameters for the initial profile of the subtracted dilaton field (42) analyzed in this work.

As discussed before, the initial data $\rho_{0}$ sets the initial value of the charge density of the medium, $\hat{\rho}(\tau_{0})=\rho_{0}/\tau_{0}$ . On the other hand, $a_{2}(\tau_{0})$ together with $\phi_{2}(\tau_{0})=\phi_{s}(\tau_{0},u=0)$ specify the initial value of the energy density of the fluid, $\hat{\epsilon}(\tau_{0})$ , according to Eq. (25a). In section IV, we will present our results for the physical observables of the 1RCBH plasma undergoing Bjorken flow, using as initial data the profiles for $B_{s}(\tau_{0},u)$ given in Table 1 with variations on the initial charge density of the medium, keeping fixed its initial energy density and setting $\phi_{s}(\tau_{0},u)=0$ . In section V, we vary instead the initial energy density of the medium, while keeping fixed its initial charge density and setting again $\phi_{s}(\tau_{0},u)=0$ . In section VI, we shall analyze the results obtained with the different profiles for $\phi_{s}(\tau_{0},u)$ given in Table 2, while keeping fixed the remaining initial data.

Regarding the physical observables we shall consider in this work, they will be expressed through the following dimensionless ratios: $\Delta\hat{p}/\hat{\epsilon}=(\hat{p}_{T}-\hat{p}_{L})/\hat{\epsilon}$ , which is obtained from Eqs. (25a) — (25c); $\langle\hat{O}_{\phi}\rangle/\hat{\epsilon}^{1/2}$ , which is obtained from Eqs. (25a) and (25e); $\hat{\rho}^{4/3}/\hat{\epsilon}$ , which is obtained from Eqs. (25a) and (25d); while for the entropy, we will use the ratio $\hat{s}_{\textrm{AH}}^{4/3}/\hat{\epsilon}$ , which is obtained from Eqs. (25a) and (40), only in the late time physical consistency analysis of the numerical solutions to be discussed in appendix A. Indeed, whilst $\hat{s}_{\textrm{AH}}^{4/3}/\hat{\epsilon}$ is important to check the late time convergence of the entropy density to the thermodynamically stable branch of equilibrium black hole solutions previously discussed in section II.1, since the energy density is a nontrivial function of time, one needs a different dimensionless normalization for the entropy density $s_{\textrm{AH}}(\tau)$ in order to be able to follow the actual time evolution of the entropy function $S_{\textrm{AH}}(\tau)$ , which is directly related to the area of the apparent horizon $A_{\textrm{AH}}(\tau)$ through Eq. (39). This is important in order to check the validity of the second law of thermodynamics during the time evolution of the fluid, which is associated here with the fact that the area of the apparent horizon should not decrease as a function of time during the evolution of the system; moreover, it is also important to directly track the time evolution of the area of the apparent horizon because flat regions for this observable are a direct measure of regions in the time evolution of the fluid with zero entropy production [90, 91]. Thus, for the main physical analyses present in this paper we are going to describe the time evolution of the entropy of the medium through the following dimensionless ratio,

\displaystyle\frac{\tau\hat{s}_{\textrm{AH}}(\tau)}{\Lambda^{2}}=\frac{\hat{S}_{\textrm{AH}}(\tau)}{\mathcal{A}\Lambda^{2}}=\frac{2\pi A_{\textrm{AH}}(\tau)}{\mathcal{A}\Lambda^{2}}=\frac{2\pi|\Sigma(\tau,u_{\textrm{AH}})|^{3}}{\Lambda^{2}},

(43)

where $\Lambda$ is an energy scale, extracted for each initial condition, through a fit of the late time result for the full numerical energy density to its corresponding analytical hydrodynamic Navier-Stokes result [85]. This energy scale shall also be used to define the following effective dimensionless time measure, $\omega_{\Lambda}(\tau)\equiv\tau T_{\textrm{eff}}(\tau)$ , where we follow [10, 90, 91] and take the ‘effective non-equilibrium temperature’ (at zero density), $T_{\textrm{eff}}(\tau)$ , to be given by the third-order hydrodynamic truncation for the energy density of the pure thermal SYM plasma [140],

T_{\textrm{3rd}}^{(0)}(\tau)=\frac{\Lambda}{(\Lambda\tau)^{1/3}}\left[1-\frac{1}{6\pi(\Lambda\tau)^{2/3}}+\frac{-1+\ln 2}{36\pi^{2}(\Lambda\tau)^{4/3}}+\,\frac{-21+2\pi^{2}+51\ln 2-24\ln^{2}2}{1944\pi^{3}(\Lambda\tau)^{2}}\right].

(44)

All the aforementioned dimensionless ratios for the physical observables of the 1RCBH plasma undergoing Bjorken flow will be numerically interpolated as functions of the dimensionless time measure $\omega_{\Lambda}(\tau)\equiv\tau T_{\textrm{3rd}}^{(0)}(\tau)$ in order to plot our main results in sections IV and V.

Regarding the value of $\mu/T$ for each initial condition, as discussed before, it is not an initial data in the Bjorken flow of the 1RCBH plasma, instead it is extracted from the latest time result for $[\hat{\rho}^{4/3}/\hat{\epsilon}](\tau)$ , which is matched to the corresponding thermodynamically stable equilibrium result for $[\hat{\rho}^{4/3}/\hat{\epsilon}]_{\textrm{(eq.)}}$ discussed in section II.1. Clearly, for a reliable estimate of the value of $\mu/T$ , this matching procedure needs to be done when $[\hat{\rho}^{4/3}/\hat{\epsilon}](\tau)$ is almost stabilized in a constant value at late times, and the result will be more precise the larger is the end time $\tau_{\textrm{end}}$ of the numerical simulations.

Closely related to the above discussion, several important analytical consistency checks of our numerical results can be performed. Since all the near-equilibrium hydrodynamic results and also all the equilibrium thermodynamic results are functions of $\mu/T$ , and not of the specific initial conditions of the far from equilibrium fluid, once the value of $\mu/T$ for a given initial condition, extracted from the late time analysis of $[\hat{\rho}^{4/3}/\hat{\epsilon}](\tau)$ , is plugged into the analytical expressions for $[\langle\hat{O}_{\phi}\rangle/\hat{\epsilon}^{1/2}]_{\textrm{(eq.)}}$ and $[\hat{s}_{\textrm{AH}}^{4/3}/\hat{\epsilon}]_{\textrm{(eq.)}}$ in the thermodynamically stable branches of equilibrium black hole solutions, one has the asymptotic values that should be attained at late times by the dynamical observables $[\langle\hat{O}_{\phi}\rangle/\hat{\epsilon}^{1/2}](\tau)$ and $[\hat{s}_{\textrm{AH}}^{4/3}/\hat{\epsilon}](\tau)$ ; moreover, for $[\Delta\hat{p}/\hat{\epsilon}](\tau)$ , one should obtain at late times convergence of the numerical results to the corresponding analytical hydrodynamic Navier-Stokes (NS) results coming from Eq. (13),

\displaystyle\left[\frac{\Delta\hat{p}}{\hat{\epsilon}}\right]_{\textrm{NS}}(\tau,\mu/T)=2+\frac{3}{2}\,\tau\frac{\partial_{\tau}\hat{\epsilon}_{\textrm{NS}}(\tau,\mu/T)}{\hat{\epsilon}_{\textrm{NS}}(\tau,\mu/T)},

(45)

with $\hat{\epsilon}_{\textrm{NS}}(\tau,\mu/T)$ given by Eq. (12) of [85] evaluated at the same value of $\mu/T$ obtained from the late time analysis of $[\hat{\rho}^{4/3}/\hat{\epsilon}](\tau)$ . We indeed consistently see, for all the initial conditions, the late time convergence of all these physical observables to their corresponding analytical hydrodynamic or thermodynamic results evaluated at the same value of $\mu/T$ . These outcomes constitute independent and highly nontrivial physical consistency checks of our numerical simulations by analytical results in the late time regime of the system.

We close this section by briefly reviewing the weak energy condition (WEC) and the dominant enegy condition (DEC) for a conformal fluid expanding according to the Bjorken flow dynamics, which will be important in our physical analyses in the course of the next sections. These classical energy conditions are commonly postulated in general relativity in order to constrain the content of matter’s energy-momentum tensor used in Einstein’s equations and ensure energy positiveness, even though some quantum effects are known to violate such classical energy conditions [122, 123]. The WEC states that $\langle\hat{T}_{\mu\nu}\rangle t^{\mu}t^{\nu}\geq 0$ for any timelike vector $t^{\mu}$ . It implies the following inequalities for a conformal fluid undergoing Bjorken flow [141, 90, 91],

\displaystyle\hat{\epsilon}(\tau)\geq 0\,;\qquad\{\partial_{\tau}\hat{\epsilon}(\tau)\leq 0,\,\,\tau\partial_{\tau}\ln[\hat{\epsilon}(\tau)]\geq-4\}\,\,\,\Rightarrow\,\,\,-4\leq\left[\frac{\Delta\hat{p}}{\hat{\epsilon}}\right](\tau)\leq 2.

(46)

Therefore, one may violate the WEC by either having a negative energy density and/or by having a normalized pressure anisotropy assuming values outside the aformentioned interval. All the cases with transient, far from equilibrium violations of the WEC analyzed here, and also in Refs. [90, 91], were related to violations on the bounds for the pressure anisotropy in (46), while the energy density was always positive during the time evolution of the medium. We further remark that for a conformal fluid, the strong energy condition (SEC), stating that $\langle\hat{T}_{\mu\nu}\rangle t^{\mu}t^{\nu}\geq-\langle\hat{T}_{\mu}^{\mu}\rangle/2$ , is equivalent to the WEC, since $\langle\hat{T}_{\mu}^{\mu}\rangle=0$ in a conformal system. On the other hand, the DEC states that for any future-directed timelike vector $t^{\mu}$ , $X^{\mu}\equiv-\langle\hat{T}^{\mu\nu}\rangle t_{\nu}$ must also be a future-directed timelike or null vector. This is a sufficient but not a necessary condition to establish causal propagation of matter [121]. It was shown in [90, 91] that for a conformal fluid the DEC implies the following inequalities, which are more restrictive than the WEC (46),

\displaystyle\hat{\epsilon}(\tau)\geq 0\,;\qquad-1\leq\left[\frac{\Delta\hat{p}}{\hat{\epsilon}}\right](\tau)\leq 2.

(47)

IV Results for variations of the initial charge density

In this section, we analyze the time evolution of several different far from equilibrium initial conditions of the Bjorken expanding 1RCBH plasma, where for each profile for the initial subtracted metric anisotropy specified in Eq. (41) and in Table 1, we consider variations of the initial charge density of the medium (25d), $\hat{\rho}(\tau_{0})=\rho_{0}/\tau_{0}$ , while keeping fixed its initial energy density (25a), $\hat{\epsilon}(\tau_{0})=-3a_{2}(\tau_{0})$ (we set $\phi_{s}(\tau_{0},u)=0$ throughout this section). The corresponding results are shown in Figs. 2 — 12. We remark that, due to the presence of several initial data which can be independently varied in the 1RCBH plasma, in the present and in the following sections, we plot the time evolution of many initial conditions in order to explore in details several qualitatively different possibilities for the dynamic evolution of the relevant physical observables of the 1RCBH plasma undergoing Bjorken flow. The several pictures considered illustrate the main general features and observations we summarize in the text.

One notices that by increasing the initial charge density of the 1RCBH plasma (by increasing $\rho_{0}$ ), while keeping its initial energy density fixed, the value of $\mu/T$ in the medium is enhanced, producing the following effects on the physical observables analyzed as functions of the dimensionless time measure $\omega_{\Lambda}$ :

a)

$[\Delta\hat{p}/\hat{\epsilon}](\omega_{\Lambda})$ : the hydrodynamization of the pressure anisotropy of the medium, as measured by its convergence, within small relative tolerances, to the Navier-Stokes regime, is generally delayed as $\mu/T$ increases (in line with what was reported in [85]). The enhancement of $\mu/T$ also leads to different physical possibilities for the maxima and the minima of the pressure anisotropy (not always both extrema are simultaneously present): while for all the initial conditions considered which present a minimum, we observed an increase of its magnitude, for initial conditions presenting a maximum, its magnitude can either increase or decrease depending on the chosen profile for the initial subtracted metric anisotropy $B_{s}(\tau_{0},u)$ ; moreover, in the cases where far from equilibrium transient violations of the DEC (47) and of the WEC (46) are observed, such violations increase with increasing $\mu/T$ , clearly indicating that they become more and more relevant as the system approaches criticality at large charge densities.
b)

$\hat{S}_{\textrm{AH}}(\omega_{\Lambda})/\mathcal{A}\Lambda^{2}$ : the magnitude of this dimensionless ratio always decreases as $\mu/T$ is increased. For different choices of the profile for the initial subtracted metric anisotropy $B_{s}(\tau_{0},u)$ , one may realize different physical possibilities for the time evolution of the entropy of the medium, with the presence or absence of early time (quasi)plateau structures. When single or double plateaus are produced in the time evolution of the entropy, indicating the existence of far from equilibrium time windows with no entropy production in the Bjorken expanding fluid, it is later observed almost or effective violations of DEC from below with $\Delta\hat{p}/\hat{\epsilon}\sim-1$ or $\Delta\hat{p}/\hat{\epsilon}<-1$ , respectively (in line with the observations reported in [90, 91] for the particular case corresponding to the pure thermal SYM plasma with $\mu/T=0$ — here we also observe the aforementioned correlations under much more general situations with $\mu/T\geq 0$ ).¹⁰¹⁰10We remark that, in general, such correlations do not hold in reverse order, i.e. there are evolutions for the 1RCBH plasma with transient violations of DEC which present no far from equilibrium plateaus in the entropy of the medium (see e.g. Fig. 10). In particular, we notice the following general trend regarding the deformation of such plateau structures as $\mu/T$ is increased, leading to progressively lower dips in the pressure anisotropy: first a single plateau is formed, which then becomes more spread out in the $\omega_{\Lambda}$ direction, and next it is progressively deformed into double plateaus, then becoming double quasiplateaus, until the plateau structure is finally lost for sufficiently strong violations of the DEC from below (in such a progression, the lower plateau structure is undone first than the higher plateau). We also observe that it is possible to form a single quasiplateau (i.e. with a time derivative close but not actually zero) in the far from equilibrium entropy, which leads to no posterior violations of the DEC (see e.g. Fig. 7; in particular, for $\rho_{0}=0$ , this has been already seen in [90, 91] for the SYM plasma at $\mu/T=0$ , see the initial condition $\#20$ in those works, corresponding to the curve in cyan with the highest initial entropy). We also remark that we have not observed the final two steps in the aforementioned general trend of progressive deformations of the plateau structures in cases with $\mu/T=0$ (i.e., the progressive deformation of the double plateaus into double quasiplateaus and the final loss of the plateau structure was only observed in association with later, progressively stronger violations of the DEC from below with nonzero values of the initial charge density, which leads to a medium with $\mu/T>0$ — this will be further illustrated in Fig. 28 in appendix A).
c)

$[\hat{\rho}^{4/3}/\hat{\epsilon}](\omega_{\Lambda})$ : the magnitude of this observable always increases as $\mu/T$ is increased, but its qualitative behavior at early times may be very different depending on the chosen profile for the initial subtracted metric anisotropy $B_{s}(\tau_{0},u)$ : while in some cases this observable monotonically decreases in time, in other cases it presents the formation of extrema at early times, with the peculiar feature that the position of these extrema in the axis of the dimensionless time measure $\omega_{\Lambda}$ displays little to almost no variation as $\mu/T$ is enhanced.
d)

$[\langle\hat{O}_{\phi}\rangle/\hat{\epsilon}^{1/2}](\omega_{\Lambda})$ : the effective thermalization of the scalar condensate, as measured by its convergence, within small relative tolerances, to the thermodynamically stable equilibrium, is generally delayed as $\mu/T$ increases, and it only occurs for time scales much larger than the ones observed for the hydrodynamization of the pressure anisotropy. The enhancement of $\mu/T$ always leads to an increase in the magnitude of the scalar condensate, with the formation of at least one maximum as function of the dimensionless time measure $\omega_{\Lambda}$ , with other later extrema possessing smaller magnitudes being also eventually observed for some high values of $\mu/T$ in the medium.

We close this section with an extra technical information regarding the numerical simulations used to obtain the results displayed in Figs. 2 — 12. For almost all the initial conditions considered in those results, it is enough to use $N\sim 20$ collocations points in the radial grid in order to obtain convergent and physically reliable results for all the observables considered. However, specifically for $B_{s}11$ in Table 1 with $a_{2}(\tau_{0})=-7.1$ and $\rho_{0}=1$ , we noted spurious numerical oscillations at early times for the scalar condensate, with such an issue being completely fixed by increasing the number of collocation points to $N\sim 30$ .

V Results for variations of the initial energy density

One notices that by decreasing the initial energy density of the 1RCBH plasma (by decreasing $|a_{2}(\tau_{0})|$ ), while keeping its initial charge density fixed, the value of $\mu/T$ in the medium is enhanced, producing for the physical observables analyzed as functions of the dimensionless time measure $\omega_{\Lambda}$ , generally the same physical possibilities as discussed in section IV. However, for a given profile for the initial subtracted metric anisotropy $B_{s}(\tau_{0},u)$ , considering the enhancement of $\mu/T$ of the medium caused specifically by increasing its initial charge density, or specifically by reducing its initial energy density, it is something that may lead to qualitatively different outcomes for some physical observables. For instance, by comparing Figs. 2 and 13 for $B_{s}1$ in Table 1, one notices that the peak of the normalized pressure anisotropy is reduced (increased) by increasing $\mu/T$ associated to increasing (decreasing) the initial charge (energy) density of the medium.

We close this section with an extra technical information regarding the numerical simulations used to obtain the results displayed in Figs. 13 — 23. For almost all the initial conditions considered in those results, it is enough to use $N\sim 20$ collocations points in the radial grid in order to obtain convergent and physically reliable results for all the observables considered. However, specifically for $B_{s}8$ in Table 1 with $\rho_{0}=0.4$ and $a_{2}(\tau_{0})=-6.4$ , and also for $B_{s}11$ with $\rho_{0}=0.8$ and $a_{2}(\tau_{0})=-6.67$ , we noted spurious numerical oscillations at early times for the scalar condensate, with such issues being completely fixed by increasing the number of collocation points to $N\sim 30$ .

VI Results for variations of the initial dilaton profile

In this section, we analyze the time evolution of some different far from equilibrium initial conditions of the Bjorken expanding 1RCBH plasma, where we consider variations of the initial profile for the subtracted dilaton field specified in Eq. (42) and in Table 2, while keeping fixed all the other initial data. One notices from Eq. (25a) that when the boundary value of the initial profile for the subtracted dilaton field is non-vanishing, $\phi_{2}(\tau_{0})=\phi_{s}(\tau_{0},u=0)\neq 0$ , it contributes to the initial energy density of the medium, $\hat{\epsilon}(\tau_{0})=-3a_{2}(\tau_{0})-\phi_{2}(\tau_{0})^{2}/6$ . The corresponding results are shown in Figs. 24 and 25.

In the case of time evolutions with zero charge density, as displayed in Fig. 24, one typically has relatively small variations of the pressure anisotropy and the entropy at early times in terms of variations of the initial dilaton (keeping fixed the remaining initial data). On the other hand, the relative variations of the scalar condensate before thermalization may be large depending on the chosen profiles for the initial dilaton field. The fact that the scalar condensate may be nonzero when the 1RCBH fluid is far from equilibrium, even for time evolutions with zero charge density, implies that there are far from equilibrium solutions with zero charge density in the 1RCBH model that are different from pure thermal SYM solutions far from equilibrium. Indeed, even though the asymptotic equilibrium state is the same for all initial conditions with $\rho_{0}=0$ , since in those cases $\mu/T=0$ and the equilibrium state in the 1RCBH model only depends on the value of $\mu/T$ , at early times, when $\phi_{s}(\tau_{0},u)\neq 0$ , the time evolutions with zero charge density in the 1RCBH plasma develop transiently nontrivial profiles for the scalar condensate (which vanish when the medium thermalizes). One also notices from Figs. 24 (b) and (d) that the thermalization time associated to the effective equilibration of the scalar condensate is considerably larger than the hydrodynamization time of the pressure anisotropy of the medium.

In the case of time evolutions at finite charge density, as shown in Fig. 25, again there are typically relatively small variations of the pressure anisotropy, entropy, and charge density at early times, while the relative variations of the scalar condensate may be large. The value of $\mu/T$ is typically almost insensitive to variations only on the initial dilaton profile.

Interestingly, and completely different from all the cases considered with $\phi_{s}(\tau_{0},u)=0$ , when the initial dilaton profile is nontrivial, the scalar condensate typically acquires negative values with pronounced dips when the medium is still far from equilibrium.

We close this section by observing in Fig. 24 that, just as it happens for the pressure anisotropy and the scalar condensate, also the normalized entropy converges to a single curve at late times when $\mu/T$ is kept fixed (in the case of Fig. 24, it is zero). In two previous works of ours [90, 91], it was not very clear that the entropy does converge to a single curve (it seemed instead that it was converging to a tiny band of values for different initial conditions at zero chemical potential) because the numerical fits done to obtain the energy scale $\Lambda$ (see below Eq. (43)) were not very precise, since they were calculated at not very late times. In the present work, since we evolved each initial condition for much longer times than in Refs. [90, 91], we were able to obtain precise results for $\Lambda$ , such that the normalized entropy $\hat{S}_{\textrm{AH}}/\mathcal{A}\Lambda^{2}=\tau\hat{s}_{\textrm{AH}}/\Lambda^{2}$ clearly converges to a single curve for different initial conditions evolving into the same value of $\mu/T$ at late times.

VII Conclusions and future perspectives

In this paper, we studied the time evolution of several different far from equilibrium initial states for a hot and dense strongly coupled quantum fluid expanding according to the Bjorken flow dynamics. The corresponding medium, called the 1RCBH plasma, describes a conformal $\mathcal{N}=4$ SYM plasma charged under an Abelian $U(1)$ group of the $SU(4)$ R-symmetry, having a critical point in its phase diagram. We analyzed the time evolution of the Bjorken expanding 1RCBH plasma taking into account the behavior of several physical observables, including for the first time the calculation of the holographic non-equilibrium entropy for this model.

We observed that the value of $\mu/T$ in the medium is enhanced either by increasing its initial charge density, or by decreasing its initial energy density. We found that as $\mu/T$ is enhanced towards its critical value, the hydrodynamization of the pressure anisotropy of the medium, measured by its late time convergence to the corresponding Navier-Stokes regime, is generally delayed, in line with previous work [85]. We have also seen, for the first time, the effective thermalization of the scalar condensate, measured by its late time convergence to the corresponding thermodynamically stable equilibrium, which is also generally delayed as the R-charge density of the medium is increased. The thermalization of the scalar condensate only happens at much later times than the hydrodynamization of the pressure anisotropy, requiring the numerical simulations to run for much longer times than in previous works [85, 90, 91], which in turn demands the numerical code to be implemented using efficient programming languages for numerical calculations.

For some sets of initial data preserving all the classical energy conditions, dynamically driven transient violations of the dominant and the weak energy conditions are observed when the 1RCBH plasma is still far from the hydrodynamic regime. The far from equilibrium violations of the dominant and weak energy conditions get stronger at larger values of $\mu/T$ , indicating that such violations become more and more relevant as the strongly coupled quantum fluid approaches its critical regime. For some of these energy conditions violations, it is observed a clear correlation with a previous formation of different plateau structures in the far from equilibrium entropy of the medium, indicating the presence of transient, early time windows where the Bjorken expanding plasma has zero entropy production even while being far from equilibrium.

More specifically, when single or double plateaus are produced in the time evolution of the far from equilibrium entropy of the medium, it is later observed almost or effective violations of DEC from below with $\Delta\hat{p}/\hat{\epsilon}\sim-1$ or $\Delta\hat{p}/\hat{\epsilon}<-1$ , respectively, in line with the observations in [90, 91] for the pure thermal SYM plasma with $\mu/T=0$ , although here we also observe the aforementioned correlations under much more general situations with $\mu/T\geq 0$ . In general, such correlations do not hold in reverse order, i.e. there are evolutions for the 1RCBH plasma with transient violations of DEC which present no far from equilibrium plateaus in the entropy of the medium.

In particular, we notice the following general trend regarding the deformation of plateau structures as $\mu/T$ is increased, leading to progressively lower dips in the pressure anisotropy to energy density ratio: first a single plateau is formed (near the boundary to DEC violation from below, with $\Delta\hat{p}/\hat{\epsilon}\sim-1$ ), which then becomes more spread out in the time direction, and next it is progressively deformed into double plateaus (with $\Delta\hat{p}/\hat{\epsilon}<-1$ ), then becoming double quasiplateaus, until the plateau structure is finally lost for sufficiently strong violations of the DEC from below (in such a progression, the lower plateau structure is undone first than the higher plateau). We also observe that it is possible to form a single quasiplateau (i.e. with a time derivative close but not actually zero) in the far from equilibrium entropy, which leads to no posterior violations of the DEC. However, we have not observed the final two steps in the aforementioned general trend of progressive deformations of plateau structures in cases with $\mu/T=0$ — i.e., the progressive deformation of the double plateaus into double quasiplateaus and the final loss of the plateau structure was only observed in association with later, progressively stronger violations of the DEC from below (i.e., for $\Delta\hat{p}/\hat{\epsilon}<-1$ ) with nonzero values of the initial charge density, which leads to a medium with $\mu/T>0$ . On the other hand, no particular correlation has been noticed between peculiar features in the far from equilibrium entropy of the medium and violations of DEC and WEC from above (i.e., for $\Delta\hat{p}/\hat{\epsilon}>2$ ).

We also analyzed variations on the initial profile for the subtracted dilaton field. When the initial charge density of the medium is set to zero but the initial dilaton is nontrivial, the early time evolution of the 1RCBH medium is different from pure thermal SYM far from equilibrium states, since in the case of the 1RCBH plasma the scalar condensate develops a nontrivial time dependence.¹¹¹¹11The 1RCBH model features pure thermal SYM as a particular case, which is obtained by setting both the initial charge density and the initial subtracted dilaton profile to zero. At later times, when the 1RCBH fluid approaches equilibrium, the scalar condensate for charge neutral configurations goes to zero and all the other observables also match the pure thermal SYM equilibrium state. This asymptotic time behavior is expected, since the equilibrium state in the 1RCBH model only depends on the value of $\mu/T$ , which vanishes for charge neutral configurations. When the initial charge density and the initial dilaton are both nonzero, by varying only the initial dilaton while keeping the remaining initial data fixed, one typically notices relatively small variations on the pressure anisotropy, entropy and charge density, while the relative variations on the scalar condensate may be large when the system is still far from equilibrium. By varying only the initial dilaton profile, the impact on the value of $\mu/T$ in the medium is typically negligible. Remarkably, in the cases with a nontrivial profile for the initial dilaton field, the scalar condensate typically acquires negative values and develops a pronounced dip when the medium is still far from equilibrium. This is very different from the solutions with zero initial dilaton, for which the scalar condensate has been always observed to be non-negative in the configurations generated in the present work.

Although we have not calculated in the present work the holographic entanglement entropy, whose second order functional derivative is related in Ref. [142] to the Quantum Null Energy Condition (QNEC) [143], which is a local energy condition proposed to hold for any quantum field theory, we analyzed the behavior of the second order proper time derivative of the non-equilibrium entropy associated to the area of the apparent horizon. Due to the noisy numerical data for the proper time evolution of the area of the apparent horizon at finite density (see appendix B.2), we restricted our analysis of the second order derivative of the non-equilibrium entropy to charge neutral solutions corresponding to pure thermal SYM evolutions, which produce significantly less noisy numerical data for the proper time evolution of the area of the apparent horizon.

We looked at several different initial conditions and compared the normalized pressure anisotropy and the normalized logarithmic second order derivative of the non-equilibrium entropy in cases with no violations of the classical energy conditions, and also in cases with violations of DEC and WEC. The main conclusions we obtained are illustrated with the plots in Fig. 26.

We found that when DEC or WEC/DEC are violated from below (by having, respectively, $\Delta\hat{p}/\hat{\varepsilon}<-1$ or $\Delta\hat{p}/\hat{\varepsilon}<-4$ ) and next WEC/DEC are violated from above (by having $\Delta\hat{p}/\hat{\varepsilon}>2$ ), with no further violations afterwards, then the time evolution of (minus) the normalized logarithmic second order derivative of the non-equilibrium entropy is qualitatively similar to the time evolution of the normalized pressure anisotropy, with the former being slightly delayed in time relatively to the latter. If a third violation happens afterwards (generally corresponding to a second violation of DEC from below), then this correlation is lost with $D^{2}_{\tau}\hat{S}_{\textrm{AH}}$ developing more extrema than $\Delta\hat{p}/\hat{\varepsilon}$ (at the boundary of such a third violation, with $\Delta\hat{p}/\hat{\varepsilon}=-1$ , the extra extremum developed by $D^{2}_{\tau}\hat{S}_{\textrm{AH}}$ is an inflection point). When no violations of DEC or WEC are observed, or when just DEC is violated from below, $D^{2}_{\tau}\hat{S}_{\textrm{AH}}$ and $\Delta\hat{p}/\hat{\varepsilon}$ may generally display different behaviors depending on the chosen initial conditions.

In this regard, when DEC or WEC/DEC are violated from below and next WEC/DEC are violated from above, with no further violations afterwards, we found that the curve for $D^{2}_{\tau}\hat{S}_{\textrm{AH}}$ generally anticipates what will happen with the curve for $\Delta\hat{p}/\hat{\varepsilon}$ . Clearly, such a inferred correlation is physically nontrivial and it can be deeply investigated elsewhere.

It would be also interesting to systematically investigate the violation of energy conditions in other far from equilibrium holographic models, as well as possible correlations with early time windows with zero entropy production for the out of equilibrium medium. It may be that the correlations reported here and in [90, 91] are general for strongly coupled quantum fluids, which is then a remarkable qualitative difference with regard to the physical possibilities which may be realized in classical weakly coupled approaches for non-equilibrium media, like kinetic theory, where violations of energy conditions are not observed.

Acknowledgements.

W.B. is grateful to São Paulo Research Foundation (FAPESP) for the received support under grant No. 2022/02503-9. We are greatly grateful to Renato Critelli for developing an earlier version of the symbolic-algebraic and numerical code for the Bjorken flow dynamics of the 1RCBH model [85], from which we heavily based the development of our current code using Mathematica, Fortran and Python for different tasks. We are also grateful to Braudel Maqueira for developing a Docker container which we used to install and run our calculations on a Linux virtual machine.

Appendix A Further physical consistency checks and details

In this appendix, we discuss some extra physical consistency checks of our numerical results and also provide some clear illustrations for the previously discussed behavior regarding the formation of plateau structures in the far from equilibrium entropy of the medium for some selected initial data.

The radial location of the event horizon is determined by the solution of the outgoing radial null geodesics equation subjected to the condition that at asymptotically large times it is given by a zero of the metric coefficient $A(\tau,r)$ ,

\displaystyle\frac{dr_{\textrm{EH}}(\tau)}{d\tau}=A(\tau,r_{\textrm{EH}}(\tau)),\qquad r_{\textrm{EH}}(\tau\to\infty)=r_{\textrm{EH}}^{\textrm{(eq.)}},

(48)

where $r_{\textrm{EH}}^{\textrm{(eq.)}}$ is the largest simple root of equation $A(\tau\to\infty,r)=0$ , corresponding to the radial position of the event horizon in equilibrium. Translating the above conditions to the coordinate $u=1/r$ and making use of the subtracted field $A_{s}(\tau,u)$ defined by Eq. (26a), one obtains the following first order differential equation,

	$\displaystyle\frac{du_{\textrm{EH}}(\tau)}{d\tau}=-u^{2}_{\textrm{EH}}(\tau)A(\tau,u_{\textrm{EH}}(\tau))=-u^{4}_{\textrm{EH}}(\tau)A_{s}(\tau,u_{\textrm{EH}}(\tau))-\frac{1}{2}-u_{\textrm{EH}}(\tau)\lambda(\tau)-u^{2}_{\textrm{EH}}(\tau)\left(\frac{\lambda^{2}(\tau)}{2}-\partial_{\tau}\lambda(\tau)\right),$
	$\displaystyle u_{\textrm{EH}}(\tau\to\infty)=u_{\textrm{EH}}^{\textrm{(eq.)}}(\mu/T),$		(49)

where $u_{\textrm{EH}}^{\textrm{(eq.)}}(\mu/T)$ is the smallest simple root of the following equation,

\displaystyle A(\tau\to\infty,u)=u^{2}A_{s}(\tau\to\infty,u)+\frac{1}{2u^{2}}+\frac{\lambda(\tau\to\infty)}{u}+\frac{\lambda^{2}(\tau\to\infty)}{2}-\partial_{\tau}\lambda\biggr{|}_{\tau\to\infty}=0.

(50)

In practice, we approximate the asymptotic limit $\tau\to\infty$ in Eq. (50) by evaluating it at $\tau=\tau_{\textrm{end}}$ , where the end time of the numerical simulations needs to be sufficiently large so that the dynamical background black hole geometries are near equilibrium at $\tau=\tau_{\textrm{end}}$ . The ‘alternative non-equilibrium entropy density’ associated to the area of the dynamical event horizon can be calculated using Eq. (40) with the substitution $u_{\textrm{AH}}(\tau)\to u_{\textrm{EH}}(\tau)$ . The corresponding results comparing the time evolution of the apparent and the event horizons, and their associated entropy densities, are shown in Fig. 27 for a given far from equilibrium initial condition. One can see that, as expected, the apparent horizon is always behind the event horizon (notice that the boundary is at $u=0$ ) and that they converge at late times. Moreover, at early times, the normalized entropy density $[\hat{s}_{\textrm{AH}}^{4/3}/\hat{\epsilon}](\tau)$ associated to the apparent horizon is always less than the corresponding result involving the event horizon, while both converge to the same stable equilibrium result at late times, as they should. In fact, this is one of the important analytical consistency checks of our numerical simulations in the late time evolution of the system, as previously discussed in section III.

In Fig. 28 (a), we show for $B_{s}5$ in Table 1 at $\mu/T=0$ , a sequence with progressive reductions of the initial energy density of the medium, which progressively lower the dip in the pressure anisotropy leading from no violations to progressively stronger violations of the DEC (47) from below. One notices from Fig. 28 (b) the associated behavior for the entropy of the medium, which progressively goes from an absence of plateaus, to the formation of a single plateau, which is next deformed into double plateaus, with such plateau structures being explicitly confirmed by the calculation of the time derivative of the entropy, as shown in Fig. 28 (c).

On the other hand, we show in Fig. 28 (d) the time derivative of the entropy for $B_{s}11$ in Table 1 at fixed initial energy density, considering progressive enhancements of the initial charge density of the fluid, which lead to a richer variety of steps than in the previous analyzed case at $\mu/T=0$ . Indeed, in Fig. 28 (d) one notices that as $\mu/T$ is progressively increased (leading to stronger violations of the DEC from below, as previously displayed in Fig. 12), the double plateaus are further deformed into double quasiplateaus, until the plateau structure is finally lost for strong enough violations of the DEC at high values of $\mu/T$ .

One also notices that the time derivative of the entropy shown in Fig. 28 (d) at finite $\mu/T$ is clearly more noisy than the results displayed in Fig. 28 (c) at $\mu/T=0$ . Indeed, the results shown in Fig. 28 (d) have been already smoothed out by using a smoothing technique employed to reduce the very strong numerical noise in the original data associated to the area of the apparent horizon evaluated at finite $\mu/T$ , as we discuss next in appendix B.

Appendix B Numerical error analysis

In this appendix, we briefly discuss some details regarding the numerical error and the treatment of the inherent noise in our code.

B.1 Monitoring the convergence and the constraints

We have two constraint equations, denoted here by $\mathcal{C}_{E}$ and $\mathcal{C}_{M}$ , respectively given by Eqs. (15h) and (15a), which are used to globally monitor the time evolution of the initial data. In order to accomplish such a task, we define a root mean square (RMS) norm,¹²¹²12We denote by $u_{\textrm{IR}}$ the infrared radial cutoff used as the end point of the radial integration deep into the bulk [90, 91]. In the present work, we selected different values for $u_{\textrm{IR}}$ (typically between $1.0$ and $1.2$ ) depending on the initial condition considered, such as to have the initial value of the radial location of the apparent horizon $u_{\textrm{AH}}(\tau_{0})$ (calculated with $\lambda(\tau_{0})=0$ ) within the radial grid $u\in[0,u_{\textrm{IR}}]$ .

||\mathcal{C}||_{2}=\sqrt{\frac{1}{2}\int_{0}^{u_{\textrm{IR}}}du\,\mathcal{C}^{2}}\,,

(51)

which is our measure of error for the time evolution of any given initial condition. The integral in Eq. (51) is calculated using a Gauss-Lobatto quadrature off the collocation points. By doing so, we display in Fig. 29, as a typical example, the time evolution of the error of the constraints for some values of the number of radial collocation points for a given initial condition (as discussed in section II.4, we use $\Delta\tau=12\times 10^{-5}$ as the time step size). The constraint $\mathcal{C}_{M}$ , given by Eq. (15a), displays an apparent high and structured error. A careful study of this error led us to efficiently remove the associated noise, which is crucial to the more accurate calculation of the entropy, and particularly its time derivative, which is used to confirm the presence of plateaus with zero entropy production.

B.2 Treatment of numerical noise

As shown in Fig. 30 (a), in the course of time the numerically calculated radial location of the apparent horizon typically displays the same kind of structured noise as observed for the constraint $\mathcal{C}_{M}$ . This behavior is not clearly noticed by eye in the associated entropy,¹³¹³13Although it does become noticeable in the entropy at large times by increasing the number of collocation points $N$ . but its time derivative does manifest it in full swing. An additional test was conducted for the calculation of the apparent horizon. We used an alternative method to the Newton-Raphson algorithm for that purpose, namely, a bisection standard method. The result was numerically the same for a similar numerical tolerance.

For all time evolutions the noise was monitored and, when required by the analysis of the time derivative of the entropy, we treated it within the relevant time window, as displayed e.g. in Fig. 28 (d). In the cases when it was not enough to simply select an appropriate number of radial collocation points in order to diminish the numerical noise increasing with the action of time, we proceed to smooth out the resulting data by using standard procedures and tools. As an extra reference for the particular evolutions with $\mu/T=0$ , we also compared our present results with the outcomes from the numerical code developed for the pure thermal SYM plasma undergoing Bjorken flow in Refs. [90, 91], which is much less noisy. We display in Fig. 30 (b) a comparison for the time derivative of the entropy, where it becomes clear that our best implementation of the smoothing technique corresponds to the moving window procedure. The algorithm is quite simple: we select one window to fit a mean smooth curve which is later differentiated. Depending on the quality of the result — comparing e.g. with the reference pure thermal SYM code [90, 91] — we proceed to move the time windows considered by doing sub-samplings. For that purpose, we implemented a short script using a couple of Python libraries and functions. Particularly, the Scipy interpolating tools for a smooth spline approximation. Without smoothing out the results for the time derivative of the entropy it is not possible to extract useful information from this specific numerical data, as it is clear from Fig. 30 (b).

Appendix C Numerical code’s performance

For this work, we developed a general code which makes use of different programming languages for different tasks.

All the lengthy analytical and symbolic-algebraic manipulations were implemented using Wolfram’s Mathematica and the Riemannian Geometry $\&$ Tensor Calculus (RGTC) library developed by Sotirios Bonanos. Although Mathematica is useful for such tasks, it is extremely inefficient to deal with numerical calculations when compared to other programming languages like e.g. Fortran and C.

Indeed, for the numerical simulations considered in the present work, typically done with $N=21$ radial collocation points, a time step size of $\Delta\tau=12\times 10^{-5}$ , and with the evolutions computed within the long time interval $\tau\in[\tau_{0}=0.2,\tau_{\textrm{end}}=35]$ , it is practically unfeasible to numerically evolve a single initial data using Mathematica. By considering a much shorter end time for the simulations, like $\tau_{\textrm{end}}=7.5$ , as used in [90, 91], Mathematica takes several hours to run a single evolution on an Intel Core i5 1.8 GHz Dual-Core; however, $\tau_{\textrm{end}}=7.5$ is not large enough to allow one to see the effective thermalization of the scalar condensate as analyzed in the present work, which generally requires a much larger value for $\tau_{\textrm{end}}$ . Indeed, as mentioned in section II.4, we used here $\tau_{\textrm{end}}=35$ , which can only be simulated using more efficient programming languages for numerical purposes.

In this regard, we developed an integrated numerical code combining Fortran and Python to automate each run, including all the post-processing and plotting. The performance of this numerical code for each run with the same general configurations as aforementioned, including the calculation of $\Lambda$ and $\mu/T$ in the post-processing, is about $20$ minutes on an Intel Core i5 1.8 GHz Dual-Core.

In order to standardize the runs between different operating systems, we used a Docker container to configure a virtual Linux machine and run the initial data to produce all the results analyze in this work.

References

Arsene et al. [2005] I. Arsene et al. (BRAHMS), Nucl. Phys. A 757, 1 (2005), arXiv:nucl-ex/0410020 .
Adcox et al. [2005] K. Adcox et al. (PHENIX), Nucl. Phys. A 757, 184 (2005), arXiv:nucl-ex/0410003 .
Back et al. [2005] B. B. Back et al. (PHOBOS), Nucl. Phys. A 757, 28 (2005), arXiv:nucl-ex/0410022 .
Adams et al. [2005] J. Adams et al. (STAR), Nucl. Phys. A 757, 102 (2005), arXiv:nucl-ex/0501009 .
Aad et al. [2013] G. Aad et al. (ATLAS), JHEP 11, 183, arXiv:1305.2942 [hep-ex] .
Gyulassy and McLerran [2005] M. Gyulassy and L. McLerran, Nucl. Phys. A 750, 30 (2005), arXiv:nucl-th/0405013 .
Heinz and Snellings [2013] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
Luzum and Petersen [2014] M. Luzum and H. Petersen, J. Phys. G 41, 063102 (2014), arXiv:1312.5503 [nucl-th] .
Shuryak [2017] E. Shuryak, Rev. Mod. Phys. 89, 035001 (2017), arXiv:1412.8393 [hep-ph] .
Florkowski et al. [2018] W. Florkowski, M. P. Heller, and M. Spalinski, Rept. Prog. Phys. 81, 046001 (2018), arXiv:1707.02282 [hep-ph] .
Romatschke and Romatschke [2019] P. Romatschke and U. Romatschke, Relativistic Fluid Dynamics In and Out of Equilibrium, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
Durr et al. [2008] S. Durr et al., Science 322, 1224 (2008), arXiv:0906.3599 [hep-lat] .
Borsanyi et al. [2015] S. Borsanyi et al., Science 347, 1452 (2015), arXiv:1406.4088 [hep-lat] .
Aoki et al. [2006] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature 443, 675 (2006), arXiv:hep-lat/0611014 .
Borsanyi et al. [2010] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo (Wuppertal-Budapest), JHEP 09, 073, arXiv:1005.3508 [hep-lat] .
Borsanyi et al. [2014a] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Phys. Lett. B 730, 99 (2014a), arXiv:1309.5258 [hep-lat] .
Borsanyi et al. [2016] S. Borsanyi et al., Nature 539, 69 (2016), arXiv:1606.07494 [hep-lat] .
Borsanyi et al. [2013] S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 111, 062005 (2013), arXiv:1305.5161 [hep-lat] .
Borsanyi et al. [2014b] S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 113, 052301 (2014b), arXiv:1403.4576 [hep-lat] .
Bazavov et al. [2017] A. Bazavov et al., Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
Borsányi et al. [2021] S. Borsányi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pásztor, C. Ratti, and K. K. Szabó, Phys. Rev. Lett. 126, 232001 (2021), arXiv:2102.06660 [hep-lat] .
Philipsen [2013] O. Philipsen, Prog. Part. Nucl. Phys. 70, 55 (2013), arXiv:1207.5999 [hep-lat] .
Meyer [2011] H. B. Meyer, Eur. Phys. J. A 47, 86 (2011), arXiv:1104.3708 [hep-lat] .
Ratti [2018] C. Ratti, Rept. Prog. Phys. 81, 084301 (2018), arXiv:1804.07810 [hep-lat] .
Alexandru et al. [2022] A. Alexandru, G. Basar, P. F. Bedaque, and N. C. Warrington, Rev. Mod. Phys. 94, 015006 (2022), arXiv:2007.05436 [hep-lat] .
Maldacena [1998] J. M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998), arXiv:hep-th/9711200 .
Gubser et al. [1998] S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Phys. Lett. B 428, 105 (1998), arXiv:hep-th/9802109 .
Witten [1998a] E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998a), arXiv:hep-th/9802150 .
Witten [1998b] E. Witten, Adv. Theor. Math. Phys. 2, 505 (1998b), arXiv:hep-th/9803131 .
Gursoy and Kiritsis [2008] U. Gursoy and E. Kiritsis, JHEP 02, 032, arXiv:0707.1324 [hep-th] .
Gursoy et al. [2008] U. Gursoy, E. Kiritsis, and F. Nitti, JHEP 02, 019, arXiv:0707.1349 [hep-th] .
Gubser and Nellore [2008] S. S. Gubser and A. Nellore, Phys. Rev. D 78, 086007 (2008), arXiv:0804.0434 [hep-th] .
Gubser et al. [2008] S. S. Gubser, A. Nellore, S. S. Pufu, and F. D. Rocha, Phys. Rev. Lett. 101, 131601 (2008), arXiv:0804.1950 [hep-th] .
DeWolfe et al. [2011a] O. DeWolfe, S. S. Gubser, and C. Rosen, Phys. Rev. D 83, 086005 (2011a), arXiv:1012.1864 [hep-th] .
DeWolfe et al. [2011b] O. DeWolfe, S. S. Gubser, and C. Rosen, Phys. Rev. D 84, 126014 (2011b), arXiv:1108.2029 [hep-th] .
Finazzo et al. [2015] S. I. Finazzo, R. Rougemont, H. Marrochio, and J. Noronha, JHEP 02, 051, arXiv:1412.2968 [hep-ph] .
Rougemont et al. [2016a] R. Rougemont, A. Ficnar, S. Finazzo, and J. Noronha, JHEP 04, 102, arXiv:1507.06556 [hep-th] .
Rougemont et al. [2015] R. Rougemont, J. Noronha, and J. Noronha-Hostler, Phys. Rev. Lett. 115, 202301 (2015), arXiv:1507.06972 [hep-ph] .
Finazzo and Rougemont [2016] S. I. Finazzo and R. Rougemont, Phys. Rev. D 93, 034017 (2016), arXiv:1510.03321 [hep-ph] .
Rougemont et al. [2017] R. Rougemont, R. Critelli, J. Noronha-Hostler, J. Noronha, and C. Ratti, Phys. Rev. D 96, 014032 (2017), arXiv:1704.05558 [hep-ph] .
Critelli et al. [2017a] R. Critelli, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, Phys. Rev. D 96, 096026 (2017a), arXiv:1706.00455 [nucl-th] .
Rougemont et al. [2018] R. Rougemont, R. Critelli, and J. Noronha, Phys. Rev. D 98, 034028 (2018), arXiv:1804.00189 [hep-ph] .
Grefa et al. [2021] J. Grefa, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, Phys. Rev. D 104, 034002 (2021), arXiv:2102.12042 [nucl-th] .
Grefa et al. [2022] J. Grefa, M. Hippert, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, (2022), arXiv:2203.00139 [nucl-th] .
Rougemont et al. [2016b] R. Rougemont, R. Critelli, and J. Noronha, Phys. Rev. D 93, 045013 (2016b), arXiv:1505.07894 [hep-th] .
Finazzo et al. [2016] S. I. Finazzo, R. Critelli, R. Rougemont, and J. Noronha, Phys. Rev. D 94, 054020 (2016), [Erratum: Phys.Rev.D 96, 019903 (2017)], arXiv:1605.06061 [hep-ph] .
Critelli et al. [2016] R. Critelli, R. Rougemont, S. I. Finazzo, and J. Noronha, Phys. Rev. D 94, 125019 (2016), arXiv:1606.09484 [hep-ph] .
Rougemont [2020] R. Rougemont, Phys. Rev. D 102, 034009 (2020), arXiv:2002.06725 [hep-ph] .
Järvinen [2022] M. Järvinen, Eur. Phys. J. C 82, 282 (2022), arXiv:2110.08281 [hep-ph] .
Demircik et al. [2021] T. Demircik, C. Ecker, and M. Järvinen, (2021), arXiv:2112.12157 [hep-ph] .
Hoyos et al. [2022] C. Hoyos, N. Jokela, M. Järvinen, J. G. Subils, J. Tarrio, and A. Vuorinen, Phys. Rev. D 105, 066014 (2022), arXiv:2109.12122 [hep-th] .
Knaute et al. [2018] J. Knaute, R. Yaresko, and B. Kämpfer, Phys. Lett. B 778, 419 (2018), arXiv:1702.06731 [hep-ph] .
Cai et al. [2022] R.-G. Cai, S. He, L. Li, and Y.-X. Wang, (2022), arXiv:2201.02004 [hep-th] .
Casalderrey-Solana et al. [2014] J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal, and U. A. Wiedemann, Gauge/String Duality, Hot QCD and Heavy Ion Collisions (Cambridge University Press, 2014) arXiv:1101.0618 [hep-th] .
DeWolfe et al. [2014] O. DeWolfe, S. S. Gubser, C. Rosen, and D. Teaney, Prog. Part. Nucl. Phys. 75, 86 (2014), arXiv:1304.7794 [hep-th] .
Kovtun et al. [2003] P. Kovtun, D. T. Son, and A. O. Starinets, JHEP 10, 064, arXiv:hep-th/0309213 .
Buchel and Liu [2004] A. Buchel and J. T. Liu, Phys. Rev. Lett. 93, 090602 (2004), arXiv:hep-th/0311175 .
Kovtun et al. [2005] P. Kovtun, D. T. Son, and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005), arXiv:hep-th/0405231 .
Arnold et al. [2003] P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 01, 030, arXiv:hep-ph/0209353 .
Kurkela and Lu [2014] A. Kurkela and E. Lu, Phys. Rev. Lett. 113, 182301 (2014), arXiv:1405.6318 [hep-ph] .
Almaalol et al. [2020] D. Almaalol, A. Kurkela, and M. Strickland, Phys. Rev. Lett. 125, 122302 (2020), arXiv:2004.05195 [hep-ph] .
Chesler and Yaffe [2009] P. M. Chesler and L. G. Yaffe, Phys. Rev. Lett. 102, 211601 (2009), arXiv:0812.2053 [hep-th] .
Chesler and Yaffe [2010] P. M. Chesler and L. G. Yaffe, Phys. Rev. D 82, 026006 (2010), arXiv:0906.4426 [hep-th] .
Chesler and Yaffe [2011] P. M. Chesler and L. G. Yaffe, Phys. Rev. Lett. 106, 021601 (2011), arXiv:1011.3562 [hep-th] .
Heller et al. [2012a] M. P. Heller, R. A. Janik, and P. Witaszczyk, Phys. Rev. Lett. 108, 201602 (2012a), arXiv:1103.3452 [hep-th] .
Heller et al. [2012b] M. P. Heller, R. A. Janik, and P. Witaszczyk, Phys. Rev. D 85, 126002 (2012b), arXiv:1203.0755 [hep-th] .
van der Schee [2013] W. van der Schee, Phys. Rev. D 87, 061901 (2013), arXiv:1211.2218 [hep-th] .
Casalderrey-Solana et al. [2013] J. Casalderrey-Solana, M. P. Heller, D. Mateos, and W. van der Schee, Phys. Rev. Lett. 111, 181601 (2013), arXiv:1305.4919 [hep-th] .
Chesler and Yaffe [2014] P. M. Chesler and L. G. Yaffe, JHEP 07, 086, arXiv:1309.1439 [hep-th] .
van der Schee [2014] W. van der Schee, Gravitational collisions and the quark-gluon plasma, Ph.D. thesis, Utrecht U. (2014), arXiv:1407.1849 [hep-th] .
Jankowski et al. [2014] J. Jankowski, G. Plewa, and M. Spalinski, JHEP 12, 105, arXiv:1411.1969 [hep-th] .
Fuini and Yaffe [2015] J. F. Fuini and L. G. Yaffe, JHEP 07, 116, arXiv:1503.07148 [hep-th] .
Chesler [2015] P. M. Chesler, Phys. Rev. Lett. 115, 241602 (2015), arXiv:1506.02209 [hep-th] .
Bellantuono et al. [2015] L. Bellantuono, P. Colangelo, F. De Fazio, and F. Giannuzzi, JHEP 07, 053, arXiv:1503.01977 [hep-ph] .
Pedraza [2014] J. F. Pedraza, Phys. Rev. D 90, 046010 (2014), arXiv:1405.1724 [hep-th] .
DiNunno et al. [2017] B. S. DiNunno, S. Grozdanov, J. F. Pedraza, and S. Young, JHEP 10, 110, arXiv:1707.08812 [hep-th] .
Attems et al. [2017a] M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana, and M. Zilhão, JHEP 01, 026, arXiv:1604.06439 [hep-th] .
Casalderrey-Solana et al. [2016] J. Casalderrey-Solana, D. Mateos, W. van der Schee, and M. Triana, JHEP 09, 108, arXiv:1607.05273 [hep-th] .
Grozdanov and van der Schee [2017] S. Grozdanov and W. van der Schee, Phys. Rev. Lett. 119, 011601 (2017), arXiv:1610.08976 [hep-th] .
Attems et al. [2017b] M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana, and M. Zilhão, JHEP 06, 154, arXiv:1703.09681 [hep-th] .
Romatschke [2018] P. Romatschke, Phys. Rev. Lett. 120, 012301 (2018), arXiv:1704.08699 [hep-th] .
Spaliński [2018] M. Spaliński, Phys. Lett. B 776, 468 (2018), arXiv:1708.01921 [hep-th] .
Critelli et al. [2017b] R. Critelli, R. Rougemont, and J. Noronha, JHEP 12, 029, arXiv:1709.03131 [hep-th] .
Casalderrey-Solana et al. [2018] J. Casalderrey-Solana, N. I. Gushterov, and B. Meiring, JHEP 04, 042, arXiv:1712.02772 [hep-th] .
Critelli et al. [2019] R. Critelli, R. Rougemont, and J. Noronha, Phys. Rev. D 99, 066004 (2019), arXiv:1805.00882 [hep-th] .
Attems et al. [2018] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos, M. Triana, and M. Zilhão, Phys. Rev. Lett. 121, 261601 (2018), arXiv:1807.05175 [hep-th] .
Kurkela et al. [2020] A. Kurkela, W. van der Schee, U. A. Wiedemann, and B. Wu, Phys. Rev. Lett. 124, 102301 (2020), arXiv:1907.08101 [hep-ph] .
Cartwright and Kaminski [2019] C. Cartwright and M. Kaminski, JHEP 09, 072, arXiv:1904.11507 [hep-th] .
Cartwright [2021] C. Cartwright, JHEP 01, 041, arXiv:2003.04325 [hep-th] .
Rougemont et al. [2021] R. Rougemont, J. Noronha, W. Barreto, G. S. Denicol, and T. Dore, Phys. Rev. D 104, 126012 (2021), arXiv:2105.02378 [nucl-th] .
Rougemont et al. [2022] R. Rougemont, W. Barreto, and J. Noronha, Phys. Rev. D 105, 046009 (2022), arXiv:2111.08532 [nucl-th] .
Cartwright et al. [2022a] C. Cartwright, M. Kaminski, and B. Schenke, Phys. Rev. C 105, 034903 (2022a), arXiv:2112.13857 [hep-ph] .
Ecker et al. [2021] C. Ecker, J. Erdmenger, and W. van der Schee, SciPost Phys. 11, 047 (2021), arXiv:2103.10435 [hep-th] .
Ghosh et al. [2021] J. K. Ghosh, S. Grieninger, K. Landsteiner, and S. Morales-Tejera, Phys. Rev. D 104, 046009 (2021), arXiv:2105.05855 [hep-ph] .
Cartwright et al. [2022b] C. Cartwright, M. Kaminski, and M. Knipfer, (2022b), arXiv:2207.02875 [hep-th] .
Bjorken [1983] J. D. Bjorken, Phys. Rev. D 27, 140 (1983).
Gubser [1999] S. S. Gubser, Nucl. Phys. B 551, 667 (1999), arXiv:hep-th/9810225 .
Behrndt et al. [1999] K. Behrndt, M. Cvetic, and W. A. Sabra, Nucl. Phys. B 553, 317 (1999), arXiv:hep-th/9810227 .
Kraus et al. [1999] P. Kraus, F. Larsen, and S. P. Trivedi, JHEP 03, 003, arXiv:hep-th/9811120 .
Cai and Soh [1999] R.-G. Cai and K.-S. Soh, Mod. Phys. Lett. A 14, 1895 (1999), arXiv:hep-th/9812121 .
Cvetic and Gubser [1999a] M. Cvetic and S. S. Gubser, JHEP 04, 024, arXiv:hep-th/9902195 .
Cvetic and Gubser [1999b] M. Cvetic and S. S. Gubser, JHEP 07, 010, arXiv:hep-th/9903132 .
Tian et al. [2014] Y. Tian, X.-N. Wu, and H.-B. Zhang, JHEP 10, 170, arXiv:1407.8273 [hep-th] .
Bekenstein [1973] J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
Hawking [1975] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
Ryu and Takayanagi [2006a] S. Ryu and T. Takayanagi, Phys. Rev. Lett. 96, 181602 (2006a), arXiv:hep-th/0603001 .
Ryu and Takayanagi [2006b] S. Ryu and T. Takayanagi, JHEP 08, 045, arXiv:hep-th/0605073 .
Nishioka et al. [2009] T. Nishioka, S. Ryu, and T. Takayanagi, J. Phys. A 42, 504008 (2009), arXiv:0905.0932 [hep-th] .
Rangamani and Takayanagi [2017] M. Rangamani and T. Takayanagi, Holographic Entanglement Entropy, Vol. 931 (Springer, 2017) arXiv:1609.01287 [hep-th] .
Faulkner et al. [2013] T. Faulkner, A. Lewkowycz, and J. Maldacena, JHEP 11, 074, arXiv:1307.2892 [hep-th] .
Engelhardt and Wall [2015] N. Engelhardt and A. C. Wall, JHEP 01, 073, arXiv:1408.3203 [hep-th] .
Braga and da Rocha [2017] N. R. F. Braga and R. da Rocha, Phys. Lett. B 767, 386 (2017), arXiv:1612.03289 [hep-th] .
Braga and Junqueira [2021] N. R. F. Braga and O. C. Junqueira, Phys. Lett. B 814, 136082 (2021), arXiv:2010.00714 [hep-th] .
Dore et al. [2022] T. Dore, J. M. Karthein, I. Long, D. Mroczek, J. Noronha-Hostler, P. Parotto, C. Ratti, and Y. Yamauchi, (2022), arXiv:2207.04086 [nucl-th] .
Figueras et al. [2009] P. Figueras, V. E. Hubeny, M. Rangamani, and S. F. Ross, JHEP 04, 137, arXiv:0902.4696 [hep-th] .
Friess et al. [2007] J. J. Friess, S. S. Gubser, G. Michalogiorgakis, and S. S. Pufu, JHEP 04, 080, arXiv:hep-th/0611005 .
Buchel et al. [2016] A. Buchel, M. P. Heller, and J. Noronha, Phys. Rev. D 94, 106011 (2016), arXiv:1603.05344 [hep-th] .
Müller et al. [2020] B. Müller, A. Rabenstein, A. Schäfer, S. Waeber, and L. G. Yaffe, Phys. Rev. D 101, 076008 (2020), arXiv:2001.07161 [hep-ph] .
Engelhardt and Wall [2018] N. Engelhardt and A. C. Wall, Phys. Rev. Lett. 121, 211301 (2018), arXiv:1706.02038 [hep-th] .
Hawking and Ellis [1975] S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Space-Time (Cambridge Monographs on Mathematical Physics) (Cambridge University Press, 1975) p. 404.
Wald [2010] R. M. Wald, General Relativity (University of Chicago Press, 2010).
Visser and Barcelo [2000] M. Visser and C. Barcelo, in 3rd International Conference on Particle Physics and the Early Universe (2000) pp. 98–112, arXiv:gr-qc/0001099 .
Costa and Matsas [2022] B. A. Costa and G. E. A. Matsas, Phys. Rev. D 105, 085016 (2022), arXiv:2112.08881 [gr-qc] .
Kurkela et al. [2019] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, and D. Teaney, Phys. Rev. Lett. 122, 122302 (2019), arXiv:1805.01604 [hep-ph] .
Groot [1980] S. R. D. Groot, Relativistic Kinetic Theory. Principles and Applications (Amsterdam, Netherlands: North-holland, 1980).
Baym and Kadanoff [1961] G. Baym and L. P. Kadanoff, Phys. Rev. 124, 287 (1961).
Arnold et al. [2014] P. Arnold, P. Romatschke, and W. van der Schee, JHEP 10, 110, arXiv:1408.2518 [hep-th] .
Finazzo et al. [2017] S. I. Finazzo, R. Rougemont, M. Zaniboni, R. Critelli, and J. Noronha, JHEP 01, 137, arXiv:1610.01519 [hep-th] .
York [1972] J. W. York, Jr., Phys. Rev. Lett. 28, 1082 (1972).
Gibbons and Hawking [1977] G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2752 (1977).
Bianchi et al. [2002] M. Bianchi, D. Z. Freedman, and K. Skenderis, Nucl. Phys. B 631, 159 (2002), arXiv:hep-th/0112119 .
Skenderis [2002] K. Skenderis, Class. Quant. Grav. 19, 5849 (2002), arXiv:hep-th/0209067 .
de Haro et al. [2001] S. de Haro, S. N. Solodukhin, and K. Skenderis, Commun. Math. Phys. 217, 595 (2001), arXiv:hep-th/0002230 .
Ebrahim et al. [2019] H. Ebrahim, M. Asadi, and M. Ali-Akbari, JHEP 09, 023, arXiv:1811.12002 [hep-th] .
Ebrahim and Nafisi [2020] H. Ebrahim and G.-M. Nafisi, Phys. Rev. D 102, 106007 (2020), arXiv:2002.09993 [hep-th] .
Amrahi et al. [2021] B. Amrahi, M. Ali-Akbari, and M. Asadi, Phys. Rev. D 103, 086019 (2021), arXiv:2101.03994 [hep-th] .
Asadi et al. [2021] M. Asadi, H. Soltanpanahi, and F. Taghinavaz, JHEP 05, 287, arXiv:2102.03584 [hep-th] .
Boyd [2001] J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed., Dover Books on Mathematics (Dover Publications, Mineola, NY, 2001).
van der Schee [2015] W. van der Schee, Lectures on numerics in dynamical AdS spacetimes (Personal Website, 2015).
Booth et al. [2009] I. Booth, M. P. Heller, and M. Spalinski, Phys. Rev. D 80, 126013 (2009), arXiv:0910.0748 [hep-th] .
Janik and Peschanski [2006] R. A. Janik and R. B. Peschanski, Phys. Rev. D 73, 045013 (2006), arXiv:hep-th/0512162 .
Ecker et al. [2018] C. Ecker, D. Grumiller, W. van der Schee, and P. Stanzer, Phys. Rev. D 97, 126016 (2018), arXiv:1710.09837 [hep-th] .
Bousso et al. [2016] R. Bousso, Z. Fisher, S. Leichenauer, and A. C. Wall, Phys. Rev. D 93, 064044 (2016), arXiv:1506.02669 [hep-th] .

Holographic entropy production in a Bjorken expanding hot and dense strongly coupled quantum fluid