Stairway to equilibrium entropy

In this section we briefly review the main points regarding the thermodynamics of the 1RCBH model required for the purposes of the present work. Thermodynamic equilibrium is only reached at late times in the evolution of the homogeneous isotropization dynamics starting from far from equilibrium anisotropic initial states.

The 1RCBH model is a conformal field theory with no intrinsic energy scale, therefore, its hot and dense phase diagram is one dimensional, being effectively described in terms of a single control parameter, the dimensionless ratio $\mu/T$, instead of independent $(T,\mu)$. The phase diagram of the 1RCBH model has the peculiar feature of being a limited line segment, $\mu/T\in[0,\pi/\sqrt{2}]$, where for $\mu/T=0$ (and zero scalar condensate) the theory reduces to the purely thermal SYM plasma, while at $\mu/T=\pi/\sqrt{2}$ there is a critical point where second (and higher) order derivatives of the pressure of the medium diverge [54, 62]. However, since the phase diagram of the model ends at the critical point, there is no actual phase transition. For each value of $\mu/T$ there are two different solutions, one of which is thermodynamically unstable and another one which is stable, with the latter being of physical interest for us. In which branch of black hole solutions the system lies within is something controlled by the dimensionless ratio between the black hole charge and the radial position of its equilibrium event horizon, $Q/r_{\textrm{EH}}$ — see Fig. 1 of [62].

Within the thermodynamically stable branch of equilibrium black hole solutions, we shall use in this work two results that will serve as analytical consistency checks of the late time numerical evolution of the scalar condensate and the entropy density, besides serving also to characterize the effective equilibration times of these observables for different initial conditions. As in Ref. [11], we fix here the equilibrium temperature as $T=1/\pi$ and then measure the chemical potential of the medium with respect to this scale, setting $x\equiv\mu/T=\pi\mu$. By solving Eq. (4.10) of [11] for the radial position of the event horizon in equilibrium and then substituting the result into Eq. (4.8) with the aforementioned choice for $T=1/\pi$, one finds the following result for the black hole charge in equilibrium as a function of $\mu/T$,

$\displaystyle Q(x\equiv\mu/T=\pi\mu)=\sqrt{4-\frac{(x/x_{c})^{2}}{2}-\frac{2\,(x/x_{c})^{2}}{1-\sqrt{1-(x/x_{c})^{2}}}},$

(3)

where $x_{c}=(\mu/T)_{c}=\pi/\sqrt{2}=\pi\mu_{c}$ gives the critical chemical potential. Letting $X$ be any physical observable and taking $\hat{X}\equiv\kappa_{5}^{2}X=4\pi^{2}X/N_{c}^{2}$, by means of Eq. (4.24) of [11] one finds that the thermodynamically stable equilibrium value of the normalized scalar condensate is given by,

$\displaystyle\frac{\langle\hat{O}_{\phi}\rangle_{\textrm{eq}}}{T^{2}}=\pi^{2}\,\sqrt{\frac{2}{3}}\,\,Q^{2}(x\equiv\mu/T=\pi\mu),$

(4)

with $Q(x\equiv\mu/T=\pi\mu)$ given by Eq. (3). Moreover, by means of Eq. (4.11) of [11] one obtains the following result for the thermodynamically stable equilibrium value of the normalized entropy density of the medium,

$\displaystyle\frac{\hat{s}_{\textrm{eq}}}{T^{3}}=\frac{\pi^{4}}{4}\left[3-\sqrt{1-\left(\frac{x}{x_{c}}\right)^{2}}\,\right]^{2}\left[1+\sqrt{1-\left(\frac{x}{x_{c}}\right)^{2}}\,\right].$

(5)

The general EMD equations of motion obtained by extremizing the bulk action (1) are [11, 10],

	$\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,$		(6a)
	$\displaystyle\nabla_{\mu}(f(\phi)F^{\mu\nu})=\frac{1}{\sqrt{-g}}\partial_{\mu}(f(\phi)F^{\mu\nu})$	$\displaystyle=0,$		(6b)
	$\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,$		(6c)

while the ansatze for the bulk fields compatible with the symmetries of the homogeneous isotropization dynamics read as follows in generalized infalling Eddington-Finkelstein (EF) coordinates [11],

$$ds^{2}=2dv\left[dr-A(v,r)dv\right]+\Sigma(v,r)^{2}\left[e^{B(v,r)}(dx^{2}+dy^{2})+e^{-2B(v,r)}dz^{2}\right],\,\,\,\,A_{\mu}dx^{\mu}=\Phi(v,r)dv,\,\,\,\,\phi=\phi(v,r),$$

(7)

where $v$ is the EF time defined by the relation,

$$dv=dt+\sqrt{-\frac{g_{rr}}{g_{tt}}}dr.$$

(8)

