Paths to equilibrium in non-conformal collisions

TL;DR

This paper investigates the sequence of relaxation processes in non-conformal holographic heavy ion collisions, revealing that different relaxation times can occur in various orders and that the velocity field remains nearly boost-invariant at hydrodynamization.

Contribution

It provides a detailed numerical analysis of multiple relaxation times in non-conformal theories and explores their ordering and implications for heavy ion collision modeling.

Findings

Condensate can remain far from equilibrium after hydrodynamization.

Rapidity distribution width decreases with increasing non-conformality.

Velocity field at hydrodynamization is nearly boost-invariant regardless of non-conformality.

Abstract

We extend our previous analysis of holographic heavy ion collisions in non-conformal theories. We provide a detailed description of our numerical code. We study collisions at different energies in gauge theories with different degrees of non-conformality. We compare four relaxation times: the hydrodynamization time (when hydrodynamics becomes applicable), the EoSization time (when the average pressure approaches its equilibrium value), the isotropization time (when the longitudinal and transverse pressures approach each other) and the condensate relaxation time (when the expectation value of a scalar operator approaches its equilibrium value). We find that these processes can occur in several different orderings. In particular, the condensate can remain far from equilibrium even long after the plasma has hydrodynamized and EoSized. We also explore the rapidity distribution of the energy density at hydrodynamization. This is far from boost-invariant and its width decreases as the non-conformality increases. Nevertheless, the velocity field at hydrodynamization is almost exactly boost-invariant regardless of the non-conformality. This result may be used to constrain the initialization of hydrodynamic fields in heavy ion collisions.

Figures

Equations

$$S=\frac{2}{{\kappa }_5^2}\int d^{5}x\sqrt{-g}\left[\frac{1}{4}\mathcal{R}-\frac{1}{2}{\left(\nabla \varphi \right)}^2-V(\varphi )\right].$$

$$R_{\mu \nu }-\frac{R}{2}{g}_{\mu \nu }=8\pi T_{\mu \nu },$$

$$\square \varphi =\frac{\partial V}{\partial \varphi },$$

$$8\pi T_{\mu \nu }$$

$$L^{2}V(\varphi )=-3-\frac{3}{2}{\varphi }^2-\frac{1}{3}{\varphi }^4+\left(\frac{1}{2{\varphi }_{\mathrm{M}}^4}+\frac{1}{3{\varphi }_{\mathrm{M}}^2}\right){\varphi }^6-\frac{1}{12{\varphi }_{\mathrm{M}}^4}{\varphi }^8,$$

$$L_{\text{IR}}=\sqrt{-\frac{3}{V\left({\varphi }_{\text{M}}\right)}}=\frac{1}{1+\frac{1}{6}{\varphi }_{\text{M}}^2}L.$$

$$m_{\text{IR}}^2=\frac{12}{L^2}\left(1+\frac{1}{9}{\varphi }_{\text{M}}^2\right)=\frac{12}{L_{\text{IR}}^2}\frac{\left(1+\frac{1}{9}{\varphi }_{\text{M}}^2\right)}{{\left(1+\frac{1}{6}{\varphi }_{\text{M}}^2\right)}^2}.$$

$${\Delta }_{\text{IR}}=2+2\sqrt{1+\frac{m_{\text{IR}}^2{L}_{\text{IR}}^2}{4}}=6\left(1+\frac{{\varphi }_{\text{M}}^2}{9}\right){\left(1+\frac{{\varphi }_{\text{M}}^2}{6}\right)}^{-1}.$$

$$d{s}^2=\frac{L^2}{u_{\text{FG}}^2}d{u}_{\text{FG}}^2+e^{2a_{\text{FG}}(u_{\text{FG}})}{\eta }_{\mu \nu }d{x}^{\mu }d{x}^{\nu },$$

$$V(\varphi )=-\frac{4}{3}W{\left(\varphi \right)}^2+\frac{1}{2}{W}^{\prime }{\left(\varphi \right)}^2,$$

$$LW\left(\varphi \right)=-\frac{3}{2}-\frac{{\varphi }^2}{2}+\frac{{\varphi }^4}{4{\varphi }_{\text{M}}^2}.$$

$$u_{\text{FG}}\frac{d{a}_{\text{FG}}}{d{u}_{\text{FG}}}=\frac{2}{3}W,u_{\text{FG}}\frac{d\varphi }{d{u}_{\text{FG}}}=-\frac{\partial W}{\partial \varphi },$$

$$e^{2a_{\text{FG}}}$$

$$\varphi$$

$$\Lambda ={\varphi }_0.$$

$$d{s}^2=\frac{L^2}{u_{\text{FG}}^2}\left(d{u}_{\text{FG}}^2+g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }\right),$$

$$g_{\mu \nu }$$

$$\varphi$$

$$⟨T_{\mu \nu }⟩$$

$$⟨\mathcal{O}⟩$$

$$⟨T_{\mu }^{\mu }⟩=-\Lambda ⟨\mathcal{O}⟩,$$

$$\Big{(}\mathcal{E},J_{\mathcal{E}},P_{x^{i}},\mathcal{V}\Big{)}=\frac{\kappa_{5}^{2}}{2L^{3}}\,\Big{(}-T^{t}_{t},T^{z}_{t},T^{x^{i}}_{x^{i}},\mathcal{O}\Big{)}\,.$$

$$\mathcal{E}-3\bar{P}=\Lambda \mathcal{V},$$

$$\bar{P}=\frac{1}{3}\sum _i{P}_{x^i}$$

$$\bar{P}=P_{\mathrm{eq}}(\mathcal{E}),$$

$$P_{\rm eq}(\mathcal{E})=\frac{1}{3}\Big{[}\mathcal{E}-\Lambda\mathcal{V}_{\rm eq}(\mathcal{E})\Big{]}\,.$$

$$s_R=\frac{{\kappa }_5^2}{2{\pi }^4{L}^3}\frac{s}{T^3},$$

$$\frac{s}{T^3}=\frac{{\pi }^2}{2}{N}^2,$$

$$\frac{s}{T^3}=\frac{2{\pi }^4{L}^3}{{\kappa }_5^2}.$$

$$\frac{\zeta }{\eta }=4{\left(\frac{d\log s}{d{\varphi }_{\text{H}}}\right)}^{-2}.$$

$$d{s}^2=\frac{L^2}{u_{\text{FG}}^2}d{u}_{\text{FG}}^2+f(u_{\text{FG}})h(x_{\pm })d{x}_{\pm }^2+e^{2a_{\text{FG}}(u_{\text{FG}})}\left(-d{x}_{+}d{x}_{-}+d{\mathbf{x}}_{\perp }^2\right)\text{,}$$

$$-f\left[2\left(u_{\text{FG}}^2\frac{{\partial }^2{a}_{\text{FG}}}{\partial u_{\text{FG}}^2}+u_{\text{FG}}\frac{\partial a_{\text{FG}}}{\partial u_{\text{FG}}}\right)+4{\left(-u_{\text{FG}}\frac{\partial a_{\text{FG}}}{\partial u_{\text{FG}}}\right)}^2\right]+u_{\text{FG}}^2\frac{{\partial }^{2}f}{\partial u_{\text{FG}}^2}+u_{\text{FG}}\frac{\partial f}{\partial u_{\text{FG}}}=0.$$

Full text

††preprint: ICCUB-17-009

1 Introduction

“Holographic Heavy Ion Collisions”, namely shockwave collisions in an asymptotically AdS spacetime, have provided interesting insights into the far-from-equilibrium properties of hot, strongly-coupled, non-Abelian plasmas that are potentially relevant for the quark-gluon plasma (QGP) created in heavy ion collision experiments (see e.g. CasalderreySolana:2011us for a review). Until recently, all such holographic studies (see e.g. Chesler:2010bi ; Casalderrey-Solana:2013aba ; Casalderrey-Solana:2013sxa ; Chesler:2015wra ; Chesler:2015bba ; Chesler:2015lsa ; Chesler:2016ceu ) were performed in models dual to conformal field theories (CFTs). One notable lesson of this body of work is that “hydrodynamization”, the process by which the plasma comes to be well described by hydrodynamics, can occur before “isotropization”, the process by which all pressures become approximately equal to one another in the local rest frame.

We have recently begun the study of holographic collisions in non-conformal theories Attems:2016tby ; Attems:2017ezz based on the set of models introduced in Attems:2016ugt .111Some second-order transport coefficients Kleinert:2016nav and the entanglement entropy Rahimi:2016bbv have been computed for these models. One crucial difference between the conformal and the non-conformal cases is that in the latter the equation of state, namely the relation between the energy density and the average pressure, is not fixed by symmetry, and hence it needs not be obeyed out of equilibrium. The relaxation process therefore involves an additional channel, namely the evolution of the energy density and the average pressure towards asymptotic values related by the equation of state. This process was dubbed “EoSization” in Attems:2016tby , and once it has taken place we say that the system has “EoSized”. The main result of Attems:2016tby was that EoSization and hydrodynamization can occur in any order.

The models of Attems:2016ugt are dual to CFTs deformed by a source $\Lambda$ for a dimension-three operator. The source breaks scale invariance explicitly and triggers a non-trivial Renormalization Group (RG) flow. In this paper we will examine the relaxation process by which the expectation value (the condensate) of this scalar operator approaches its equilibrium value. We refer to the time at which this happens as the “condensate relaxation time”, $t_{\textrm{cond}}$. It is particularly interesting to compare this relaxation time to the hydrodynamization, EoSization and isotropization times, $t_{\textrm{hyd}}$, $t_{\textrm{EoS}}$ and $t_{\textrm{iso}}$. The reason is that the latter three times refer to the approach to equilibrium of conserved charges (energy and momentum), whereas the former refers to the relaxation of a non-conserved quantity (the expectation value of the scalar operator). In all the collisions that we have examined we find that isotropization happens last, reinforcing the intuition from conformal collisions that this process is extremely slow. For this reason, in most of the paper we will focus on the other three times and we will come back to $t_{\textrm{iso}}$ in section 6. In contrast, we find that the other three times can occur in several different orderings. In particular, $t_{\textrm{cond}}$ can be much longer than $t_{\textrm{hyd}}$ and $t_{\textrm{EoS}}$. This shows that one-point functions of non-conserved operators can remain far from equilibrium long after a plasma has hydrodynamized and EoSized.