Since $g_{rr}$ and $g_{tt}$ are the holographic radial and temporal diagonal components of an asymptotically AdS₅ spacetime, as one approaches the boundary of the bulk geometry at $r\rightarrow\infty$, one obtains the time coordinate of the dual quantum field theory living at the boundary, $v\rightarrow t$. Infalling radial null geodesics satisfy $v=\textrm{constant}$ and outgoing radial null geodesics obey $dr/dv=A(v,r)$ [17].³³3Note that $A(v,r)$ here is half the corresponding metric function in the convention adopted in Ref. [17]. There is also a residual diffeomorphism invariance for the metric in (7) associated to the radial shift $r\mapsto r+\lambda(v)$, with $\lambda(v)$ being an arbitrary function of the EF time [24]. Close to the boundary the bulk metric approaches AdS₅, the dilaton field approaches zero and the Maxwell field gives at asymptotically large times the R-charge chemical potential of the strongly coupled quantum fluid living at the boundary. The precise form of the boundary conditions for the bulk functions to be integrated will be specified in subsection II.3.

By substituting the particular ansatze (7) in the general EMD equations (6a) — (6c), one gets the following set of coupled $1+1$ partial differential equations [11, 10],⁴⁴4We employed Eq. (9b) to write down the constraint (9a), taking also into account that, from the definitions of $d_{+}$ and $\mathcal{E}$, one can rewrite $(d_{+}\Phi)^{\prime}=-\partial_{v}\mathcal{E}-A^{\prime}\mathcal{E}-A\mathcal{E}^{\prime}$. In writing down the constraint (9h) we employed the Hamiltonian constraint (9g).

	$\displaystyle\partial_{v}\mathcal{E}+A\mathcal{E}^{\prime}+\left(3\frac{d_{+}\Sigma}{\Sigma}+\frac{\partial_{\phi}f}{f}d_{+}\phi\right)\mathcal{E}$	$\displaystyle=0,$		(9a)
	$\displaystyle\frac{3\Sigma^{\prime}}{\Sigma}+\frac{\partial_{\phi}f}{f}\phi^{\prime}+\frac{\mathcal{E}^{\prime}}{\mathcal{E}}$	$\displaystyle=0,$		(9b)
	$\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,$		(9c)
	$\displaystyle(d_{+}\Sigma)^{\prime}+\frac{2\Sigma^{\prime}}{\Sigma}d_{+}\Sigma+\frac{\Sigma}{12}\left(2V+f\mathcal{E}^{2}\right)$	$\displaystyle=0,$		(9d)
	$\displaystyle\Sigma\,(d_{+}B)^{\prime}+\frac{3}{2}(B^{\prime}d_{+}\Sigma+\Sigma^{\prime}d_{+}B)$	$\displaystyle=0,$		(9e)
	$\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,$		(9f)
	$\displaystyle\Sigma^{\prime\prime}+\frac{\Sigma}{6}\left(3\left(B^{\prime}\right)^{2}+\left(\phi^{\prime}\right)^{2}\right)$	$\displaystyle=0,$		(9g)
	$\displaystyle d_{+}(d_{+}\Sigma)+\frac{\Sigma}{2}(d_{+}B)^{2}-A^{\prime}d_{+}\Sigma+\frac{\Sigma}{6}(d_{+}\phi)^{2}$	$\displaystyle=0,$		(9h)

where $X^{\prime}\equiv\partial_{r}X$ is the directional derivative along infalling radial null geodesics, $d_{+}X\equiv[\partial_{v}+A(v,r)\partial_{r}]X$ is the directional derivative along outgoing radial null geodesics, and $\mathcal{E}\equiv-\Phi^{\prime}$. Eqs. (9a) and (9b) are the nontrivial components of Maxwell’s equation, Eq. (9c) is the dilaton equation, and Eqs. (9d) — (9h) are the nontrivial components of Einstein’s equations. There are five unknown functions, $\{A(v,r),\Sigma(v,r),B(v,r),\phi(v,r),\mathcal{E}(v,r)\}$, to be determined by the five dynamical equations of motion (9b) — (9f), besides three constraints given by Eqs. (9a), (9g), and (9h). Eqs. (9b) — (9g) form a nested set of equations of motion which may be numerically integrated, while the constraints (9a) and (9h) may be employed to check the accuracy of such numerical solutions.