We also examine the physics away from mid-rapidity. For this purpose we compute the rapidity profile of the energy density at hydrodynamization. Just like in the conformal case Casalderrey-Solana:2013aba ; Chesler:2015fpa , this profile is not boost-invariant but Gaussian. The width of this Gaussian decreases as the degree of non-conformality increases. Although the energy profile is determined by far-from-equilibrium physics beyond hydrodynamics, this decrease seems correlated with the bulk viscosity in our models. Indeed, as the non-conformality increases the bulk viscosity grows, which reduces the longitudinal expansion and hence the width of the region where energy is deposited.

A remarkable result of our away-from-mid-rapidity analysis is the fact that, although the energy density profile is far from boost-invariant, the velocity field is almost exactly boost-invariant even for the most non-conformal collisions. For CFTs this was first observed in Chesler:2015fpa . Therefore our result implies that, although the non-conformality has a large effect on the energy density profile at hydrodynamization, it leaves the velocity field essentially unmodified.

This paper is organized as follows. In section 2 we introduce our non-conformal models, along with its thermodynamic and transport properties. In section 3 we describe the numerical procedure used to evolve the corresponding equations, and in section 4 we present tests on the numerical code we have developed to this end. In section 5 we perform a detailed study of shockwave collisions in our models. We conclude with a general discussion in section 6.

2 Setup

2.1 The model

We will consider dynamics in a five-dimensional holographic model consisting of gravity coupled to a scalar field with a non-trivial potential. The action for our Einstein-scalar model is

[TABLE]

The dynamic equations resulting from it read

[TABLE]

where

[TABLE]

and $\kappa_{5}$ is the five-dimensional Newton constant. The potential $V(\phi)$ encodes the details of the dual gauge theory. We choose a simple potential characterised by a single parameter, $\phi_{\textrm{\tiny M}}$, which reads

[TABLE]

where $L$ is a length scale. Note that $V(\phi)$ is negative, possesses a maximum at $\phi=0$ and a minimum at $\phi=\phi_{\textrm{\tiny M}}>0$. A detailed study of this model’s thermodynamics and near-equilibrium properties was presented in Attems:2016tby ; here we will briefly recall the most important points.

The motivation for choosing the potential (5) is that it has three important properties. First, the resulting vacuum solution is asymptotically AdS5 in the UV with radius $L$, since $V(0)=-3/L^{2}$. Second, the second derivative of the potential at $\phi=0$ implies that the scalar field has mass $m^{2}=-3/L^{2}$ therein. This means that, in the UV, this field is dual to an operator in the gauge theory, $\mathcal{O}$, with dimension $\Delta_{\textrm{UV}}=3$. Third, the solution near $\phi=\phi_{\textrm{\tiny M}}$ is again AdS5 with a different radius

[TABLE]

In this region the effective mass of the scalar field differs from its UV value and it is given by

[TABLE]

As a consequence, the operator $\mathcal{O}$ at the IR fixed point has dimension

[TABLE]

To compute the vacuum state of these theories, one needs to first set an ansatz for the solution. In Fefferman-Graham (FG) coordinates, the solution with translation invariance and no horizon can be written in the following form,

[TABLE]

with $a_{\textrm{\tiny FG}}(u_{\textrm{\tiny FG}})$ and $\phi(u_{\textrm{\tiny FG}})$ the non-trivial fields characterising the solution and $u_{\textrm{\tiny FG}}$ the holographic coordinate. The computation of the vacuum state can be simplified when the potential is derived from a super-potential as

[TABLE]

which for the potential selected (5) will be

[TABLE]

In this case, the scalar profile $\phi(u_{\textrm{\tiny FG}})$ and the metric coefficient $a_{\textrm{\tiny FG}}(u_{\textrm{\tiny FG}})$ can be obtained from the equations

[TABLE]

and normalizability boundary conditions. Luckily enough, the equations have an analytic solution for the super-potential chosen,222Note that with respect to our conventions in Attems:2016ugt we have $L\,\phi_{0}^{\textrm{\tiny[here]}}=\phi_{0}^{\textrm{\tiny[there]}}$.

[TABLE]

where $\phi_{0}$ is an arbitrary constant with dimensions of mass that controls the magnitude of the non-normalizable mode of the scalar field. As we will see below, $\phi_{0}$ is equal to the source of the dimension-three operator $\mathcal{O}$ in the dual gauge theory:

[TABLE]

The presence of this source breaks conformal invariance explicitly. Throughout the paper we will use a redundant notation since we will use $\phi_{0}$ when we wish to emphasize the gravitational description and $\Lambda$ when we wish to emphasize the gauge theory scale.

2.2 Gauge theory quantities

Noticing that the small field behaviour of the superpotential (11) is identical to that of the GPPZ flow Girardello:1998pd , we can readily determine the expectation values of the stress tensor and the scalar operator. We begin by expanding the metric and the scalar field in powers of $u_{\textrm{\tiny FG}}$ in the $u_{\textrm{\tiny FG}}\rightarrow 0$ limit. Following Bianchi:2001kw , we write the 5-dimensional metric for asymptotically AdS geometries in generic FG form

[TABLE]

and we write the power expansions of the metric and the scalar field as333Note that with respect to our conventions in Attems:2016ugt we have $\phi^{(2)}_{\textrm{\tiny[here]}}=\Lambda\phi^{(2)}_{\textrm{\tiny[there]}}$.

[TABLE]

The expectation values of the field theory operators are then given by

[TABLE]

As expected, equations (19) and (20) imply the Ward identity for the trace of the stress tensor

[TABLE]

and we adopt a renormalization scheme such that $\left<T_{\mu\nu}\right>=\left<\mathcal{O}\right>=0$ in the vacuum. Henceforth we will omit the expectation value signs and work with the rescaled quantities

[TABLE]

In these variables the Ward identity takes the form

[TABLE]

where

[TABLE]

is the average pressure. Out of equilibrium the average pressure is not determined by the energy density because the scalar expectation value V fluctuates independently. In equilibrium, however, $\mathcal{V}$ is determined by the energy density and the Ward identity becomes the equation of state

[TABLE]

with

[TABLE]

2.3 Thermodynamics and transport

To explore the thermal physics of our model, we search for static black brane solutions of the action (1) following the approach of Gubser:2008ny . Since for these solutions the scalar field is a monotonic function of $u_{\textrm{\tiny FG}}$, we may use it as a coordinate when solving the dynamic equations. The value of $\phi$ at the black brane horizon, $\phi_{\textrm{\tiny H}}$, univocally characterises the black brane solution. Therefore, by imposing the appropriate “horizon” boundary conditions at different $\phi_{\textrm{\tiny H}}$ values one can compute all the equilibrium geometries. Finding the thermodynamics then amounts to finding a family of black brane solutions parametrized by $\phi_{\textrm{\tiny H}}$, and obtaining their Hawking temperatures $T$ and entropy densities $s$. This construction is done is detail in Attems:2016ugt , to where we refer the interested reader.

For our purposes here, it is enough to note that we find a set of values $(\phi_{\textrm{\tiny H}},T,s)$ for each model, i.e. for each $\phi_{\textrm{\tiny M}}$. With these, one can compute all thermodynamic quantities of interest as well as the bulk viscosity $\zeta$. In figure 1 we plot the dimensionless quantity

[TABLE]

as a function of the inverse temperature for two different values of $\phi_{\textrm{\tiny M}}$. Since the theory is conformal both at the UV and at the IR, the high and low temperature behaviour of the entropy density must coincide with that of a relativistic conformal theory and scale as $T^{3}$. In the intermediate region, this scaling is not fulfilled and therefore we can interpret this quantity as a measure of the non-conformality of the gauge theory.

For a relativistic CFT, $s/T^{3}$ is proportional to the number of degrees of freedom in the theory, which for an $SU(N)$ gauge theory with matter in the adjoint representation scales as $N^{2}$. For example, for $\mathcal{N}=4$ SYM

[TABLE]

but the precise coefficient depends on the specific theory. In terms of the parameters of the dual gravity description this quantity becomes

[TABLE]

In our bottom-up setup, the above argument allows us to define the number of degrees of freedom at the fixed points in terms of the effective AdS radius. In particular, the quantity $s_{R}$ should approach 1 at high temperature and $(L_{\textrm{\tiny IR}}/L)^{3}$ at low temperature, which is confirmed by the plots in figure 1.

Another quantity that one can compute from $T$ and $s$ is $P_{\rm eq}(\mathcal{E})$, introduced in (26), also known as the equation of state. This quantity gives another measure of the degree of non-conformality of the gauge theory, and will also be necessary later on for the hydrodynamic estimations. For the representative cases of $\phi_{\textrm{\tiny M}}=2\,,3\,,5\,,20$, this quantity can be seen in figure 2. As expected, both at high and low energies the physics becomes approximately conformal and $P_{\rm eq}$ asymptotes to $\mathcal{E}/3$.

The transport properties of the dual gauge theory plasma also reflect the non-conformal behaviour observed in the equation of state. Due to the isotropy of the plasma, at leading order in gradients transport phenomena are controlled by only two coefficients: the shear viscosity $\eta$ and the bulk viscosity $\zeta$. Because of the universality of the shear viscosity to entropy ratio Kovtun:2004de in all theories with a two-derivative gravity dual, we are ensured that this ratio in our model takes the same value as in the conformal $\mathcal{N}=4$ theory, i.e. $\eta/s=1/4\pi$. On the other hand, the bulk viscosity (which would vanish identically in a CFT) is non-zero in our model. Following Eling:2011ms we determine the bulk viscosity by studying the dependence of the entropy on the value of the scalar field at the horizon

[TABLE]

The temperature dependence of this ratio is shown in figure 3 for different values of $\phi_{\textrm{\tiny M}}$.

2.4 Shockwave metric

In the Fefferman-Graham frame it is possible to find a quasi-analytic solution for a single travelling shockwave on a vacuum background. The metric form will simply correspond to the vacuum metric (9) plus the addition of the term $f(u_{\textrm{\tiny FG}})h(x_{\pm})dx_{\pm}^{2}$:

[TABLE]