The form of the ultraviolet near-boundary expansions of the bulk fields for the 1RCBH model undergoing homogeneous isotropization dynamics may be written as follows [11],

	$\displaystyle A(v,r)$	$\displaystyle=\frac{\left[r+\lambda(v)\right]^{2}}{2}-\partial_{v}\lambda(v)+\sum_{n=1}^{\infty}\frac{A_{n}(v)}{r^{n}},$		(10a)
	$\displaystyle\Sigma(v,r)$	$\displaystyle=r+\lambda(v)+\sum_{n=1}^{\infty}\frac{\Sigma_{n}(v)}{r^{n}},$		(10b)
	$\displaystyle B(v,r)$	$\displaystyle=\sum_{n=1}^{\infty}\frac{B_{n}(v)}{r^{n}},$		(10c)
	$\displaystyle\phi(v,r)$	$\displaystyle=\sum_{n=2}^{\infty}\frac{\phi_{n}(v)}{r^{n}},$		(10d)
	$\displaystyle\Phi(v,r)$	$\displaystyle=\Phi_{0}(v)+\sum_{n=2}^{\infty}\frac{\Phi_{n}(v)}{r^{n}}.$		(10e)

As detailed discussed in [11], one may derive from the Maxwell’s equation (9b) and from the ultraviolet near-boundary analysis of the bulk fields the following relation,

$$\mathcal{E}(v,r)=2\Phi_{2}\Sigma(v,r)^{-3}e^{2\sqrt{\frac{2}{3}}\phi(v,r)}.$$

(11)

Moreover, the physical observables related to the renormalized one-point functions of the boundary energy-momentum tensor and the charge current operator read as follows [11],

	$\displaystyle\hat{\varepsilon}$	$\displaystyle\equiv\kappa_{5}^{2}\,\varepsilon=\kappa_{5}^{2}\,\langle T_{tt}\rangle=\kappa_{5}^{2}\,\varepsilon_{\textrm{eq}}(x=\mu/T)=\frac{3}{32}\left[3-\sqrt{1-\left(\frac{x}{x_{c}}\right)^{2}}\,\right]^{3}\left[1+\sqrt{1-\left(\frac{x}{x_{c}}\right)^{2}}\,\right],$		(12a)
	$\displaystyle\hat{\rho}$	$\displaystyle\equiv\kappa_{5}^{2}\,\rho=\kappa_{5}^{2}\,\langle J^{t}\rangle=-\Phi_{2}=\kappa_{5}^{2}\,\rho_{\textrm{eq}}(x=\mu/T)=\frac{x}{4\pi}\left[3-\sqrt{1-\left(\frac{x}{x_{c}}\right)^{2}}\,\right]^{2},$		(12b)
	$\displaystyle\Delta\hat{p}(v)$	$\displaystyle\equiv\kappa_{5}^{2}\,(p_{T}-p_{L})=\kappa_{5}^{2}\,\left[\langle T_{xx}\rangle-\langle T_{zz}\rangle\right]=6B_{4}(v),$		(12c)
	$\displaystyle\langle\hat{O}_{\phi}\rangle(v)$	$\displaystyle\equiv\kappa_{5}^{2}\,\langle O_{\phi}\rangle=-\phi_{2}(v),$		(12d)

where $\kappa_{5}^{2}=4\pi^{2}/N_{c}^{2}$ for a strongly coupled SYM plasma (as the 1RCBH model), $\varepsilon$ is the internal energy density and $\rho$ is the R-charge density, both of which are conserved in the homogeneous isotropization dynamics, and $\Delta p$ is the pressure anisotropy of the medium.

As discussed in [11], one may fix the vast majority of the ultraviolet expansion coefficients in Eqs. (10a) — (10e) in terms of just a few undetermined coefficients and their time derivatives. This is accomplished by substituting these expansions in the EMD equations of motion and then solving the obtained algebraic equations order by order in powers of $r$. By considering ultraviolet expansions up to order $n=8$ there remain five undetermined coefficients in such analysis: $\{\Phi_{0}(v),\Phi_{2}(v),A_{2}(v),B_{4}(v),\phi_{2}(v)\}$. The coefficient $\Phi_{0}(v)$ is fixed by the Dirichlet boundary condition for the (nonzero time component of the) Maxwell field,

$\displaystyle\lim_{v\rightarrow\infty}\lim_{r\to\infty}\Phi(v,r)=\lim_{v\rightarrow\infty}\Phi_{0}(v)=\mu,$

(13)

which gives the $U(1)$ R-charge chemical potential at the boundary quantum field theory. The coefficient $\Phi_{2}(v)$ is actually constant and it is related to the conserved R-charge density of the fluid as given in Eq. (12b). The coefficient $A_{2}(v)$ may be fixed in terms of the constant $H\equiv 18A_{2}(v)+\phi_{2}^{2}(v)=-\hat{\varepsilon}/3$ if one knows the coefficient $\phi_{2}(v)$. In fact, the two remaining undetermined ultraviolet coefficients $\{B_{4}(v),\phi_{2}(v)\}$ are dynamical quantities related to the scalar condensate and the pressure anisotropy according to Eqs. (12c) and (12d). The values of these two coefficients at the initial time slice can be freely chosen since they are given by the boundary values of the initial profiles for the metric anisotropy function $B(v,r)$ and the dilaton field $\phi(v,r)$, which are two of the three initial data that must be chosen for the 1RCBH model undergoing isotropization dynamics (the third initial data is the value of $\mu/T$), as we shall discuss in a moment; and once their initial values are specified, their subsequent time evolutions are determined by numerically solving the nested set of partial differential equations previously obtained.

Schematically, the numerical integration of the nested set of $1+1$ partial differential equations of motion describing the homogeneous isotropization dynamics of the 1RCBH model proceeds as follows:

1. a)

   On the hypersurface at the initial time slice $v_{0}$ (which we set to be zero here, $v_{0}=0$), choose the initial profiles for the metric anisotropy function $B(v_{0},r)$ and for the dilaton field $\phi(v_{0},r)$, besides also the value for the dimensionless ratio $\mu/T$ defining the chemical potential of the medium at the boundary, which is associated to the charge of the black hole solution within the bulk;⁵⁵5If one chooses to work with nonzero radial shift function $\lambda(v)$, also its initial value must be chosen, and we set here $\lambda(v_{0}=0)=0$; its time evolution can be obtained by requiring that the radial position of the apparent horizon of the black hole solution remains fixed during the time evolution of the system [24].

2. b)

   Next radially solve the Hamiltonian constraint (9g) to obtain $\Sigma(v_{0},r)$, which at this step fixes the value of $\mathcal{E}(v_{0},r)$ through Eq. (11);

3. c)

   Next radially solve Eq. (9d) to obtain $d_{+}\Sigma(v_{0},r)$;

4. d)

   Next radially solve Eq. (9e) to obtain $d_{+}B(v_{0},r)$;

5. e)

   Next radially solve the dilaton Eq. (9c) to obtain $d_{+}\phi(v_{0},r)$;

6. f)

   Next radially solve Eq. (9f) to obtain $A(v_{0},r)$;

7. g)

   At this step, from the definition of the directional derivative along outgoing radial null geodesics, $d_{+}\equiv\partial_{v}+A(v,r)\partial_{r}$, one has $\{\partial_{v}B(v_{0},r),\partial_{v}\phi(v_{0},r)\}$, which together with the initial profiles chosen for the metric anisotropy and the dilaton field, comprise the set of initial conditions required to evolve $\{B(v_{0},r),\phi(v_{0},r)\}$ to the next time slice $v_{0}+\Delta v$ using discrete numerical integration techniques (here we employ the pseudospectral method [69] to perform the radial integrations, while the time integrations are performed with the 4th order Adams-Bashforth method);

8. h)

   Repeat the previous steps to obtain all the bulk functions in the current time slice and iterate the procedure until reaching any desired end time $v_{\textrm{end}}$ for the numerical simulations (for the calculations performed in the present work, we set $v_{\textrm{end}}=13$, which gives the dimensionless time measure $v_{\textrm{end}}T=13/\pi$).

In Ref. [11] it was considered the formulation of the homogeneous isotropization dynamics of the 1RCBH model with $\lambda(v)=0$. The physics does not depend on the choice of $\lambda(v)$ (since it works like a gauge function), however, depending on the system under consideration, for numerical stability the introduction of this function in the formalism may be required [24]. This was the case for the Bjorken flow dynamics of the 1RCBH model [10]. For the homogeneous isotropization dynamics one does not really need to work with a nonzero radial shift function, nonetheless, for completeness, in the present work we consider a nontrivial $\lambda(v)$. We explicitly checked that the results for all the physical observables at the boundary quantum field theory are the same obtained with vanishing $\lambda(v)$, as it should be.

By considering the ultraviolet near-boundary expansions of the bulk fields with nonzero $\lambda(v)$ and by introducing a new compact radial coordinate $u\equiv 1/r$ suited for numerical integration, we introduce below subtracted bulk fields with the purpose of obtaining radial constants as the boundary values of the subtracted fields to be numerically integrated. We define $u^{p}X_{s}(v,u)\equiv X(v,u)-X_{\textrm{UV}}(v,u)$, where $p\in\mathbb{Z}$ and $X_{\textrm{UV}}(v,u)$ is some ultraviolet truncation of the field $X(v,u)$ such that $X_{s}(v,u=0)$ gives a radial constant. The subtracted fields to be numerically integrated are defined here as follows (note that in the equations below we also provide the boundary conditions for the subtracted fields),⁶⁶6One substitutes the expansions (10a) — (10e) into the equations of motion (9b) — (9g), eliminates all possible coefficients in favor of the others, and then passes from the original radial coordinate $r$ to the new compact radial coordinate $u=1/r$.

	$\displaystyle u^{2}A_{s}(v,u)$	$\displaystyle\equiv A(v,u)-\frac{1}{2u^{2}}-\frac{\lambda(v)}{u}-\frac{\left[\lambda^{2}(v)-2\partial_{v}\lambda(v)\right]}{2},$		(14a)
	$\displaystyle\Rightarrow A_{s}(v,u\to 0)$	$\displaystyle\to\frac{18H-\phi_{2}^{2}(v)}{18}-\frac{\left[\phi_{2}(v)\partial_{v}\phi_{2}(v)+\lambda(v)\left(36H-2\phi_{2}^{2}(v)\right)\right]u}{18}+\mathcal{O}\left(u^{2}\right),$		(14b)
	$\displaystyle\Rightarrow\partial_{u}A_{s}(v,u=0)$	$\displaystyle=-\frac{\phi_{2}(v)\left[\partial_{u}\phi_{s}(v,u=0)+2\lambda(v)\phi_{2}(v)\right]+\lambda(v)\left(36H-2\phi_{2}^{2}(v)\right)}{18},$		(14c)