where $x_{\pm}=z\pm t$, $z$ is the direction of propagation of the shockwave, and ${\bf x}_{\perp}$ are the perpendicular directions to it. The function $h(x_{\pm})$ is an arbitrary function for the waveform.

The propagation of the shockwave at the speed of light does not alter the vacuum profiles of $a_{\textrm{\tiny FG}}$ and $\phi$, thus the only remaining function to be determined is be $f$. The equation for $f(u_{\textrm{\tiny FG}})$ is a second-order differential equation coming from the Einstein’s equations whose solution can only be obtained numerically:

[TABLE]

From the differential equation one can derive the equivalent integral expression

[TABLE]

An additional difficulty for the computation of the function $f(u_{\textrm{\tiny FG}})$ is that it grows exponentially with $u_{\textrm{\tiny FG}}$. However, inspection of (33) shows that this can be circumvented by computing the redefined function

[TABLE]

which takes values between 0 and 1.

Solving equation (33) order by order, we see that $f(u_{\textrm{\tiny FG}})$ behaves as

[TABLE]

With this expression, the metric (31), and the vacuum profile of the scalar field (14), one obtains from (19)-(20) the dual gauge theory quantities of such a shockwave, namely

[TABLE]

where $P_{L}$ is the longitudinal pressure (along the $z$ direction), and $P_{T}$ the transverse pressure (along the transverse directions ${\bf x}_{\perp}$).

3 Numerical procedure

In this section we set $L=1$ for notational simplicity.

3.1 Evolution equations

We follow the notation of Chesler:2010bi and begin by writing the following 5D metric ansatz in Eddington-Finkelstein (EF) coordinates

[TABLE]

where $A$, $B$, $\Sigma$, and $F$ are functions of the radial coordinate $r$, time $t$ and $z$. The shocks will be propagating along $z$, and $\bm{x}_{\perp}$ denotes the two perpendicular directions $\bm{x}_{\perp}=x_{1},x_{2}$. Note that $t$ is a null time coordinate (usually called $v$ in EF coordinates), i.e. constant-$t$ surfaces are not spacelike but null.

Written in this form, the metric is invariant under the following transformation

[TABLE]

Upon plugging the metric (37) in (2) the resulting system conveniently obeys a particular nested structure, consisting of a sequence of radial ODEs at each $t=\text{const}$ null slice that can be solved in order, see e.g. Winicour:2005ge and references therein.

The equations of motion for our present case are given by

[TABLE]

where, for any function $g$, we define

[TABLE]

Note that these equations are all of the general form

[TABLE]

where $g=\Sigma,~{}F,~{}d_{+}\Sigma,~{}d_{+}B,~{}d_{+}\phi,~{}A,~{}d_{+}F$. These are solved imposing reflecting boundary conditions at the AdS boundary $u=1/r=0$, which take the form

[TABLE]

The subleading coefficient of the scalar field in EF coordinates $\phi_{2}$, introduced in equation (44e), is related to its FG counterpart, $\phi^{(2)}$, through

[TABLE]

The function $\xi(t,z)$ encodes our residual gauge freedom, whereas the functions $a_{4}(t,z)$ and $f_{2}(t,z)$ are constrained to obey

[TABLE]

with $b_{4}$ read off from $B$ through (44b) and both $\phi_{2}$ and $\partial_{t}\phi_{2}$ read off from $\phi$ through (44e).

To solve the resulting system we follow the general approach of Chesler:2010bi ; Chesler:2013lia , with some important differences that we will outline below.

3.2 Expectation values from evolution variables

With the near-boundary behaviours above, together with the Fefferman-Graham expansions (17) and (18), one finds the coordinate transformation relating the fall-off coefficients in each frame. With these, and the expectation values (19) and (20), one can write the expressions for the gauge theory values in terms of our evolution variables $(b_{4},a_{4},f_{2},\phi_{2})$ as

[TABLE]

where $P_{L}$ and $P_{T}$ are the longitudinal and transverse pressures.

3.3 Gauge fixing

We start with the procedure to fix the residual gauge freedom (38). A convenient choice is treating $\xi(t,z)$ as another evolved variable and choosing its evolution equation by requiring that the position of the apparent horizon lie at some constant radial coordinate $r=r_{h}$. We thus want to impose

[TABLE]

at all times, where $\Theta$ is the expansion of outgoing null geodesics for the metric (37). At surfaces $r=\text{const}$, $\Theta$ is given by

[TABLE]

A simple way to impose the conditions (52) numerically is the following

[TABLE]

where $\kappa$ is a positive parameter typically chosen to be 1. The advantage of imposing such a condition is that it is constructed to drive the $\Theta=0$ surface back to $r=r_{h}$ whenever numerical errors accumulate. This turns out to work very well in practice.

Equation (54), when expanded, gives us an equation for $\partial_{t}\xi$ of the form

[TABLE]

to be evaluated at $r=r_{h}$. This is a second-order, linear ODE in the coordinate $z$, which we solve imposing periodicity in $z$.

3.4 Field redefinitions and evolution algorithm

To integrate the resulting system subject to the boundary conditions (44), it is very convenient to introduce $u=1/r$ as our radial coordinate and redefine the evolved variables so that the divergent pieces at $u=0$ are absent.

Motivated by (44), we make the following definitions

[TABLE]

Our equations are then rewritten in terms of the “$g_{1}$” and “$g_{2}$” variables above. $g_{1}$ variables are adapted to the AdS boundary $u=0$. The corresponding resulting equations, however, are extremely long and carry terms with huge powers of the coordinate $u$. Upon trying to solve this system in the whole grid, we were finding that numerical errors would accumulate very early on in the evolution, quickly spoiling the convergence of the solution. We then decided to make use of the system $g_{1}$ only in the vicinity of $u\sim 0$ (grid1, spanning $u\in[0,u_{0}]$)—where a much simpler series expanded version of the aforementioned equations was used—and another grid (grid2, spanning $u\in[u_{0},u_{h}]$) was introduced where the much simpler system of equations $g_{2}$ was used instead.

Our numerical grid thus consists of a double grid in the $u$ direction $u\in[0,u_{0}]\cup[u_{0},u_{h}]$, where $u_{0}$ is typically chosen to be 0.1, and $u_{h}=1/r_{h}$ is typically chosen to be 2 or 3. We integrate the $g_{1}$ equations with boundary conditions given by (44) in grid1; we then read off the integrated values at $u=u_{0}$ and use these as boundary conditions for integrating the $g_{2}$ equations in grid2. Note, however, that we also need to deal with the junction point $u_{0}$ in our $u$-dependent hyperbolic equations $\partial_{t}B(u,t,z)$ and $\partial_{t}\phi(u,t,z)$, given by equation (42). We explain this procedure in appendix A.

We are now in possession of all the necessary equations for the evolution procedure. The evolution algorithm is then as follows:

at any given time $t_{n}$ (which can be the initial time after having performed the transformation (38) that puts the apparent horizon at constant $u$) we know $B(u,t_{n},z)$, $\phi(u,t_{n},z)$, $\xi(t_{n},z)$, $a_{4}(t_{n},z)$ and $f_{2}(t_{n},z)$; 2. 2.

successively solve the elliptic equations (39) (or rather, the corresponding system obtained in terms of the redefined “$g_{1}$” and “$g_{2}$” functions) in the order $\Sigma_{g_{1,2}}$, $F_{g_{1,2}}$, $\dot{\Sigma}_{g_{1,2}}$, $\dot{B}_{g_{1,2}}$, $\dot{\phi}_{g_{1,2}}$, $A_{g_{1,2}}$, which are a sequence of radial ODEs subjected to the boundary conditions (44); 3. 3.

equation (55) is solved to get $\partial_{t}\xi(t_{n},z)$ and afterwards $\partial_{t}B_{g_{1,2}}(t_{n},u,z)$ and $\partial_{t}\phi_{g_{1,2}}(t_{n},u,z)$ can be obtained through equation (42) with (56g) and (56h) (see also appendix A); 4. 4.

obtain $\partial_{t}a_{4}(t_{n},z)$ and $\partial_{t}f_{2}(t_{n},z)$ through (46) and, together with the already obtained $\partial_{t}\xi(t_{n},z)$, $\partial_{t}B_{g_{1,2}}(u,t_{n},z)$, $\partial_{t}\phi_{g_{1,2}}(u,t_{n},z)$, advance all these quantities to time $t_{n+1}$ with a Runge-Kutta procedure or equivalent. 5. 5.

GOTO 1.

3.5 Discretization

Equations (39) are written in a form that decouples the coordinates $u$ and $z$ (the collision axis) and can therefore be solved as ODEs in the $u$ direction for each point in $z$. For this reason, both coordinates can be treated separately. The $z$ direction is discretized on a uniform grid where periodic boundary conditions are imposed, while along the $u$ direction we make use of two grids, grid1 spanning $[0,u_{0}]$ and grid2 spanning $[u_{0},u_{h}]$. Both $u$ grids are Lobatto-Chebyshev grids with $N_{u}+1$ points. The collocation points, given by

[TABLE]

are defined in the range $[-1:+1]$, and can be mapped to our physical grid by

[TABLE]

where $u_{L}$ and $u_{R}$ are the limits of each of the grids.

As the differential equations are solved in $u$ for each $z$ point, the only important operation performed in the $z$ direction are the partial derivatives present in the equations (39). To evaluate these we use a fourth-order accurate (central) finite difference approximation. Also in this direction, we find spurious high-frequency noise common to any finite differencing schemes. In order to remove it we add numerical dissipation to damp these modes. We have therefore implemented the usual Kreiss-Oliger dissipation operator of order 6 Kreiss1973 whereby, after each time step, all our evolved quantities $f\in\{B_{g_{1,2}},\phi_{g_{1,2}},a_{4},f_{2},\xi\}$ are added a term of the form

[TABLE]

where $i$ labels the grid point in the $z$ direction and $\sigma$ is a tuneable dissipation parameter which must be smaller than 1 for stability, and which we have typically fixed to be 0.2. This procedure effectively works as a low-pass filter.