where we made use of Eq. (17c) in obtaining Eq. (14c) (which is used as an extra boundary condition in the radial integration of $A_{s}$) from Eq. (14b),

	$\displaystyle u^{4}B_{s}(v,u)$	$\displaystyle\equiv B(v,u),$		(15a)
	$\displaystyle\Rightarrow B_{s}(v,u\to 0)$	$\displaystyle\to B_{4}(v)+\left[\partial_{v}B_{4}(v)-4\lambda(v)B_{4}(v)\right]u+\mathcal{O}\left(u^{2}\right),$		(15b)
	$\displaystyle\Rightarrow\partial_{v}B_{4}(v)$	$\displaystyle=\partial_{u}B_{s}(v,u=0)+4\lambda(v)B_{4}(v),$		(15c)

	$\displaystyle u^{2}\Sigma_{s}(v,u)$	$\displaystyle\equiv\Sigma(v,u)-\frac{1}{u}-\lambda(v),$		(16a)
	$\displaystyle\Rightarrow\Sigma_{s}(v,u\to 0)$	$\displaystyle\to-\frac{\phi_{2}^{2}(v)u}{18}-\frac{\left[3\phi_{2}(v)\partial_{v}\phi_{2}(v)-5\lambda(v)\phi_{2}^{2}(v)\right]u^{2}}{30}+\mathcal{O}\left(u^{3}\right),$		(16b)

	$\displaystyle u^{2}\phi_{s}(v,u)$	$\displaystyle\equiv\phi(v,u),$		(17a)
	$\displaystyle\Rightarrow\phi_{s}(v,u\to 0)$	$\displaystyle\to\phi_{2}(v)+\left[\partial_{v}\phi_{2}(v)-2\lambda(v)\phi_{2}(v)\right]u+\mathcal{O}\left(u^{2}\right),$		(17b)
	$\displaystyle\Rightarrow\partial_{v}\phi_{2}(v)$	$\displaystyle=\partial_{u}\phi_{s}(v,u=0)+2\lambda(v)\phi_{2}(v),$		(17c)

$\displaystyle\mathcal{E}_{s}(v,u)\equiv\mathcal{E}(v,u),$

(18)

	$\displaystyle u^{2}\left(d_{+}\Sigma\right)_{s}(v,u)$	$\displaystyle\equiv\left(d_{+}\Sigma\right)(v,u)-\frac{1}{2u^{2}}-\frac{\lambda(v)}{u}-\frac{\lambda^{2}(v)}{2},$		(19a)
	$\displaystyle\Rightarrow\left(d_{+}\Sigma\right)_{s}(v,u\to 0)$	$\displaystyle\to H+\frac{\phi_{2}^{2}(v)}{36}+\mathcal{O}\left(u\right),$		(19b)

	$\displaystyle u^{3}\left(d_{+}B\right)_{s}(v,u)$	$\displaystyle\equiv\left(d_{+}B\right)(v,u),$		(20a)
	$\displaystyle\Rightarrow\left(d_{+}B\right)_{s}(v,u\to 0)$	$\displaystyle\to-2B_{4}(v)+\mathcal{O}\left(u\right),$		(20b)

	$\displaystyle u\left(d_{+}\phi\right)_{s}(v,u)$	$\displaystyle\equiv\left(d_{+}\phi\right)(v,u),$		(21a)
	$\displaystyle\Rightarrow\left(d_{+}\phi\right)_{s}(v,u\to 0)$	$\displaystyle\to-\phi_{2}(v)+\mathcal{O}\left(u^{2}\right).$		(21b)

The equations of motion to be numerically solved as functions of the coordinates $(v,u)$ are obtained by rewriting the original equations of motion in terms of the subtracted fields, whose boundary values were derived above. As detailed discussed in section 5.4 of [11], by discretizing the radial part of these continuum differential equations of motion using the pseudospectral method, one obtains an eigenvalue problem where one needs to invert a diagonal $(N-1)\times(N-1)$ matrix for each of the bulk fields, with $N$ being the number of collocation points of the radial Chebyshev-Gauss-Lobatto grid. These matrices correspond to the homogeneous part of the discretized differential equations of motion evaluated at each radial grid point, with the exception of the boundary grid point. The multiplication of the inverse matrices by the column vectors corresponding to the inhomogeneous part of the equations of motion provides the numerical solutions for the bulk fields. At the boundary grid point one must impose the boundary conditions derived above for the subtracted bulk fields. Then one must join to the $(N-1)$-dimensional eigenvectors obtained as solutions of the aforementioned eigenvalue problem the values of the respective bulk fields determined at the boundary grid point by the associated boundary conditions. In this way one obtains the complete $N$-dimensional eigenvectors providing the numerical solutions for the radial part of the EMD equations of motion at a hypersurface defined at any given time slice (the components of these $N$-dimensional eigenvectors are the values of the bulk fields at each of the $N$ collocation points of the radial grid).

The set of initial data needed to be chosen in order to start the time evolution of the system of partial differential equations of motion is $\{B_{s}(v_{0},u),\phi_{s}(v_{0},u),\mu/T\}$, besides also the value of $\lambda(v_{0})$ if one wants to use a nontrivial radial shift function $\lambda(v)$. As aforementioned, we set here $\lambda(v_{0}=0)=0$. The equation of motion for $\lambda(v)$ comes from Eq. (14a) evaluated at the radial location of the apparent horizon, $u=u_{\textrm{AH}}=1/r_{\textrm{AH}}$, which for any metric of the form shown in Eq. (7) is obtained as the outermost solution of the equation $\left(d_{+}\Sigma\right)(v,r_{\textrm{AH}})=0$ [24], where one imposes that the radial location of the apparent horizon stays fixed during the time evolution of the system by requiring that $\partial_{v}r_{\textrm{AH}}=0\,\Rightarrow\,d_{+}\left[d_{+}\Sigma\right](v,r_{\textrm{AH}})=A(v,r_{\textrm{AH}})\,\partial_{r}\left[d_{+}\Sigma\right](v,r_{\textrm{AH}})$, which leads to the following condition when substituted in the constraint (9h),

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

(22)

Substituting (22) in Eq. (14a) one obtains the equation of motion for the radial shift function,

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

(23)

Eq. (23) evolves in time the initial condition $\lambda(v_{0})$ by shifting the radial coordinate $u$ at each time slice such as to keep $u_{\textrm{AH}}=\textrm{constant}$. The apparent horizon is the outermost trapped null surface inside the event horizon and it (usually) converges to the latter at late times, when the black hole geometry approaches thermodynamic equilibrium. Since the apparent horizon is inside the event horizon, by cutting off the radial integration of the bulk equations of motion at some position inside the apparent horizon, one guarantees that the radial domain of the bulk spacetime causally connected to observers at the boundary is being properly taken into account and, consequently, no physical information is lost in this integration procedure.

In order to evolve in time the initial data $\{B_{s}(v_{0},u),\phi_{s}(v_{0},u)\}$ one needs to obtain the time derivatives $\partial_{v}B_{s}$ and $\partial_{v}\phi_{s}$, what can be done by taking the expressions for $d_{+}B=\partial_{v}B+A\partial_{r}B$ and $d_{+}\phi=\partial_{v}\phi+A\partial_{r}\phi$ rewritten in terms of the subtracted bulk fields and the compact radial coordinate $u=1/r$. One then obtains the following results,

	$\displaystyle\partial_{v}B_{s}(v,u)$	$\displaystyle=\frac{[d_{+}B]_{s}}{u}+\frac{2B_{s}}{u}+\frac{B_{s}^{\prime}}{2}+4u^{3}A_{s}B_{s}+u^{4}A_{s}B_{s}^{\prime}+\left(4B_{s}+uB_{s}^{\prime}\right)\lambda+\left(2uB_{s}+\frac{u^{2}B_{s}^{\prime}}{2}\right)\lambda^{2}$
		$\displaystyle\,\,\,\,\,\,\,-\left(4uB_{s}+u^{2}B_{s}^{\prime}\right)\partial_{v}\lambda,$		(24a)
	$\displaystyle\partial_{v}B_{s}(v,u=0)$	$\displaystyle=\partial_{v}B_{4}(v),$		(24b)

	$\displaystyle\partial_{v}\phi_{s}(v,u)$	$\displaystyle=\frac{[d_{+}\phi]_{s}}{u}+\frac{\phi_{s}}{u}+\frac{\phi_{s}^{\prime}}{2}+2u^{3}A_{s}\phi_{s}+u^{4}A_{s}\phi_{s}^{\prime}+\left(u\phi_{s}^{\prime}+2\phi_{s}\right)\lambda+\left(u\phi_{s}+\frac{u^{2}\phi_{s}^{\prime}}{2}\right)\lambda^{2}$
		$\displaystyle\,\,\,\,\,\,\,-\left(2u\phi_{s}+u^{2}\phi_{s}^{\prime}\right)\partial_{v}\lambda,$		(25a)
	$\displaystyle\partial_{v}\phi_{s}(v,u=0)$	$\displaystyle=\partial_{v}\phi_{2}(v),$		(25b)