In the radial direction $u$, the use of the Chebyshev-Lobatto grid allow us to use pseudo-spectral collocation methods Boyd2001 . These methods are based in the approximation of our solutions in a basis of known functions, Chebyshev polynomials $T_{n}(X)$ in our case, but, in addition to the spectral basis, we have an additional physical representation and therefore we can perform operations in one basis or the other depending on our needs. Discretization using the pseudo-spectral method consists in the exact imposition of our equations at the collocation points of the Lobatto-Chebyshev grid. Thanks to the trigonometric representation of the Chebyshev polynomials, we can use the Fast Fourier Algorithm (FFT) for changing from one basis to the other. One of the uses of these method is high-accuracy interpolation of any function $f$ to values of $u$ not present in our grid. This can be computed using the standard spectral representation of the function

[TABLE]

where ${\hat{f}}_{k}$ are the coefficients of the spectral basis that are computed from the values of the function in the collocation points through the FFT. The cost of the FFT algorithm scales as $O\left(N_{u}\log N_{u}\right)$, in contrast with the matrix transformation from the physical and spectral representations, which scales as $O\left(N_{u}^{2}\right)$.

As we mentioned previously, cf. equation (43), the radial equations for solving the metric coefficients can be written in the form

[TABLE]

where, again, $g$ represents the metric coefficients previously mentioned. Once our coordinate is discretized, the differential operator becomes an algebraic one acting over the values of the functions in the collocation points taking the form

[TABLE]

where $\mathcal{D}_{uu}$, $\mathcal{D}_{u}$ represent the derivative operator for a Lobatto-Chebyshev grid in the physical representation and $i$, $j$ indices in the $u$ coordinate. We now construct the operator defined inside the brackets and then invert it to solve the function $g$. Boundary conditions are imposed by replacing full rows in this operator by the values we need to fix. In the general case, for a second order operator we replace the lines $j=0$, $j=N$ by the value of the function and its derivative at $u=0$ in the case of grid1 and at $u=u_{0}$ in the case of grid2. At grid1, we obtain the boundary conditions from (44); at grid2 these are read off from the obtained values at grid1.

Another useful feature of the spectral methods is the possibility of filtering. As we did with the dissipation in the direction $z$, we can damp high order modes but in this case directly in the spectral representation. After each time step, we apply an exponential filter to the spectral coefficients of our $u$-dependent evolved quantities $\hat{f}\in\{\hat{B}_{g_{1,2}},\hat{\phi}_{g_{1,2}}\}$. The complete scheme is

[TABLE]

where $\alpha$ and $\gamma$ are tuneable parameters which we typically fix to $\alpha=36.0437$, $\gamma=8$. This effectively dampens the coefficients of the higher-order Chebyshev polynomials.

3.6 Initial data

Our chosen formulation of Einstein’s equations, known as the characteristic formulation, allows one to specify the initial data needed for an evolution through freely setting the functions $B(u,z)$, $\phi(u,z)$, $\xi(z)$, $a_{4}(z)$ and $f_{2}(z)$. For our intended applications, we wish to have initial data resembling an ultra-relativistic projectile, such as the shockwave metric in AdS. The starting point to construct such initial data is thus the shockwave metric in FG coordinates (31). Once the function $f(u_{\textrm{\tiny FG}})$ therein is computed, one can proceed to transform the metric to the EF frame (37) in which the numerical integration is performed. Owing to the fact that both the FG and the EF metrics have an explicit Killing vector, one can use the following ansatz for the coordinate transformation between the two frames

[TABLE]

for a left-moving shock Chesler:2013lia . The differential equations for the transformation functions $\lambda_{1}(u,z)$, $\lambda_{2}(u,z)$, and $\lambda_{3}(u,z)$ are obtained by simply taking the slots $g_{uu}$, $g_{ut}$, and $g_{uz}$ from the equation

[TABLE]

Equivalently, one might use the fact that the EF coordinate $u$ is a non-affine parameter for ingoing null geodesics

[TABLE]

where $k^{\mu}(u)$ is the parametrized geodesic, and $F(u)=\frac{-2}{25u}$ is a non-affinity function set to meet the desired EF frame with $g_{tr}=1$. The geodesic equation has the advantage of being explicitly dependent on $t+z$ and therefore its solution reduces to a set of ODEs parametrised by the boundary point $z$ for $t=0$. We thus write our initial data for a left-moving shock as follows:

[TABLE]

where $u_{\textrm{\tiny FG}}$, $\lambda_{1,2,3}$ are functions of $u$ and $z$ obtained from (62). Recall that the function $h$ enters the metric (31) and specifies the energy density, the longitudinal pressure and the energy flux in the initial state according to (36). The choice (65) corresponds to a Gaussian profile with width $\omega$ and height $\mu^{3}/\omega\sqrt{2\pi}$. $\mathcal{E}$ is the energy density per unit volume of the boundary field theory and $\mu^{3}$ is the energy density per unit transverse area. As usual Chesler:2010bi ; Casalderrey-Solana:2013aba ; Casalderrey-Solana:2013sxa ; Casalderrey-Solana:2016xfq ; Attems:2016tby we have added a “regulator” $\mathcal{E}_{0}$, namely a background thermal bath with energy density much smaller than all other scales of interest, in order to avoid the large gradients that develop in the deep IR. Given $\mathcal{E}_{0}$, we know the solution in the absence of shocks, which has $B=0$. In particular, we know the subleading coefficient $\phi_{2}$ of the scalar operator as a function of $\mathcal{E}_{0}$, and this is the value that features in equation (68). An important point is that the $z$-independent equilibrium value $\phi(u)$ in FG coordinates is only known numerically. Equation (68) is a good approximation to this numerical solution. The advantage of having an analytic approximation is that, in order to locate the apparent horizon in EF coordinates in the presence of the shocks, it is necessary to know the value of the scalar field in FG coordinates slightly beyond the position of the horizon in those coordinates. Equation (68) provides a good approximation to this value simply by declaring that it applies beyond the horizon. We have verified that the analytic form (68) quickly relaxes upon time evolution and therefore that this way of initializing our code has no effect whatsoever on the collision dynamics. We choose the initial value for $a_{4}$ by comparing (47) and (66). Finally, the function $\xi(z)$ is initialized by imposing that the apparent horizon lie at a constant value of the $u$ coordinate.

4 Code tests

We implement the above construction in a standalone C code, where we use the GNU Scientific Library GSL to solve the linear system (39), the FFTW3 library FFTW05 for FFTs, and use a fourth-order Adams-Bashforth method to integrate the functions $B(u,z)$, $\phi(u,z)$, $a_{4}(z)$, $f_{2}(z)$ and $\xi(z)$ forward in time, using the procedure outlined in section 3.4. The code is trivially parallelized with OpenMP. The resulting code is quite fast, being able to evolve a configuration with

[TABLE]

with 12 + 48 $u$-points and $\phi_{0}\,\Delta z=1/20$ (400 $z$-points) from $t=0$ to $\phi_{0}\,t=1$ in 3 minutes using two cores Intel i7-4820K CPU @ 3.70GHz.

4.1 Quasi-normal modes

In order to test the code and our numerical implementation we have recovered some quasi-normal frequencies reported in Attems:2016ugt . For these tests, we evolved a $\phi_{\rm M}=10$ $z$-independent configuration where the energy density was set to $\mathcal{E}/\phi_{0}^{4}=0.379686$. $a_{4}$ and $\phi_{2}$ were initialised to their corresponding equilibrium values, whereas $B$ and $\phi$ were set to

[TABLE]

Since this configuration is not in equilibrium, $b_{4}$ and $\phi_{2}$ will oscillate and relax, allowing us to compute the quasi-normal modes (QNM) of the system.

Gravitational set-ups containing a single scalar field will typically show two scalar, independent, gauge invariant types of perturbations, each one with its own tower of modes. Hence, the system will have two independent channels to relax to equilibrium. In the model studied in this work, the two channels control independently the fluctuations of the anisotropy and the trace of the stress-energy tensor of the dual plasma respectively. Since $b_{4}$ only contributes to the anisotropy and $\phi_{2}$ only to the trace, their fluctuations will be governed by different towers of modes. Therefore, the frequencies extracted from $b_{4}$ should match the anisotropy tower frequencies’ and the ones from $\phi_{2}$ should match the trace, or “bulk”, tower Attems:2016ugt .

In figure 4, we have fitted numerical data with damped sinusoidals of the form

[TABLE]