with $X_{s}^{\prime}(v,u)\equiv\partial_{u}X_{s}(v,u)$ being calculated at any fixed time slice by applying the pseudospectral finite differentiation matrix [69] to the numerical solution $X_{s}(v,u)$.

Now we derive the holographic formula for the non-equilibrium entropy density of the strongly coupled quantum fluid living at the boundary of the bulk geometry. The famous Bekenstein-Hawking’s equation [70, 71] relates the area of the event horizon of a black hole in equilibrium with its entropy and, through the holographic dictionary, this gives the thermodynamic entropy of the dual quantum field theory also in equilibrium. However, for strongly coupled media out of equilibrium, it has been argued e.g. in Ref. [72] that the corresponding holographic non-equilibrium entropy should be calculated from the area of the apparent horizon instead of the area of the event horizon. This is so because one expects that entropy production is local in time, therefore, associating the out of equilibrium entropy to the area of a dynamic event horizon seems rather unnatural because in such a context the event horizon can only be determined with the knowledge of the entire future evolution of the black hole geometry, consequently, it is a global rather than a local observable.⁷⁷7Ref. [72] also provides a strong counterexample to further justify such a proposal: for the conformal soliton flow [73], corresponding to an ideal fluid, entropy production must be identically zero at all times. The entropy calculated from the area of the apparent horizon in such a case is indeed constant, however, the area of the event horizon diverges showing that it is an inadequate measure of the non-equilibrium entropy of the dual medium at least in this case. In fact, for the conformal soliton flow the system does not approach a stationary state at late times, and the apparent horizon does not converge to the event horizon, differently from what happens in systems where dissipation drives the evolution of the medium, as in the 1RCBH model analyzed in the present work. Also in several other works in the literature [18, 20, 21, 26, 25, 34, 74, 4, 75] the holographic non-equilibrium entropy has been associated to the area of the apparent horizon, and here we adopt the same approach. Note that since the apparent horizon is inside the event horizon and, for the 1RCBH model it does converge to the event horizon at late times, for sufficiently long times the areas of both horizons coincide and provide the same result for the entropy of the dual medium in thermodynamic equilibrium.

The area of the black hole apparent horizon in infalling EF coordinates is given by,

$\displaystyle A_{\textrm{AH}}(v)=\int_{\textrm{AH}}d^{3}x\,\sqrt{|\gamma_{\textrm{AH}}|}\,\biggr{|}_{u=u_{\textrm{AH}}}=\int_{\textrm{AH}}dxdydz\,\sqrt{g_{xx}g_{yy}g_{zz}}\,\biggr{|}_{u=u_{\textrm{AH}}}=|\Sigma(v,u_{\textrm{AH}})|^{3}\,V_{3},$

(26)

where $V_{3}=\int_{\textrm{AH}}dxdydz$ is the volume of the planar horizon parallel to the boundary. Resorting to an analogous expression to the Bekenstein-Hawking’s relation, one sets for the holographic non-equilibrium entropy of the fluid,

$\displaystyle\hat{s}_{\textrm{AH}}(v)\equiv\kappa_{5}^{2}\,s_{\textrm{AH}}(v)=\kappa_{5}^{2}\,\frac{S_{\textrm{AH}}(v)}{V_{3}}=\kappa_{5}^{2}\,\frac{A_{\textrm{AH}}(v)/4G_{5}}{V_{3}}=2\pi|\Sigma(v,u_{\textrm{AH}})|^{3},$

(27)

where, from Eq. (16a),

$\displaystyle\Sigma(v,u_{\textrm{AH}})=u^{2}_{\textrm{AH}}\Sigma_{s}(v,u_{\textrm{AH}})+\frac{1}{u_{\textrm{AH}}}+\lambda(v).$

(28)

Therefore, setting $T=1/\pi$ as before, one has for the normalized non-equilibrium entropy density of the medium,

$\displaystyle\frac{\hat{s}_{\textrm{AH}}(v)}{T^{3}}=2\pi^{4}|\Sigma(v,u_{\textrm{AH}})|^{3},$

(29)

which must converge to the analytical result (5) at late times for any chosen value of $\mu/T$. This in fact happens for all the initial data we analyzed, as will be shown in section III.

Now, by following a similar approach as the one devised in [76], we derive the weak energy condition (WEC) and the dominant energy condition (DEC) for the conformal 1RCBH model undergoing homogeneous isotropization dynamics. These classical energy conditions are usually postulated in general relativity to restrict the content of the energy-momentum tensor of matter employed in Einstein’s equations with the aim of enforcing energy positiveness, although some quantum effects are known to violate these energy conditions [77, 78].

The WEC states that $\langle\hat{T}_{\mu\nu}\rangle t^{\mu}t^{\nu}\geq 0$ for any timelike vector $t^{\mu}$,⁸⁸8We note that for a conformal fluid, like the 1RCBH model, the strong energy condition (SEC), which states 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$ for a conformal medium. while the DEC posits that for any future-directed timelike vector $t^{\mu}$ (i.e. with $t^{v}>0$), $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 [79]).

For the homogeneous isotropization dynamics of the 1RCBH model, one has from Eqs. (4.40) — (4.42) and (4.45) of [11] the following form for the boundary energy-momentum tensor,

$\displaystyle\langle\hat{T}_{\mu\nu}\rangle=\textrm{diag}\left(\hat{\varepsilon},\hat{p}_{T},\hat{p}_{T},\hat{p}_{L}\right)=\textrm{diag}\left(\hat{\varepsilon},\frac{\hat{\varepsilon}+\Delta\hat{p}}{3},\frac{\hat{\varepsilon}+\Delta\hat{p}}{3},\frac{\hat{\varepsilon}-2\Delta\hat{p}}{3}\right).$

(30)

On the other hand, the most general timelike or null vector at the flat boundary may be written as, $t^{\mu}\equiv\left(\sqrt{s^{2}+2\omega^{2}+\xi^{2}},\omega,\omega,\xi\right)$. In fact, this implies that, $t_{\mu}t^{\mu}=-(s^{2}+2\omega^{2}+\xi^{2})+2\omega^{2}+\xi^{2}=-s^{2}\leq 0,\,\forall\,s,\omega,\xi\in\mathbb{R}$.

Now we analyze the WEC,

$\displaystyle 0\leq\langle\hat{T}_{\mu\nu}\rangle t^{\mu}t^{\nu}=\hat{\varepsilon}s^{2}+\left(4\hat{\varepsilon}+\Delta\hat{p}\right)\frac{2\omega^{2}}{3}+\left(2\hat{\varepsilon}-\Delta\hat{p}\right)\frac{2\xi^{2}}{3},$

(31)

and since $\{s,\omega,\xi\}$ are arbitrary real numbers, (31) can only be satisfied $\forall\,s,\omega,\xi\in\mathbb{R}$ if,

$$\boxed{\hat{\varepsilon}\geq 0}\,;\qquad\left.\begin{array}[]{lr}4\hat{\varepsilon}+\Delta\hat{p}&\geq 0\\ 2\hat{\varepsilon}-\Delta\hat{p}&\geq 0\end{array}\right\}\Rightarrow\boxed{-4\leq\frac{\Delta\hat{p}}{\hat{\varepsilon}}\leq 2}\,.$$

(32)

For the DEC one has that,

	$\displaystyle 0$	$\displaystyle\leq X^{v}=-\langle\hat{T}^{vv}\rangle t_{v}=\hat{\varepsilon}\sqrt{s^{2}+2\omega^{2}+\xi^{2}}\Rightarrow\boxed{\hat{\varepsilon}>0}\,,$		(33)
	$\displaystyle 0$	$\displaystyle\geq X_{\mu}X^{\mu}=-s^{2}\left[\hat{\varepsilon}^{2}\right]-2\omega^{2}\left[\hat{\varepsilon}^{2}-\left(\frac{\hat{\varepsilon}+\Delta\hat{p}}{3}\right)^{2}\right]-\xi^{2}\left[\hat{\varepsilon}^{2}-\left(\frac{\hat{\varepsilon}-2\Delta\hat{p}}{3}\right)^{2}\right],$		(34)

and since (34) must be valid $\forall\,s,\omega,\xi\in\mathbb{R}$, it follows that,

	$\displaystyle 0$	$\displaystyle\leq 1-\left(\frac{1+\Delta\hat{p}/\hat{\varepsilon}}{3}\right)^{2}\Rightarrow\left|1+\frac{\Delta\hat{p}}{\hat{\varepsilon}}\right|\leq 3\Rightarrow-4\leq\frac{\Delta\hat{p}}{\hat{\varepsilon}}\leq 2,$		(35)
	$\displaystyle 0$	$\displaystyle\leq 1-\left(\frac{1-2\Delta\hat{p}/\hat{\varepsilon}}{3}\right)^{2}\Rightarrow\left|1-\frac{2\Delta\hat{p}}{\hat{\varepsilon}}\right|\leq 3\Rightarrow\boxed{-1\leq\frac{\Delta\hat{p}}{\hat{\varepsilon}}\leq 2}\,,$		(36)

where (36) is clearly more restrictive than (35).