In order to recover the frequencies we employed the following strategy. First, we look for the lowest frequency mode. For that, we set $A_{2}=0$ in equation (73) and fit this function to our numerical data. We perform a series of fits to the data, each fit starting at a later time: we start by using the whole signal, then use only the portion $\phi_{0}t\in[1,\infty[$ (say) of the signal, then only the portion $\phi_{0}t\in[2,\infty[$ and so on. The frequencies $\omega^{(1)}$ thus obtained in each fit eventually converge to some value, the longest lived mode, which we are able to isolate through this process. We then fix the $C$, $A_{1}$, $\omega^{(1)}_{r}$, $\omega^{(1)}_{i}$, $\varphi_{1}$ fitting parameters obtained; the corresponding fit is labelled “fit1” in figure 4. Having fixed these parameters we then repeat the process using equation (73), where this time we only allow for the $A_{2}$, $\omega^{(2)}_{r}$, $\omega^{(2)}_{i}$, $\varphi_{2}$ parameters to vary. We thus obtain the frequencies $\omega^{(2)}$; the final resulting fit is labelled “fit2” in figure 4.

The results obtained with this procedure are displayed in figure 4. For the non-conformal mode (top panel) we have obtained

[TABLE]

which are to be compared with

[TABLE]

obtained in Attems:2016ugt . For the anisotropic mode (bottom panel), we have obtained

[TABLE]

which are to be compared with

[TABLE]

obtained in Attems:2016ugt .

We emphasise that the numbers from Attems:2016ugt and those of this section were obtained in a completely independent way, and the excellent agreement between them (of up to $0.004\%$ for the lowest frequency) validates both the code presented herein as well as the method of Attems:2016ugt .

4.2 Convergence analysis

Numerical simulations using finite differencing techniques typically approximate the continuum solution of the problem with an error that depends polynomially on the grid spacing $h$,

[TABLE]

Different numerical implementations will give different convergence orders $n$. In our case, since we make use of fourth-order finite difference stencils, we expect to see $n=4$. One simple way to check for consistency of a code is evolving the same configuration with coarse, medium and fine resolution, $h_{c}$, $h_{m}$ and $h_{f}$. One can then compute a convergence factor given by

[TABLE]

where $f_{h}$ is a chosen evolved variable obtained with numerical resolution $h$. Since in the radial direction we make use of pseudo-spectral methods, our error will be dominated by the resolution used in the $z$ direction, to which the grid spacing $h$ alludes to. For the analysis done in this section we therefore always make use of the same resolution in the radial direction.

We show in figure 5 the convergence properties of our code obtained for a “typical” shockwave collision with physical parameters (cf. section 3.6)

[TABLE]

This configuration was evolved with $\phi_{0}\,h_{c}=40$, $\phi_{0}\,h_{m}=60$ and $\phi_{0}\,h_{f}=80$; the expected convergence factor expected for fourth order convergence would therefore be $Q\approx 5.94$. Plotted in the figure are the results obtained for the energy density at $\phi_{0}\,t=20$, where the differences $|f_{h_{m}}-f_{h_{f}}|$ have been amplified by $Q=5.94$. The results show fourth-order convergence. We have further verified that the values obtained for our medium resolution run are within $\sim 0.4\%$ of the fourth-order Richardson-extrapolated ones, giving us an estimate of the error incurred in the simulation.

5 Non-conformal collisions

5.1 Time evolution

Using the numerical procedure described in section 3 we are now ready to explore and characterise shockwave collisions in different non-conformal theories. As in the analysis of conformal shockwave collisions in Chesler:2010bi ; Casalderrey-Solana:2013aba ; Casalderrey-Solana:2016xfq , we employ Gaussian energy density profiles in the longitudinal direction, equation (65). We choose $t=0$ as the time at which the two incoming shocks would exactly overlap in the absence of interactions.

In a CFT scale invariance guarantees that the physics can only depend on the dimensionless product of the transverse energy scale and the width of the shock, $\mu\omega$. In contrast, in a non-conformal theory with an intrinsic scale $\Lambda$ the physics will also depend on the ratio $\mu/\Lambda$. We will see that, by varying this last ratio for a fixed shock profile ($\mu\omega=\mathrm{const}$), we can study the collision dynamics from low to high energies. Indeed, our model is specified by the value of the parameter $\phi_{\textrm{\tiny M}}$, which controls the degree of non-conformality of the dual gauge theory. For any value of $\phi_{\textrm{\tiny M}}$, when $\mu$ and $\Lambda$ are of the same order, the formation and relaxation of the plasma happens in the most non-conformal region, while for large $\mu/\Lambda$, the early time evolution is approximately as that in a CFT. We will consider two different values $\mu\omega\simeq 0.30$ and $\mu\omega\simeq 0.12$, corresponding to what were dubbed “1/2-shocks” and “1/4-shocks” in Casalderrey-Solana:2013aba .

Using equations (22) and (47)-(50) we extract the stress tensor from our numerical evolution for different values of $\mu/\Lambda$. Following the standard Landau matching procedure, we define the local energy density and a velocity field by determining the time-like eigenvalue of the stress tensor.

[TABLE]

As a consequence of $z$-reflection symmetry, at $z=0$ the local and collision frames coincide and the local energy density is given by $\mathcal{E}$. Given the energy density, we can assign a value of the transport coefficients $\zeta(\mathcal{E})$ and $\eta(\mathcal{E})$ to each spacetime point after the collisions. Since $\zeta$ vanishes in a conformal theory, we can use the assigned ratio $\zeta/\eta$ as a measure of non-conformality.

In figure 6 we plot444For this plot and for all the results in this section we use several values of the energy density $\mathcal{E}_{0}$, ranging between $\mathcal{E}_{0}=\frac{\mu^{3}}{\sqrt{2\pi}\omega}(0.005,0.02)$, and we check that the effects of this regulator are small and in the linear regime. Furthermore, we extrapolate all physical results to $\mathcal{E}_{0}=0$ checking that first and second order extrapolations converge to the same value. the time dependence of this ratio of viscosities at $z=0$ for several representative values of $\mu/\Lambda$. For large $\mu/\Lambda$ values, the energy deposited by the collision in this central region is also large and the system is close to conformality. As the system expands, the energy decreases proving regions of larger and larger $\zeta/\eta$. For smaller values of $\mu/\Lambda$ the system stays in the non-conformal region from an earlier time. Although after a collision the energy at $z=0$ is continuously decreasing, we have not extended our simulation long enough to recover conformal dynamics at late time, as we would expect to happen from the IR behaviour of our model.

The assigned values of transport coefficients also control the dynamics of the stress tensor soon after the collision. In other words, hydrodynamics becomes applicable. At first order in the gradient expansion the hydrodynamic stress tensor may be expressed as

[TABLE]

where $g$ is the Minkowski metric and $P_{\rm eq}\left(\mathcal{E}_{\rm loc}\right)$ is the equilibrium pressure, $\sigma_{\mu\nu}$ and $\Pi$ are the shear and bulk tensors constructed from gradients of the velocity field, and $\Delta_{\mu\nu}$ is the projector on the fluid rest frame. As we will see, at sufficiently late times this expression approximates well the evolution of the full stress tensor.

To illustrate the non-conformal nature of the collision dynamics, in figure 7 we show the time evolution of the transverse (top) and the longitudinal (bottom) pressures at $z=0$ for a collision with $\mu/\Lambda=0.77$ and $\mu\omega=0.30$ in the $\phi_{\textrm{\tiny M}}=20$ model (black solid lines). We compare these evolutions with the first-order hydrodynamic prediction (84) turning on sequentially the two non-conformal properties in the hydrodynamic approximation, namely the non-conformal EoS and the non-zero $\zeta$. As represented by the solid red curve $P_{L,T}^{\rm hyd}(\zeta=0,\mbox{EoS})$, we see that assuming a conformal EoS and $\zeta=0$ fails to reproduce the time evolution. The inclusion of the correct equation of state, represented by the dashed red curve $P_{\{L,T\}}^{\rm hyd}(\zeta=0)$, brings the hydrodynamic prediction closer to the true evolution. Finally, the inclusion of bulk viscosity, represented by the blue solid curve $P_{L,T}^{\rm hyd}$, increases the convergence of the first order hydrodynamic prediction to the evolution of the pressures. As stated, the post collision dynamics in this regime is intrinsically non-conformal.

5.2 Hydrodynamization and EoSization

Inspection of figure 7 indicates that hydrodynamics provides a good description of the evolution of the stress tensor even when the difference between the longitudinal and the transverse pressures is large, which signals the presence of large gradient corrections. This fact, first noted for conformal systems in Heller:2011ju ; Chesler:2010bi , led to the concept of “hydrodynamization”, i.e. the process by which hydrodynamics comes to describe the dynamics of an interacting system, even if the system is far from local thermal equilibrium. In this section we systematically explore this process for different collision energies in four different non-conformal theories, parametrized by four values of the parameter $\phi_{\textrm{\tiny M}}$.

As is common in the literature, we define the hydrodynamization time as the time beyond which both pressures are described by hydrodynamics within a given accuracy. However, in contrast to conformal dynamics, where the tracelessness of the stress tensor fixes the relation between the longitudinal and the transverse pressure, in an non-conformal theory the evolution of these two quantities is unconstrained. For this reason, we introduce independent hydrodynamization criteria for each of the pressure components. As in Attems:2016tby we define the hydrodynamization time $t_{\textrm{hyd}}$ as the time beyond which the difference between the true pressures and the first-order hydrodynamics prediction is less than $10\%$,

[TABLE]

Note that we have used the average pressure $\bar{P}$ as the characteristic scale of the stress tensor, which agrees with the criterium used in Casalderrey-Solana:2013aba in the conformal case.

As noted in Attems:2016tby hydrodynamization is only sensitive to particular combinations of the shear and bulk contributions to the pressure, which we denote $P_{\eta}$ and $P_{\zeta}$. Since the shear tensor is traceless and the bulk tensor is diagonal in the local rest frame, we can write

[TABLE]

Form this decomposition it is clear that the average hydrodynamic pressure is only sensitive to bulk gradients. Furthermore, as discussed in section 2.2, the equation of state relates the average pressure of the system to the energy density via the equilibrium value of $\mathcal{V}$, the thermal expectation value of the dimension-three operator which deforms the dual gauge theory. For this reason in Attems:2016tby we introduced the EoSization time $t_{\rm EoS}$ as the time beyond which the average pressure agrees with the equilibrium pressure with a $10\%$ accuracy,

[TABLE]

In figure 8 we plot the hydrodynamization time, $t_{\textrm{hyd}}$, (red dashed line with stars for each run) and the EoSization time, $t_{\rm EoS}$, (purple full line with dots for each run) for the different non-conformal theories. There are two observable effects. First, the hydrodynamization time increases with the non-conformality. Second, hydrodynamization can happen before EoSization. The conformal value of the hydrodynamization time is indicated in each panel of figure 8 with a horizontal line at $t_{\textrm{hyd}}T_{\textrm{hyd}}=0.56$. For a slightly non-conformal theory the increase of $t_{\textrm{hyd}}$ with respect to the conformal value is minimal, as illustrated in the $\phi_{\textrm{\tiny M}}=2$ temperature scan. In this case the maximal increase is just a factor of $1.43$ larger than the conformal value. For $\phi_{\textrm{\tiny M}}=3$ ($\phi_{\textrm{\tiny M}}=5$) the increase of the hydrodynamization time is a factor of $2.05$ ($2.38$), and it takes place for a collision for which the ratio $\zeta/\eta$ at the time of hydrodynamization is $0.21$ ($0.31$). The expected increase of the hydrodynamization time is maximal in the $\phi_{\textrm{\tiny M}}=20$ temperature scan. In this case the maximum occurs for a collision with $T_{\textrm{hyd}}/\Lambda=0.184$ and the increase is a factor of $2.6$ with respect to the conformal result. This maximal hydrodynamization time is reached with a bulk viscosity over entropy density ratio of $\zeta/\eta\approx 0.36$. We have verified that the $\phi_{\textrm{\tiny M}}=20$ results are almost identical to those corresponding to $\phi_{\textrm{\tiny M}}=30$ or $\phi_{\textrm{\tiny M}}=100$. This is consistent with the fact that thermodynamic and transport properties such as the bulk viscosity and the speed of sound squared saturate with big positive values of $\phi_{\textrm{\tiny M}}$.

We see that for sufficiently large $\mu/\Lambda$ the EoSization time becomes negative, meaning that the average and the equilibrium pressures differ by less than 10% even before the shocks collide. The reason is simply that in these cases the energy density in the Gaussian tails in front of the shocks, which start to overlap at negative times, becomes much higher than $\Lambda$. At these energy densities the physics becomes approximately conformal and the equation of state becomes approximately valid as a consequence of this symmetry.

The equilibrium pressure and the average pressure are not within $10\%$ of one another for a wide range of runs with $\phi_{\textrm{\tiny M}}=\{20,5,3,2\}$. For runs for which the EoSization criterion is fulfilled at all post-collision times the extracted EoSization time is either null or negative. Those specific runs show negligible non-conformal effects for the created plasma. The reason for the sharp rise of the EoSization times at low temperatures is due to a cut-off effect of the fixed $10\%$ criterion. Runs with a slightly higher temperature do easily reach $>5\%$ non-conformal effects, but do not yet trigger a later EoSization time extraction. The shockwave literature Chesler:2013lia ; Casalderrey-Solana:2013aba typically uses a hydrodynamization criterion between $15\%$ and $20\%$, whereas here we have settled for $10\%$. We found that changing this number implies no qualitative changes to our conclusions.

Reference Attems:2016tby showed that, in the model with $\phi_{\textrm{\tiny M}}=10$, hydrodynamization precedes EoSization for collisions in a certain range of hydrodynamization temperatures. Figure 8 shows that this also happens in the models with $\phi_{\textrm{\tiny M}}=\{20,5,3\}$. For $\phi_{\textrm{\tiny M}}=20$ this ordering is maintained up to the highest hydrodynamization temperature, $T_{\textrm{hyd}}/\Lambda\approx 0.37$. Since the models with $\phi_{\textrm{\tiny M}}=20$ and $\phi_{\textrm{\tiny M}}=5$ have approximately the same thermodynamic properties Attems:2016ugt at $T/\Lambda\approx 0.4$, at the crossing of the EoSization and the hydrodynamization times the bulk viscosity-to-entropy ratio is also approximately the same, $\zeta/\eta\approx 0.22$. For $\phi_{\textrm{\tiny M}}=3$ one notices that at the high-temperature crossing between the two times the ratio is $\zeta/\eta\approx 0.25$. Models with $\phi_{\textrm{\tiny M}}=2$ and $\phi_{\textrm{\tiny M}}=1$ (not shown explicitly) show no crossing. For $\phi_{\textrm{\tiny M}}=2$ the maximal ratio is $\zeta/\eta\approx 0.18$. For $\phi_{\textrm{\tiny M}}=5$, the low-temperature crossing has $\zeta/\eta\approx 0.32$ and for $\phi_{\textrm{\tiny M}}=3$ one gets $\zeta/\eta\approx 0.23$, but the lower crossing is not yet reached with $\phi_{\textrm{\tiny M}}=20$ at the minimal temperature with $\zeta/\eta\approx 0.31$. These differences are explained by the accumulating effects of the bulk viscosity along the entire evolution of the collision. We therefore confirm our prior conservative estimate Attems:2016tby of $\zeta/s\gtrsim 0.025$ in order to have hydrodynamization before EoSization.

Furthermore, the maximal value for the ratio $t_{\rm EoS}/t_{\textrm{hyd}}\approx\{0.684,1.10,2.07,2.87\}$ with $\phi_{\textrm{\tiny M}}=\{2,3,5,20\}$ is reached at $T_{\textrm{hyd}}/\Lambda\approx\{0.275,0.279,0.322,0.346\}$. The resulting temperature of the maximal values, when comparing the different non-conformal theories, increases. This again shows evidence for the accumulating effect of the bulk viscosity. It is important to stress that the non-conformal equation of state has to be taken into account for more than twice the hydrodynamization time with a bulk viscosity over entropy ratio of $\zeta/s\geq 0.025$.

5.3 Dynamics of the scalar condensate

As we have seen, the dynamics of the longitudinal, transverse and average pressures provides information about different process in the evolving plasma. To clarify further the process of EoSization we focus here on the evolution of the scalar expectation value $\mathcal{V}$, since inspection of equation (23) and equation (26) shows that EoSization is in part controlled by how $\mathcal{V}$ approaches its equilibrium value. In figure 9 we show the spacetime evolution of the condensate for a characteristic collision.

In analogy with the hydrodynamization and EoSization times, we define the condensate relaxation time $t_{\textrm{cond}}$ as the time beyond which the normalized difference of the true expectation value of the scalar operator, $\mathcal{V}$, and its equilibrium value $\mathcal{V}_{\rm eq}$, is less than 10%:

[TABLE]

This time is a measure of how fast this one-point function reaches its equilibrium value. We explore $t_{\textrm{cond}}$ for the different collision configurations studied in the previous section.

The scalar condensate is fully out of equilibrium for most of the studied shockwave collisions. In most cases, the condensate $\mathcal{V}$ takes a much longer time to equilibrate than the system takes to hydrodynamize. With three different relaxation times—hydrodynamization, EoSization and condensate relaxation—one can in theory find six possible orderings between these times. However, we have found no configuration in which hydrodynamization comes last. Since there seems to be no obstacle of principle for this, the reason is presumably that our collisions do not generate a sufficiently large anisotropy; we will come back to this point in section 6. Thus, our models give rise to the following four orderings:

EoSization $\to$ Hydrodynamization $\to$ Condensate relaxation, 2. 2.

Hydrodynamization $\to$ EoSization $\to$ Condensate relaxation, 3. 3.

Hydrodynamization $\to$ Condensate relaxation $\to$ EoSization, 4. 4.

Condensate relaxation $\to$ Hydrodynamization $\to$ EoSization.

Each of these cases is illustrated by one of the plots in figure 10.

The fact that hydrodynamization and EoSization can happen in any order was the main result of Attems:2016tby . Here we see that the situation is richer once condensate relaxation is included. We note in figure 10 that the energy of the collision, $\mu/\Lambda$, or equivalently the hydrodynamization temperature, $T_{\textrm{hyd}}/\Lambda$, decrease monotonically from case 1 to case 4. The reason is that at $T\gg\Lambda$ the condensate grows as $\Lambda\mathcal{V}\sim\Lambda^{2}T^{2}$ Attems:2016ugt , whereas the stress tensor grows as $T^{4}$. As a consequence the relative magnitude of the $\mathcal{V}$-induced correction in the average pressure through the Ward identity (23) decreases and the dynamics of the condensate decouples from the dynamics of the stress tensor. Indeed, using equations (23), (25) and (26) we see that

[TABLE]

For EoSization to take place the left-hand side must be small in units of $\bar{P}$. Dividing this equation by $\bar{P}$ and using the scalings above we find that at high temperature

[TABLE]

In this regime the scalar condensate can still very far from its equilibrium value according to our criterion (89), since all terms in (89) scale as $\Lambda T^{2}$. In conclusion, at high temperature the equation of state approaches the conformal equation of state and the Ward identity is no impediment for the system to EoSize (and hydrodynamize) while the scalar condensate is still far from its equilibrium value. This possibility is realized in the first two plots of figure 10. Note that the hydrodynamization temperature in these cases is certainly not asymptotically high, but it is higher than the temperature at which the non-conformal effects are maximal, which for $\phi_{\textrm{\tiny M}}=20$ is $T\sim 0.2\Lambda$, as indicated in the caption of figure 3. This seems to suffice for the asymptotic argument above to apply.

In contrast, at $T\sim 0.2\Lambda$ the product $\Lambda\mathcal{V}$ can be numerically larger than $3\bar{P}$. For this reason it is possible for the right-hand side of (90) to be smaller than 0.1 in units of $\Lambda\mathcal{V}$ while the left-hand side is larger than 0.1 in units of $\bar{P}$. This is why at temperatures at which non-conformal effects are sufficiently large scalar relaxation can precede EoSization (and also hydrodynamization, since the latter can precede EoSization). This is illustrated by the last two plots in figure 10, for which the hydrodynamization temperature is close to the value at which non-conformal effects are maximal.

In figure 11 we explore the relaxation dynamics for two collisions with the same incident transverse energy, $\mu/\Lambda=0.62$, but with different widths, $\mu\omega=0.12$ (left) and $\mu\omega=0.30$ (right). In both cases, at late times $\mathcal{V}$ approaches its equilibrium value (green dashed) from above. As in figure 10, we see that the equilibrium value $\mathcal{V}(\mathcal{E})$ begins to rise before $t=0$ and reaches its maximum shortly after $t=0$. This is simply because this value tracks the energy density, which begins to rise before $t=0$ because of the forward tails of the Gaussian shocks. Instead, the true condensate would be exactly undisturbed by a single shock, and therefore it begins to respond only once a significant amount of collision dynamics has taken place. For this reason, the true condensate begins to rise almost exactly at $t=0$. Shortly after the collision the spike in the equilibrium value reflects the initial large energy density of the passing shocks, which is larger in the narrower shocks since the transverse energy density is fixed. In contrast, the peak in the true condensate is very similar in both collisions. As in the conformal case Casalderrey-Solana:2013aba , the final hydrodynamization temperature is mostly determined by the transverse energy scale, and therefore $T_{\textrm{hyd}}$ is almost identical for the two collisions. This is remarkable, since $\zeta/\eta$ at that $T_{\textrm{hyd}}$ is almost maximal, indicating large non-conformal effects. We observe that a similar statement holds true for the EoSization time, which is essentially the same in both cases, and less accurately but still approximately so for the relaxation times of the scalar condensate. It may be possible to understand these effects as finite-resolution effects, as discussed in Casalderrey-Solana:2013sxa .

On general grounds, in a CFT one would expect the time at which the true condensate reaches its peak value, $t_{\textrm{peak}}$, to be given by

[TABLE]

with $c$ an order-one constant. The intuitive reason on the gravity side is that it takes a time of order $1/\pi T_{\textrm{hyd}}$ for the effects of the dynamics near the horizon that forms deep in the bulk when the shocks collide to reach the boundary. This delay is also observed in e.g. the true drag force on a quark compared to the force that it would experience in an equilibrium plasma with the same instantaneous energy density Chesler:2013urd ; Lekaveckas:2013lha . In a non-conformal theory one would expect $c$ to be constant for high-energy collisions in which non-conformal effects are small but to deviate from a constant for collisions in which non-conformal effects are significant. These expectations are confirmed in our model, as illustrated by figure 14, where we plot $\pi t_{\textrm{peak}}T_{\textrm{hyd}}$ and, for comparison, also $\pi t_{\textrm{peak}}T_{\textrm{cond}}$ and $\pi t_{\textrm{peak}}T_{\textrm{EoS}}$.

We see that the latter two vary significantly as a function of the collision energy (represented here by its proxy, $T_{\textrm{hyd}}$) and do not become constant at high energies. Also, in these cases one must bear in mind that $T_{\textrm{cond}}$ and $T_{\textrm{EoS}}$ are only well defined when hydrodynamization precedes scalar relaxation and EoSization, respectively. In contrast, we see that $\pi t_{\textrm{peak}}T_{\textrm{hyd}}$ does approach a constant of order $c\simeq 0.8$ at high energies, and that it deviates slightly from it at low energies.

In figure 12 we compare the hydrodynamization time and the condensate relaxation time, in units of $\Lambda^{-1}$, as a function of the hydrodynamization temperature. Both times attain their maximum values at the lowest temperatures we were able to probe, where non-conformal effects are large. Comparing different theories, we see that the maximal $t_{\textrm{hyd}}\Lambda$ happens for $\phi_{\textrm{\tiny M}}=3$. For $\phi_{\textrm{\tiny M}}=20$ we observe a crossing of the scalar relaxation time and the hydrodynamization time, as illustrated above in figure 10, meaning that condensate relaxation can precede hydrodynamization or vice versa. In contrast, models with $\phi_{\textrm{\tiny M}}=\{5,3,2\}$ show scalar relaxation times that are always significantly longer than the corresponding hydrodynamization times. In particular, for small non-conformality (small $\phi_{\textrm{\tiny M}}$) and small temperatures the condensate may still be out of equilibrium at hydrodynamization. Also for small non-conformality ($\phi_{\textrm{\tiny M}}=2$) the oscillations in the scalar condensate cause jumps in the relaxation times extracted with the constant criterion (89). From this scan we can extract two characteristic numbers: the maximal value of $t_{\textrm{cond}}\Lambda\approx 18.8$ is reached at low temperatures with $\phi_{\textrm{\tiny M}}=3$, whereas the maximal value for the ratio $t_{\textrm{cond}}/t_{\textrm{hyd}}\approx 6.09$ occurs for $\phi_{\textrm{\tiny M}}=3$ and is reached at $T_{\textrm{hyd}}/\Lambda\approx 0.81$.

One conclusion of figure 12 is that both the condensate relaxation time and the hydrodynamization time, when measured in units of the intrinsic scale in the theory, decrease as the energy of the collision, or equivalently the hydrodynamization temperature. In fact, these values approach zero at asymptotically high energies, as is clear from the top-right panel in figure 12, where we have extended the range of the horizontal axis to high values in order to illustrate this effect. In figure 13 we show these times measured in units of the hydrodynamization temperature itself. These plots clearly show how at high temperatures the systems behaves effectively as a conformally invariant system. Indeed, if $T\gg\Lambda$ the temperature becomes the only relevant scale and both $t_{\textrm{cond}}T_{\textrm{hyd}}$ and $t_{\textrm{hyd}}T_{\textrm{hyd}}$ approach constant values. Furthermore, these asymptotic values are the same in all four models, which reflects the fact their UV properties are identical. Nevertheless, the temperature at which this asymptotic sets in depends on the model. As discussed around equation (90), at high temperatures the dynamics of the condensate decouples from the dynamics of the stress tensor. The fact that in this asymptotic regime $t_{\textrm{cond}}T_{\textrm{hyd}}$ is 5.18 times larger than $t_{\textrm{hyd}}T_{\textrm{hyd}}$ explicitly shows that a hydrodynamized plasma can be far from equilibrium, since between $t_{\textrm{hyd}}$ and $5t_{\textrm{hyd}}$ hydrodynamics provides a good description of the stress tensor but the expectation value of the scalar operator is still far from its equilibrium value.

5.4 Rapidity profile

Up to now we have focused on the mid-rapidity region, $z=0$. We will now study the energy deposition along the collision axis. To make contact with hydrodynamic simulations of ultra-relativistic heavy ion collisions, we explore the local energy density in the fluid rest frame, $\mathcal{E}_{\rm loc}$, at a fixed proper time $\tau=t_{\textrm{hyd}}$, with $t_{\textrm{hyd}}$ the hydrodynamization time at $z=0$, as a function of the spacetime rapidity $y$, with

[TABLE]

In figure 15 we show the rapidity distribution of the local energy density (normalised to the central energy density) for different collision energies in the $\phi_{\textrm{\tiny M}}=20$ model. For comparison, we also show the same distribution for collisions in $\mathcal{N}=4$ SYM Casalderrey-Solana:2013aba ; Chesler:2015fpa . As in that conformal case, here the deposited energy density exhibits a strong rapidity dependence which is well approximated by a Gaussian within a 1-2% accuracy. The width of the Gaussian, however, depends on the transverse energy scale $\mu/\Lambda$. For smaller values of the collision energy, the hydrodynamization temperature is also smaller and the non-conformal behaviour is more pronounced. As the collision energy increases, the rapidity width of the energy deposition grows, approaching the conformal distribution asymptotically at large collision energies. Although this energy density profile is controlled by non-hydrodynamized dynamics, the observed dependence with $T_{\textrm{hyd}}$ is consistent with the expectations from bulk viscosity. Similarly to the reduction of transverse expansion observed in hydrodynamic simulations of ultra-relativistic plasmas Ryu:2015vwa , the bulk viscosity reduces the longitudinal pressure, reducing the transport of energies at large rapidities. It is interesting that the increase in the width of the energy rapidity profile is in qualitative agreement with the rapidity distribution of matter in heavy ion collisions as a function of $\sqrt{s}$.

Despite the fact that the system is manifestly non-boost invariant, specially in the most non-conformal region, to a very good approximation the initial velocity field at hydrodynamization is. In figure 16 we show the component of the velocity field along the proper time direction,

[TABLE]

as a function of rapidity for several proper times after $\tau=t_{\textrm{hyd}}$. The fact that $u^{\tau}$ is so close to 1 in the two units of rapidity that we have plotted shows that the four-velocity field is well aligned with the proper time direction, with small deviation at the sub-percent level. This result was first observed in shockwave collisions in conformal gauge theories Chesler:2015fpa for a variety of initial Gaussian shock widths. What we are observing here is that this result survives the introduction of large non-conformal effects, even though those same effects do cause a narrowing of the energy density rapidity distribution. At later times, the fact that the rapidity deposition of energy is not boost-invariant will change the velocity field, increasing the rapidity component of the velocity; nevertheless, this change is completely predicted by hydrodynamics. Our simulations imply than even for non-conformal dynamics, in order to completely predict the stress tensor dynamics in different configurations only the rapidity distribution of energy density needs to be specified at hydrodynamization, since the initial velocity field is given, to a very good approximation, by $u^{\tau}=1$. This observation can be translated into consequences for hydrodynamic modellers of heavy ion collisions: even for configurations with significant rapidity dependence the initialization of the velocity field after the collision in a boost-invariant manner is well supported by our simulations.

6 Discussion

Following the procedure described in section 3 we have simulated 565 shockwave collisions in the gravity-plus-scalar models of Attems:2016ugt . Via holography, we have used the results to perform a thorough analysis of the out-of-equilibrium dynamics of the dual set of non-conformal gauge theories with different degrees of non-conformality.

One of the most astonishing results of this analysis is the tremendous success of hydrodynamics to describe the out-of-equilibrium evolution. This fact has been extensively studied in many settings for conformal theories in the past, where it was found that, at strong coupling, hydrodynamization typically precedes isotropization. We have verified that this is also the case in our non-conformal plasmas. In fact, in all the collisions that we have examined we have found that isotropization is always the last process to take place of the four that we have considered. To illustrate this quantitatively, we note that the ratio $P_{T}/P_{L}$ at the latest of the three equilibration times shown in each of the four panels of figure 10 is $2.0,1.9,1.9$ and $1.7$, respectively. In other words, at the latest equilibration time shown in the panels the transverse pressure is still at least 70% larger than the longitudinal one, indicating that the plasma is still significantly anisotropic.

It is remarkable that hydrodynamics works so well even with a non-trivial equation of state. In particular, in our most non-conformal models the number of degrees of freedom changes by several orders of magnitude between the high- and the low-temperature phases—three in the $\phi_{\textrm{\tiny M}}=10$ case shown in figure 1(right). Yet, the dynamics of the system is well described very soon after the collision by a hydrodynamic expansion around this non-trivial equation of state. The break-down of the different components of the hydrodynamic estimator displayed in figure 7 clearly illustrates this point. The success is such that in our extensive exploration of the parameter space of non-conformal collisions we have never encountered a case in which the hydrodynamization time exceeds the value in the conformal case by factor larger than 2.6. This is in agreement with the expectations based on the near-equilibrium analysis in terms of quasi-normal modes Buchel:2015saa ; Attems:2016ugt .

The success of hydrodynamics is even more surprising in cases in which hydrodynamization precedes all other equilibration processes. In these cases, which correspond to the panels 2 and 3 of figure 10, hydrodynamics provides an accurate description of the evolution of the plasma despite the fact that “everything else is far from equilibrium”, meaning that the average pressure and the condensate are still far from their equilibrium values and the plasma is still highly anisotropic.

Focusing on the particular ordering of hydrodynamization and EoSization, our results confirm that the former precedes the latter as long as the system is sufficiently non-conformal. What is perhaps surprising is that, as measured by the bulk viscosity-to-entropy ratio, a “sufficient” degree of non-conformality requires only a fairly moderate value $\zeta/s\gtrsim 0.025$, as estimated in Attems:2016tby . This indicates that similar phenomena may also occur in real-world heavy ion collisions, where both calculations Paech:2006st ; Arnold:2006fz ; Karsch:2007jc ; Denicol:2014vaa and data-driven parametrization Ryu:2015vwa ; Bernhard:2016tnd yield larger values than this estimate in a significant part of the time evolution of the resulting plasma. It would be interesting to extend existing phenomenological studies Torrieri:2008ip ; Monnai:2009ad ; Song:2009rh ; Rajagopal:2009yw ; Dusling:2011fd ; McDonald:2016vlt of the effect of bulk viscosity in heavy ion collisions to investigate the possibility that hydrodynamization may precede EoSization.

Although the Ward identity (23) implies that EoSization and condensate relaxation are related, we have seen that nevertheless these two processes can occur in any ordering. The reason for this is easy to understand in two limits, one in which the temperature is much higher than the intrinsic scale in the theory and another in which it is comparable to this scale. In the first case the different scalings with the temperature of $\bar{P}\sim T^{4}$ and of $\Lambda\mathcal{V}\sim\Lambda^{2}T^{2}$ imply that at high temperature the contribution of the condensate to the Ward identity is subleading, and the dynamics of the stress tensor decouples from that of the condensate. Thus, in this limit the system can EoSize and hydrodynamize while the condensate remains far from equilibrium. This is clearly illustrated by figures 8 and 13, in which we see that in the high-temperature limit

[TABLE]

Our simple example suggests that other one- or higher-point functions of non-conserved operators may take a long time to relax even in an almost-conformal, hydrodynamized and EoSized plasma—for example, a similar delay in the relaxation of fluctuations in a non-conformal plasma undergoing a process of isotropization has been analysed in Chesler:2012zk . This may have important implications for processes depending on non-hydrodynamic properties of the plasma created in heavy-ion collisions, such as emission rates and the reaction of the plasma to probes, which are typically assumed to be quantified in terms of equilibrium plasma properties. It would be interesting to explore the deviations from equilibrium of these phenomenologically relevant quantities with holography.

In the second case, when the temperature is close to the value at which the non-conformal effects are maximal, the value of the pressures and of the condensate are all parametrically the same. However, numerically we find that in some situations $\Lambda\mathcal{V}>3\bar{P}$ at $t=t_{\textrm{cond}}$. This means that at this time the condensate is within 10% of its equilibrium value but its contribution though the Ward identity still causes a larger-than-10% deviation between the average pressure and its equilibrium value.

In section 5.3 we determined the possible orderings once the three times $t_{\textrm{hyd}},t_{\textrm{EoS}}$ and $t_{\textrm{cond}}$ are simultaneously considered. We found that in our model only the four orderings illustrated in figure 10 seem to be realized. Out of the six orderings that are logically possible, the two missing ones are

EoSization $\to$ Condensate relaxation $\to$ Hydrodynamization, 2. 6.

Condensate relaxation $\to$ EoSization $\to$ Hydrodynamization,

namely the two orderings in which hydrodynamization happens last. Presumably the reason is simply that our collisions do not produce a plasma that is sufficiently anisotropic. Indeed, EoSization, and indirectly condensate relaxation through the Ward identity, is controlled by the bulk gradient corrections to the equilibrium pressure. Therefore it is conceivable that, in a dynamical situation in which shear corrections are much larger than bulk corrections, the average pressure and the condensate may relax to their equilibrium value at a time at which the difference between the pressures is still not well predicted by hydrodynamics.

Throughout the paper we have adopted a “10%” criterion to define the hydrodynamization, EoSization and condensate relaxation times in (85), (88) and (89). Since this criterion is arbitrary, it is interesting to ask what happens if the 0.1 in these equations is replaced by, say, 0.15 or 0.2. The result is summarised in figure 17, which shows the three equilibration times with a 15% criterion or a 20% criterion.

Comparison with the 10% criterion yields the following qualitative conclusions:

• •

Although not shown in figure 17, the isotropization time is still the longest.

• •

The conclusion that the two times $t_{\textrm{EoS}}$ and $t_{\textrm{cond}}$ can occur in any ordering remains true for any criterion.

• •

The three times $t_{\textrm{hyd}}$, $t_{\textrm{EoS}}$ and $t_{\textrm{cond}}$ can still occur in several different orderings, but which specific orderings are realized depends on the criterion. With the 15% criterion these orderings are 1, 4, 5 and 6, whereas with the 20% criterion we get 4, 5 and 6, and almost 1. However, it is possible that in a model with more general dynamics (larger gradients, larger bulk viscosity, etc) all possible orderings may be realized for a given criterion.

• •

Hydrodynamization can precede EoSization with a 15% criterion (as with the 10% criterion) but not with a 20% criterion. This is not surprising since the moderate bulk viscosity of our model is only able to produce moderate deviations of the average pressure from its equilibrium value. For example, for the collisions examined in Attems:2016tby this deviation at $t_{\textrm{hyd}}$ was about 18%.

• •

With the two new criteria there is no collision in which hydrodynamization precedes all other equilibration processes. In other words, the orderings 2 and 3 are only realized with a 10% criterion. Again, we expect that these orderings would be realized for less stringent criteria in a model with more general dynamics.

For simplicity, we have considered a model with a single scalar field, i.e. we have focused on the dynamics in the sector in which only the conserved stress tensor and one non-conserved scalar operator are included. In a model in which several non-conserved operators are considered, the Ward identity (23) would relate the trace of the stress tensor to the sum of the sources times the condensates of all the non-conserved operators. Therefore any other linearly independent combination of these operators would be unconstrained by the Ward identity. It would be interesting to study a model of this type, since presumably the dynamics would be even richer than in our one-field model.

We have also studied the post-collision deposition of energy as a function of rapidity. As in conformal collisions, the initial longitudinal flow field is, to a surprising degree of accuracy, boost invariant. Even for the most non-conformal collisions that we have studied the size of the longitudinal gradients is insufficient to alter the longitudinal expansion of the created matter. This, together with similar results found in conformal collisions Chesler:2015fpa , may be viewed as dynamical evidence in support of initializing hydrodynamic simulation of heavy ion collisions with a boost invariant flow field, even at relatively small collision energies.

Concerning the rapidity profile of the energy density, we have found that non-conformal effects make the rapidity distribution of the collision debris narrower than for a conformal collision with identical collision parameters. On the one hand, this is perhaps unsurprising since, at least in the hydrodynamized regime, this may be expected from the friction induced by the bulk viscosity. On the other hand, this feature highlights a main difference between the non-conformality of our model and that of QCD: In our model the theory flows at high energies to a strongly coupled fixed point, whereas QCD flows to a free fixed point. In other words, as any model that can be fully described by classical gravity Mateos:2011bs , our model fails to reproduce asymptotic freedom. In QCD this property makes the energy rapidity profile broader and broader as the collision energy increases, since in this asymptotic regime the physics mostly responsible for setting this profile is pre-hydrodynamic weakly coupled physics. In contrast, in our model the rapidity profile saturates at high energies to that of a strongly coupled conformal theory, which is known to result in a narrower profile Chesler:2015fpa ; vanderSchee:2015rta . It would be interesting to develop hybrid approaches, perhaps along the lines of Casalderrey-Solana:2014bpa ; Iancu:2014ava , able to address separately the strongly coupled regime at energies around $\Lambda$ via holography, and the weakly coupled regime at much higher energies via a different description.

Appendix A Matching the hyperbolic equations

Let us consider the evolution equations for the metric variable $B$ (the corresponding ones for $\phi$ are entirely analogous). As outlined in section 3.4, we have two grids, grid1 and grid2, where we need to evolve $B_{g_{1}}$ and $B_{g_{2}}$ (algebraically related with the metric coefficient $B$). The two grids can overlap, but we assume for simplicity that they merely touch at point $u=u_{0}$, i.e., grid1 covers the region $u\in[0,u_{0}]$ and grid2 covers $u\in[u_{0},u_{\rm max}]$, the AdS boundary being at $u=0$.

From equation (42), the evolution equation for $B_{g_{1}}$ (the case for grid2 is entirely analogous) takes the form

[TABLE]

which has the generic form

[TABLE]

with

[TABLE]

$c_{g_{1}}(u,z)$ is locally the propagation speed and in the vicinity of $u=u_{0}$ we can formally write the solution of this equation (ignoring from now on the $z$ dependence) as

[TABLE]

for any given function $f$.

Therefore, for $c_{g_{1}}>0$ ($c_{g_{1}}<0$), information is propagating from grid2 to grid1 (grid1 to grid2). In order to consistently solve this system, the procedure will then be to use equation (A) (and the corresponding one for $B_{g_{2}}$ on grid2) on all interior points (i.e., points where $u\neq u_{0}$) and for the junction point $u=u_{0}$ one checks the propagation speed at each $z$ point and copies the values according to the propagation direction:

• •

$c_{g_{1}}>0$

[TABLE]

i.e., we copy the modes leaving grid2 to grid1.

• •

$c_{g_{1}}<0$

[TABLE]

i.e., we copy the modes leaving grid1 to grid2.

Acknowledgements.

We thank A. Buchel, P. Chesler, P. Figueras, G. Horowitz, J. Rocha, P. Romatschke, W. van der Schee, B. Schenke, U. Sperhake, K. Tywoniuk, R. Venugopalan and U. Wiedemann for discussions. We thank the supercomputer MareNostrum at the Barcelona Supercomputing Center for providing computational resources (project no. ub65). The work of MA has been supported by a Marie Sklodowska-Curie Individual Fellowship of the European Commission’s Horizon 2020 Programme under contract number 658574 FastTh. JCS is a Royal Society University Research Fellow. JCS was also supported by a Ramón y Cajal fellowship, by the Marie Curie Career Integration Grant FP7-PEOPLE-2012-GIG-333786 and by the Spanish MINECO through grant FPA2013-40360-ERC. CFS and DS acknowledge the support from contracts ESP2013-47637-P and ESP2015-67234-P (Spanish Ministry of Economy and Competitivity of Spain, MINECO). MZ acknowledges support through the FCT (Portugal) IF programme, IF/00729/2015. We also acknowledge funding from grants MEC FPA2013-46570-C2-1-P, MEC FPA2013-46570-C2-2-P, MDM-2014-0369 of ICCUB, 2014-SGR-104, 2014-SGR-1474, CPAN CSD2007-00042 Consolider-Ingenio 2010, and ERC Starting Grant HoloLHC-306605.