Coupling Constant Corrections in a Holographic Model of Heavy Ion Collisions with Nonzero Baryon Number Density

TL;DR

This paper uses holography to model heavy ion collisions with nonzero baryon density, analyzing how finite gauge coupling corrections affect the distribution of baryon charge after collisions.

Contribution

It introduces four-derivative gravitational terms in holographic models to study finite coupling effects on baryon density dynamics in heavy ion collisions.

Findings

Less baryon density is stopped at mid-rapidity with reduced coupling.

More baryon density moves near the lightcone at lower coupling.

Finite coupling corrections influence charge distribution post-collision.

Abstract

Sufficiently energetic collisions of heavy ions result in the formation of a droplet of a strongly coupled liquid state of QCD matter known as quark-gluon plasma. By using gauge-gravity duality (holography), a model of a rapidly hydrodynamizing and thermalizing process like this can be constructed by colliding sheets of energy density moving at the speed of light and tracking the subsequent evolution. In this work, we consider the dual gravitational description of such collisions in the most general bulk theory with a four-derivative gravitational action containing a dynamical metric and a gauge field in five dimensions. Introducing the bulk gauge field enables the analysis of collisions of sheets which carry nonzero "baryon" number density in addition to energy density. Introducing the four-derivative terms enables consideration of such collisions in a gauge theory with finite gauge coupling, working perturbatively in the inverse coupling. While the dynamics of energy and momentum in the presence of perturbative inverse-coupling corrections has been analyzed previously, here we are able to determine the effect of such finite coupling corrections on the dynamics of the density of a conserved global charge, which we take as a model for the dynamics of nonzero baryon number density. In accordance with expectations, as the coupling is reduced we observe that after the collisions less baryon density ends up stopped at mid-rapidity and more of it ends up moving near the lightcone.

Equations

$$S=\frac{1}{2{\kappa }_5^2}\int d^{5}x\sqrt{-g}\left(R-2\Lambda -\frac{L^2}{4}{F}_{\mu \nu }{F}^{\mu \nu }\right),$$

$$S=\frac{1}{2{\kappa }_5^2}\int d^{5}x\sqrt{-g}\left(R-2\Lambda +\frac{1}{8}{\alpha }^{\prime 3}\zeta (3)\mathcal{W}\right),$$

$$\displaystyle S=\frac{1}{2\kappa_{5}^{2}}\int d^{5}x\sqrt{-g}\Bigl{(}R-2\Lambda+\alpha_{1}L^{2}R^{2}+\alpha_{2}L^{2}R_{\mu\nu}R^{\mu\nu}+\alpha_{3}L^{2}R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}\Bigr{)}.$$

$$S=\frac{1}{2{\kappa }_5^2}\int d^{5}x\sqrt{-g}\left[R-2\Lambda +\frac{{\lambda }_{GB}{L}^2}{2}\left(R^2-4R_{\mu \nu }{R}^{\mu \nu }+R_{\mu \nu \rho \sigma }{R}^{\mu \nu \rho \sigma }\right)\right],$$

$$S$$

$$\displaystyle+\beta\epsilon^{2}L^{4}\left(RF_{\mu\nu}F^{\mu\nu}-4R^{\mu\nu}F_{\mu\rho}F_{\nu}^{~{}\rho}+R^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)\biggr{]},$$

$${ϵ}^2=c_{ϵ}{\lambda }_{GB},\beta =c_{\beta }{\lambda }_{GB},$$

$$T^{\mu \nu }$$

$$J^{\mu }$$

$$T^{\mu \nu }$$

$$J^{\mu }$$

$$\mathcal{E}+p=sT+\mu \rho ,$$

$$\frac{\eta }{s}=\frac{1}{4\pi }\left(1-4{\lambda }_{GB}\right),$$

$$\frac{\eta }{s}=\frac{1}{4\pi }\left[1-4{\lambda }_{GB}+\mathcal{O}\left({\lambda }_{GB}^2\right)+\mathcal{O}\left({\lambda }_{GB}{\mu }^2/T^2\right)+\mathcal{O}\left(\beta {\mu }^2/T^2\right)+\dots \right];$$

$$\rho =\frac{1}{2}{\pi }^2{T}^2\mu \left(1+3{\lambda }_{GB}+16\beta \right),$$

$$\chi\equiv\frac{\partial\rho}{\partial\mu}\Big{|}_{\mu=0}=\frac{1}{2}\pi^{2}T^{2}\left(1+3\lambda_{\scriptscriptstyle GB}+16\beta\right).$$

$$\frac{\chi }{{\chi }_0}=\frac{1}{2}\left(1+3{\lambda }_{GB}+16\beta \right).$$

$$\beta \in \left[0.04,0.15\right].$$

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R-\frac{6}{L^2}{g}_{\mu \nu }=T_{\mu \nu }^{\mathrm{A}}+T_{\mu \nu }^{\mathrm{GB}}+T_{\mu \nu }^{\beta },$$

$$T_{\mu \nu }^{\mathrm{A}}$$

$$T_{\mu \nu }^{\mathrm{GB}}$$

$$-{\lambda }_{GB}\left(R{R}_{\mu \nu }-2R_{\mu \alpha }{R}_{\nu }^{\alpha }-2R_{\mu \alpha \nu \beta }{R}^{\alpha \beta }+R_{\mu \alpha \beta \gamma }{R}_{\nu }^{\alpha \beta \gamma }\right),$$

$$T_{\mu \nu }^{\beta }$$

$$-2{ϵ}^2{L}^4\beta [g_{\mu \nu }{R}^{\alpha \beta }{F}_{\alpha \lambda }{F}_{\beta }^{\lambda }-4R_{\mu \alpha }{F}_{\nu \beta }{F}^{\alpha \beta }-2R_{\alpha \beta }{F}_{\mu }^{\alpha }{F}_{\nu }^{\beta }-\square \left(F_{\mu \alpha }{F}_{\nu }^{\alpha }\right)$$

$$-g_{\mu \nu }{\nabla }_{\alpha }{\nabla }_{\beta }\left(F_{\lambda }^{\alpha }{F}^{\beta \lambda }\right)+{\nabla }_{\alpha }{\nabla }_{\mu }\left(F_{\nu \beta }{F}^{\alpha \beta }\right)+{\nabla }_{\alpha }{\nabla }_{\nu }\left(F_{\mu \beta }{F}^{\alpha \beta }\right)]$$

$$+\frac{{ϵ}^2{L}^4\beta }{2}\left[g_{\mu \nu }{R}^{\alpha \beta \gamma \delta }{F}_{\alpha \beta }{F}_{\gamma \delta }-6R_{\mu \alpha \beta \gamma }{F}_{\nu }^{\alpha }{F}^{\beta \gamma }-4{\nabla }^{\beta }{\nabla }^{\alpha }\left(F_{\mu \alpha }{F}_{\nu \beta }\right)\right].$$

$${\nabla }_{\nu }{F}^{\mu \nu }=4\beta L^2{\nabla }_{\nu }\left(R{F}^{\mu \nu }+R^{\mu \nu \rho \sigma }{F}_{\rho \sigma }-2R^{\mu \rho }{F}_{\rho }^{\nu }+2R^{\nu \rho }{F}_{\rho }^{\mu }\right).$$

$$L=1+\frac{1}{2}{\lambda }_{GB}+\mathcal{O}({\lambda }_{GB}^2).$$

$$d{s}^2$$

$$A$$

$$h_0(x_{\pm })=\frac{m^3}{\sqrt{2\pi w^2}}\exp \left(-\frac{x_{\pm }^2}{2w^2}\right),a_{\pm (0)}(x_{\pm })=\frac{h(x_{\pm })}{m},$$

$$d{s}^2=2dtdr-C(t,r,z)d{t}^2+2F(t,r,z)dzdt+S(t,r,z{)}^2\left(e^{-2B(t,r,z)}d{z}^2+e^{B(t,r,z)}d{\mathbf{x}}_{\perp }^2\right),$$

$$A=A_t(r,t,z)dt+A_z(r,t,z)dz.$$

$$C(r,t,z)$$

$$A_{\mu }(r,t,z)$$

Full text

††preprint: MIT-CTP/5136

1. I Introduction

2. II The holographic model

3. III Analysis of the model and details of the numerical setup

4. IV Results

5. IV.1 Colliding thick and thin sheets

6. IV.2 Effect of $\beta$ and $\lambda_{\scriptscriptstyle GB}$ on baryon number dynamics

7. IV.3 Spacetime rapidity distribution

8. IV.4 Hydrodynamization

9. IV.5 Subleading $\mu/T$ effects for thick sheets

10. IV.6 Asymmetric collisions of thin sheets

11. V Discussion and Outlook

12. A General four-derivative action and field redefinitions

I Introduction

Previous authors Casalderrey-Solana et al. (2016) have completed a holographic analysis of collisions of sheets of energy density, incident at the speed of light, that carry a nonzero density of a conserved global charge in ${\cal N}=4$ $SU(N_{c})$ supersymmetric Yang-Mills (SYM) theory in the limit of a large number of colors $N_{c}$ and infinite ’t Hooft coupling $\lambda\equiv g^{2}N_{c}$, with $g$ the gauge coupling. Although ${\mathcal{N}}=4$ SYM theory differs from QCD, there are many similarities between the strongly coupled liquid phase that it features at any nonzero $T$ and the strongly coupled liquid quark-gluon plasma phase of that is found in QCD over a range of temperatures. The range extends well above the crossover temperature at which ordinary hadrons form and includes the range of temperatures explored via heavy ion collisions at RHIC and the LHC. Because ${\mathcal{N}}=4$ SYM theory has a dual gravitational description that allows reliable insights into its properties and dynamics at strong coupling it has often been used as a toy model with which to mimic the QCD dynamics via which QGP is formed and probed. (For reviews, see Refs. Casalderrey-Solana et al. (2011); DeWolfe et al. (2014); Chesler and van der Schee (2016); Busza et al. (2018).) The nonzero density in the ${\mathcal{N}}=4$ SYM calculation, which is that associated with a global $U(1)_{R}$ symmetry of the gauge theory, can be used as a toy model for baryon number density in QCD. Hence, the collisions studied in Ref. Casalderrey-Solana et al. (2016) can be thought of as a toy model for heavy ion collisions with nonzero baryon number density. The authors of Casalderrey-Solana et al. (2016) found that the fluid energy, momentum, and baryon current all become well described by charged hydrodynamics at roughly the same time after the collision. They find very significant stopping of the baryon number, with the baryon density after the collision piled up at mid-rapidity. These results are quite different from what is seen in ultrarelativistic heavy ion collisions. In the holographic calculations of Ref. Casalderrey-Solana et al. (2016), the baryon number distribution in the hydrodynamic fluid produced after the collision has a rapidity distribution that is narrowly sharply peaked around zero, very different than what is seen in high energy heavy ion collisions in which the baryon number ends up near the lightcone, losing only about two units of rapidity rather than getting stopped at midrapidity. (For a review, see Ref. Busza et al. (2018).) Because the physics in the holographic model is strongly coupled at all length scales whereas in QCD it is weakly coupled at short length scales, meaning at the earliest moments in a collision, the baryon stopping is much greater in the holographic model. This study reveals the most striking qualitative difference between holographic collisions and high energy heavy ion collisions. Our goal in the present study is to test the hypothesis that this difference can be attributed to the weakness of the QCD coupling during the earliest moments of a collision.

Before continuing, we pause to review some basic facts known from experimental measurements, to provide some context. (This paragraph follows the discussion in Ref. Busza et al. (2018) closely.) On average, in high energy nucleus-nucleus collisions, each participating nucleon loses about two units of rapidity Busza and Ledoux (1988), which is to say 85% of its incident energy goes into the creation of the droplet of QGP that eventually ends up becoming a very large number of particles, up to 30,000 in PbPb collisions at the LHC, and their kinetic energy. These particles, and this energy, end up in a broad distribution centered at mid-rapidity. However, the net proton (number of protons minus antiprotons) rapidity distribution in ultrarelativistic $AA$ collisions (which traces where the incident baryon number ends up) is not centered at mid-rapidity. It has a double-hump structure Arsene et al. (2009), with each hump consisting of hot baryon-number-rich matter moving at a speed of about two units of rapidity below that of the incident beam, and having a net baryon density of about 5-10 times that of normal nuclear matter Busza and Goldhaber (1984). At LHC, the beam rapidity is at $y=8.5$, denoting rapidity by $y$. So, in the center-of-mass frame, this hot baryon rich matter is found at rapidities $y\sim 6.5$ and $y\sim-6.5$. It is interesting to consider how the collision looks upon boosting to the frame in which $y=6.5$ is at rest. In this frame, an incident disc that is Lorentz contracted by a factor of about $\cosh(2)$ is hit by a disc that is Lorentz contracted by about a factor of $\cosh(15)$ and brought approximately to rest, compressed by roughly $2\cosh(2)\approx 7.5$. A further consequence of these considerations is that the maximum value of the net baryon density at mid-rapidity is produced when heavy ions collide with a lower center-of-mass energy per nucleon pair given by $\sqrt{s_{NN}}\approx 7$ GeV. Above this collision energy, the mid-rapidity net baryon density, and so also the baryon chemical potential in the QGP produced at mid-rapidity, decreases with energy. By top RHIC collision energies ($\sqrt{s_{NN}}=200$ GeV), and even more so for LHC energies ($\sqrt{s_{NN}}=5.5$ TeV), both are essentially zero Abelev et al. (2013). This is all very different from ultrarelativistic holographic collisions in an infinitely strongly coupled gauge theory, in which after the collision, the baryon density ends up in a distribution that is centered on mid-rapidity and that is not dissimilar to the way in which the energy density is distributed Casalderrey-Solana et al. (2016).

To investigate the hypothesis that the distinction between how the baryon number density ends up distributed in rapidity in heavy ion collisions vs. in holographic collisions reflects the weakness of the QCD coupling early in the collision, we shall analyze the dynamics of the baryon number density in holographic collisions upon reducing the gauge coupling to a finite value, which is to say upon increasing the inverse coupling from zero to nonzero. A starting point for such a study can be found in Ref. Grozdanov and van der Schee (2017). These authors have performed a holographic study of coupling-constant-dependence in holographic collisions in a hypothetical quantum field theory at large $N_{c}$ and ’t Hooft coupling $\lambda$, including inverse-$\lambda$ corrections to the large $\lambda$ limit, albeit in collisions with no baryon number. The gauge theory here is typically assumed to be some deformation of ${\cal N}=4$ SYM theory. Such collisions can then be thought of as a toy model for heavy ion collisions with zero baryon number density that includes finite coupling corrections. For the first time, these authors analyzed the implications in collisions of the effects of corrections to the assumption of infinitely strong coupling, including inverse-coupling-constant corrections to leading order. In the dual description, this amounts to colliding sheets of energy density in a gravitational theory with curvature-squared terms. They find that at intermediate coupling, less energy density is stopped at mid-rapidity and more ends up near the lightcone. They also find that when the coupling is decreased to the point that the shear viscosity is increased by 80%, the hydrodynamization time increases, but only by 25%.

Our goal here is to extend the calculation of Ref. Casalderrey-Solana et al. (2016) to include finite coupling corrections. Equivalently, our goal is to extend the calculation of Ref. Grozdanov and van der Schee (2017) to incorporate collisions of sheets of energy and baryon number density. That is, we shall include corrections that are leading order in the inverse-coupling in the analysis of holographic heavy ion collisions with baryon number density.

We find that decreasing the coupling constant results in substantially less baryon stopping, with the baryon density stopped at mid-rapidity after the collision dropping, and with more baryon number density ending up moving near the light cones after the collision. In fact, it appears possible to dial the coupling constant down to a value such that the hydrodynamic plasma that forms at mid-rapidity has almost no baryon density, as is the case in ultrarelativistic heavy ion collisions. Although driving the mid-rapidity baryon density down this far would require an inverse-coupling that is much larger than can be treated in our perturbative approach, our result nevertheless provides substantive support to the hypothesis that the excessive baryon stopping found in holographic collisions is an artifact of the fact that in these models the coupling at very early times cannot be made small whereas the QCD coupling is weak in the earliest moments of an ultrarelativistic heavy ion collision.

A reader who is principally interested in our results, and in seeing how these results support the hypothesis noted above, should skip ahead to Sections IV and V. In Section II we describe the holographic model that we employ, and in Section III we detail how we analyze this model and how we conduct our numerical simulations.

II The holographic model

Quark-gluon plasma (QGP) is a hot, strongly coupled, liquid, deconfined, state of matter with temperature $T$. Generically, in addition it has some nonzero baryon number density $\rho$. The baryon number chemical potential $\mu$ is thermodynamically conjugate to $\rho$. In this work, we begin from a model for the formation and evolution of such a state in a non-Abelian gauge theory (${\cal N}=4$ SYM theory) with a large number of degrees of freedom $N_{c}$ by analyzing its holographically dual gravitational theory. In the simplest holographic model, the combined dynamics of the energy-momentum and baryon number density of a four-dimensional state with an infinite ’t Hooft coupling $\lambda$ is described by an asymptotically Anti-de Sitter (AdS) geometry in five-dimensional Einstein-Maxwell theory. Its bulk action is

[TABLE]

where $\Lambda\equiv-6/L^{2}$ is the negative cosmological constant with $L$ the AdS radius, and where the five-dimensional Newton’s constant $\kappa_{5}$ is inversely proportional to the number of colors in the four-dimensional gauge theory: $\kappa_{5}\propto 1/N_{c}\ll 1$ — a limit which allows for a study of the gravitational theory in the classical approximation to bulk (quantum) gravity. The four-dimensional boundary field theory energy-momentum tensor $T^{\mu\nu}$ is sourced by the five-dimensional metric perturbation and the four-dimensional conserved global current $J^{\mu}$ in the boundary theory is sourced by the gauge field in the bulk, according to the standard holographic dictionary. (See, e.g., Ref. Casalderrey-Solana et al. (2014a).) We will think of the dynamics of $J^{\mu}$ as a model for the dynamics of baryon number current.

The holographic “experiment” of a collision of two sheets of energy — and baryon number — density incident at the speed of light proceeds as follows. First, an initial state with a pair of well-separated (along the direction of their motion) sheets of energy density and baryon number density is prepared on the asymptotic AdS boundary. We shall only consider sheets which are infinite in transverse extent and translationally invariant in directions perpendicular to their direction of motion and whose profiles (in the direction of their motion) are Gaussian. (Both these simplifications can be relaxed Chesler and Yaffe (2015); Chesler (2015); Chesler et al. (2015).) The initial state is then evolved in a way such that the two sheets of energy and baryon number travel towards each other along the beam axis at the speed of light. The sheets collide, resulting in the formation of matter which rapidly hydrodynamizes and eventually thermalizes and equilibrates by forming a charged black hole in the bulk which is the holographic representation of strongly coupled QGP. In the absence of any baryon number density, this setup has been extensively studied at infinite coupling; see for example Refs. Chesler and Yaffe (2011); Grumiller and Romatschke (2008); Casalderrey-Solana et al. (2013); van der Schee et al. (2013); Chesler and Yaffe (2014); Casalderrey-Solana et al. (2014b); Chesler and Yaffe (2015); Chesler et al. (2015); van der Schee and Schenke (2015); Chesler (2015); Waeber et al. (2019) and the reviews in Refs. Casalderrey-Solana et al. (2014a); Chesler and van der Schee (2015); Heller (2016); Busza et al. (2018). Collisions of sheets of energy density with nonzero $\mu$ and consequently with nonzero baryon density $\rho$ were studied via their dual description in terms of the five-dimensional Einstein-Maxwell theory (1) in Ref. Casalderrey-Solana et al. (2016).

While the interactions between quarks and gluons in a realistic QCD plasma are strong, they are not infinite. This provides one motivation to move beyond infinite ’t Hooft coupling $\lambda$. The more important motivation for our consideration is that, regardless of how strong the QCD coupling is in the liquid plasma that forms after a heavy ion collision, at the earliest moments of an ultrarelativistic collision in QCD the coupling is weak. And, the question of where the baryon number ends up at late times depends crucially on what happens during the very earliest moments, long before hydrodynamization. For these reasons, we wish to investigate whether weakening the ’t Hooft coupling $\lambda$ in the holographic model calculation changes where the baryon number ends up at late times in a way that makes it look a little more realistic.

Inverse coupling constant corrections can be incorporated into the holographic calculation by perturbing the bulk theory with a series of higher-derivative terms. A well-understood such (top-down) example is the leading-order ’t Hooft coupling correction to the dynamics of energy and momentum in $SU(N_{c})$, $\mathcal{N}=4$ SYM theory. At zero baryon number density, the leading part of the corrected dimensionally reduced gravity action is

[TABLE]

where $\mathcal{W}$ is a contraction of Weyl tensors to the fourth power and $\zeta(z)$ is the Riemann zeta function. For details, see Refs. Gubser et al. (1998); Pawelczyk and Theisen (1998); Grozdanov et al. (2016). In type IIB string theory, the tension of the fundamental string sets the dual ’t Hooft coupling of $\mathcal{N}=4$ SYM theory, i.e. $\alpha^{\prime}\sim 1/\sqrt{\lambda}$. Thus, a perturbative supergravity expansion in $\alpha^{\prime}$, which starts at $\alpha^{\prime 3}$ gives rise to field theory observables computed in a perturbative expansion to ${\cal O}(1/\lambda^{3/2})$.

In a generic string theory, the corrections to the Einstein-Hilbert action start at ${\cal O}(\alpha^{\prime})$ with curvature-squared terms $\sim R^{2}$ instead of $\mathcal{W}\sim R^{4}$. The most general such action is

[TABLE]

The higher-derivative terms all need to be treated perturbatively. If we were doing a top-down construction starting from a known string theory with a known field theory dual, each of the $\alpha_{i}$ would represent some combinations of corrections that are $\sim\alpha^{\prime}$ and corrections that are $\sim\ell_{P}^{2}/L^{2}$, where $\ell_{P}$ is the Planck scale, with $\kappa_{5}^{2}=\ell_{P}^{3}$. While $\alpha^{\prime}$ corrections correspond to inverse coupling corrections in the field theory, $\ell_{P}^{2}/L^{2}$ corrections correspond to $1/N_{c}$ corrections. (See Refs. Grozdanov and Starinets (2017); Grozdanov and van der Schee (2017).) Here, in our bottom-up construction in which we are writing down the most general action to $R^{2}$ order without knowing the identity of its dual field theory, we will treat all the higher-derivative corrections as if they induce coupling constant corrections in the dual field theory. By performing a field redefinition of the metric $g_{\mu\nu}$ to ${\cal O}(\alpha_{i})$ (see e.g. Refs. Brigante et al. (2008); Grozdanov and Starinets (2015)), the action (3) can be brought into the form

[TABLE]

where $\lambda_{\scriptscriptstyle GB}=2\alpha_{3}$. This is the Einstein-Gauss-Bonnet action, which gives rise to second-order equations of motion. In this theory, the coupling $\lambda_{\scriptscriptstyle GB}$ can in principle be treated non-perturbatively as long as it lies in the interval $\lambda_{\scriptscriptstyle GB}\in(-\infty,1/4]$ Brigante et al. (2008); Grozdanov and Starinets (2015). In this work, we will treat $\lambda_{\scriptscriptstyle GB}\sim\alpha^{\prime}$ perturbatively and exploit the two-derivative form of the resulting differential equations to simplify the complexity of numerical simulations. Collisions of sheets of energy density in this theory were studied up to $\mathcal{O}(\lambda_{\scriptscriptstyle GB}^{2})$ in Ref. Grozdanov and van der Schee (2017).

Now, in order to study coupling-dependent collisions at finite $\mu$, we add all possible four-derivative terms involving the gauge field and its coupling to the metric with the exception of Chern-Simons-type terms odd in $A_{\mu}$, which we explicitly set to zero as their presence would violate parity in both the gravitational theory and the field theory to which it is dual. We treat all other higher-derivative couplings as being of the same (perturbative) order as $\lambda_{\scriptscriptstyle GB}$. To find the relevant terms in the bulk action, we start by writing the most general four-derivative action involving $g_{\mu\nu}$ and $A_{\mu}$ Myers et al. (2009); Cremonini et al. (2009, 2010); Grozdanov and Starinets (2017). Then, we perform the field redefinitions of both $g_{\mu\nu}$ and $A_{\mu}$ to first order in higher-derivative couplings and up to two derivatives, and expand the action to fourth order in derivatives. The final form of the four-derivative part of the action has only two remaining coupling constants: $\lambda_{\scriptscriptstyle GB}$ from Eq. (4) and a single new constant $\beta$ that characterizes the strength of the higher-derivative coupling between $g_{\mu\nu}$ and $A_{\mu}$. A particularly convenient choice for the action that will be used in all of our analyses is Grozdanov and Starinets (2017)111Note that if we wanted to include a density of axial charge or to study less symmetric collisions in which vorticity or a magnetic field could arise in the gauge theory plasma, and study consequent anomalous transport effects, we would need to extend the gravitational theory relative to the one that we shall employ in this paper, completing it by adding a Chern-Simons term to this 4+1-dimensional action Gauntlett and Gutowski (2003); Erdmenger et al. (2009); Casalderrey-Solana et al. (2016). We shall not pursue any of these directions, meaning that this term would play no role in any of our analyses. For this reason we leave it out.

[TABLE]

which, in analogy with the Einstein-Gauss-Bonnet theory in Eq. (4), gives rise to purely second-order equations of motion for all components of $g_{\mu\nu}$ and $A_{\mu}$. Since both multiply four-derivative terms in the action, we shall assume throughout that $\beta$ is of the same order as $\lambda_{\scriptscriptstyle GB}$. We furthermore assume that $\mu^{2}/T^{2}$ is perturbatively small and of the order of the $\lambda_{\scriptscriptstyle GB}$ and $\beta$; we shall discuss this further below. We have introduced a control parameter $\epsilon$ as a device for later convenience; it could evidently be scaled away by redefining $A_{\mu}$.

Since $\mu/T$ is small in the QGP produced in ultrarelativistic heavy ion collisions, we shall only consider the case where the incident sheets of energy before the collision have a baryon number density $\rho$ that is much smaller than $\mathcal{E}^{3/4}$, with $\mathcal{E}$ their energy density. Consequently, $\mu/T$ is small in the hydrodynamic fluid produced after the collision. This assumption allows additional simplifications, as we explain below. Small $\rho$ in the boundary theory means small $A_{\mu}$ in the dual Einstein-Maxwell theory, and in particular means that the backreaction of the Maxwell field $A_{\mu}$ on the five-dimensional geometry is small. We can implement this assumption by treating $\epsilon$ as small. In the boundary theory it determines the size of the baryon number density and the chemical potential. In particular we have that $J^{\mu}\propto\epsilon$ and $\mu\propto\epsilon$. To see this, let us for a moment consider the action with the normalization of the Maxwell field used in Eq. (1). Then $\rho$ and $\mu$ are set by the size of the bulk gauge field $|A_{\mu}|$ which has a mass dimension of $1/L$. If we assume that $\rho$ and $\mu$ are small and parametrically ${\cal O}(\epsilon)$, we have that $|A_{\mu}|\sim\epsilon/L$. It is then convenient to explicitly scale $A_{\mu}\rightarrow\epsilon A_{\mu}$ and take $A_{\mu}$ of order $\sim 1/L$, which introduces powers of $\epsilon$ into the action as in Eq. (II). Note that the assumption of the smallness of $A_{\mu}$ in Eq. (1) and consequently the smallness of $\epsilon$ in Eq. (II) is an assumption about the smallness of $\rho$ in the incident sheets of energy and baryon number density. Thus, $\epsilon$ in Eq. (II) is a control parameter governing the initial conditions in our calculation, not a coupling constant in the theory itself. The only coupling constants arising in the higher-derivative terms in Eq. (II) are $\lambda_{\scriptscriptstyle GB}$ and $\beta$.

Introducing $\epsilon$ as we do gives us a convenient way to take the smallness of $\rho$ into account in a combined perturbative expansion: we shall assume that $\epsilon^{2}$ is of order $\lambda_{\scriptscriptstyle GB}$ so that we can arrange our calculations in a single perturbation series in $\lambda_{\scriptscriptstyle GB}$. We will then determine $T^{\mu\nu}$ and $J^{\mu}/\epsilon$ to leading order, ${\cal O}(\lambda_{\scriptscriptstyle GB})$. Our motivation here is simply that we wish to take $\lambda_{\scriptscriptstyle GB}$, $\beta$ and $\epsilon$ (and consequently $\rho$) to all be small and work to lowest nontrivial order in a combined expansion: the particular choice of scaling $\epsilon^{2}\sim|\lambda_{\scriptscriptstyle GB}|$ is arbitrary; others could be investigated. However, if we had instead assumed $\epsilon\sim|\lambda_{\scriptscriptstyle GB}|$ this would not have captured the physics that we aim to elucidate as in that case the Maxwell field would have no backreaction at ${\cal O}(\lambda_{\scriptscriptstyle GB})$, and we would effectively be working at zero chemical potential rather than at small chemical potential. Since in our scheme \epsilon^{2}$$,\lambda_{\scriptscriptstyle GB} and $\beta$ scale in the same way we will frequently express $\epsilon^{2}$ and $\beta$ in terms of $\lambda_{\scriptscriptstyle GB}$ as

[TABLE]

where our choice of scaling means that both $c_{\epsilon}$ and $c_{\beta}$ are ${\cal O}(1)$ when we count powers of $\epsilon$.

With the scaling that we have described, the action Eq. (II) yields the most general equations of motion to linear order in $\lambda_{\scriptscriptstyle GB}$. The Maxwell ($F^{2}$) and Gauss-Bonnet ($R^{2}$) terms are of the same order (since $\epsilon^{2}\sim\lambda_{\scriptscriptstyle GB}$) while the term proportional to $\beta\epsilon^{2}$ ($RF^{2}$) is subleading. However, the factor $\epsilon^{2}$ cancels in the Maxwell equations, so the $\beta$ term provides the leading correction to the Maxwell equations. Further details regarding the construction of the action Eq. (II) and the field redefinitions involved are discussed in Appendix A. We close our discussion here by noting that if we had not assumed that $\epsilon$ (and consequently $\rho$ and $\mu/T$) is small then we would have had to include additional $F^{4}$ terms in the action (II).

In this work, we want to analyze the effects of coupling constant corrections on the dynamics of energy-momentum and baryon density, which are encoded in one-point functions of the energy-momentum tensor $T^{\mu\nu}$ and a current $J^{\mu}$, respectively. Both conserved tensors need to be calculated to first order in our combined expansion in small $\lambda_{\scriptscriptstyle GB}$, $\beta$ and $\epsilon^{2}$. This means that $T^{\mu\nu}$ and $J^{\mu}$ will take the form

[TABLE]

where we note that $\mu^{2}/T^{2}\propto\epsilon^{2}$ and where all of $\mathcal{T}^{\mu\nu}$, $\mathcal{J}^{\mu}$, and all of the $\delta\mathcal{T}^{\mu\nu}$ and $\delta\mathcal{J}^{\mu}$ are zeroth order in small quantities. We shall obtain solutions for the time evolution of these quantities by developing and solving the equations of motion for the dual gravitational Einstein-Maxwell-Gauss-Bonnet theory in Sections III and IV.

We shall also wish to compare the dynamics of $T^{\mu\nu}$ and $J^{\mu}$ that we obtain holographically to what we obtain from hydrodynamics. To this end, we shall need the hydrodynamic expansion of the two conserved tensors to first order in derivatives, which is given by

[TABLE]

where $u^{\mu}$ is the velocity field normalized to $u^{2}=-1$, $\mu$ and $T$ are treated as hydrodynamic near-equilibrium quantities, and $\Delta^{\mu\nu}=u^{\mu}u^{\nu}+\eta^{\mu\nu}$ is the projector in directions transverse to $u^{\mu}$. The thermodynamic quantities appearing in these expressions are the energy density $\mathcal{E}$, pressure $p$ and baryon number density $\rho$. Assuming local equilibrium, they satisfy the thermodynamic relation

[TABLE]

with $s$ the entropy density, and all of them vary in space and time on wavelengths that are long compared to $1/T$. Because ${\cal N}=4$ SYM theory and its deformations considered in this work are conformal, the equation of state is given by $\mathcal{E}=3p$. The two transport coefficients that enter into the hydrodynamic constitutive relations (7) and (8) are the shear viscosity $\eta$ and the baryon number conductivity $\sigma$. (Recall that $J^{\mu}$ is not the electric current, but the baryon current.) However, if $J^{\mu}$ were to be weakly gauged, with the introduction of an additional $U(1)$ gauge field, then $\sigma$ would indeed become the electrical conductivity of the strongly coupled field theory. For further details on relativistic hydrodynamics, see Ref. Kovtun (2012).

We shall compute $\mathcal{E}$ and $p$ to linear order in $\lambda_{\scriptscriptstyle GB}$, meaning to ${\cal O}(\epsilon^{2})$, and in so doing will precisely reproduce the analysis of finite-coupling corrections in the absence of any baryon number density from Ref. Grozdanov and van der Schee (2017). Note that there are no ${\cal O}(\epsilon^{3})$ terms ($\sim\mu^{3}$ or $\sim\lambda_{\scriptscriptstyle GB}\mu$) in $\mathcal{E}$ and $p$, since these thermodynamic quantities are even in $\mu$. We shall compute $\rho$, which at infinite coupling is proportional to $\epsilon$ for small $\epsilon\sim\mu$, to order $\epsilon\lambda_{\scriptscriptstyle GB}\sim\epsilon^{3}$. In terms of the gravity couplings $\lambda_{\scriptscriptstyle GB}$ and $\beta$, only $\lambda_{\scriptscriptstyle GB}$ will influence the energy density $\mathcal{E}$ and pressure $p$, while both $\lambda_{\scriptscriptstyle GB}$ and $\beta$ will affect the baryon number density $\rho$. Working to ${\cal O}(\epsilon\lambda_{\scriptscriptstyle GB})$ in the analysis of $\rho$ is necessary in order for us to find the leading-order finite-coupling corrections to the baryon number density, which is the central focus of this paper. For later convenience, we define the ${\cal O}(1)$ quantities $\bar{\rho}\equiv\rho/\epsilon$, $\bar{\mu}\equiv\mu/\epsilon$ and $\bar{J}^{\mu}\equiv J^{\mu}/\epsilon$. Note that all these quantities are zeroth order in $\epsilon$, meaning in particular that they are finite in the $\epsilon\rightarrow 0$ limit.

Before we turn to a detailed analysis of the action (II) in the next Section, we close this section by motivating a qualitative sense for what the relevant range of choices might be for the four-derivative coupling parameters $\lambda_{\scriptscriptstyle GB}$ and $\beta$ in our gravitational action. This certainly cannot be done in a quantitative way because although we do not know exactly what deformation of ${\cal N}=4$ SYM theory our bottom-up gravitational construction is dual to, we know that it is not dual to QCD, meaning that we cannot apply constraints coming from phenomenology or lattice QCD calculations quantitatively. As noted below Eq. (3), if we had derived the action (II) from a top-down stringy construction we would know precisely which gauge theory it is dual to and in such a construction both $\lambda_{\scriptscriptstyle GB}$ and $\beta$ would have been fixed in terms of $\alpha^{\prime}$ and $\ell_{P}^{2}/L^{2}$. Instead, in our bottom-up approach reasonable values for these quantities can be estimated by asking that $\lambda_{\scriptscriptstyle GB}$- and $\beta$-dependent observables take on reasonable values. We need at least some rough sense of what values of these parameters we should use when we later plot our results, and to this end we shall look at what we know about the shear viscosity to entropy density ratio $\eta/s$ and the baryon number susceptibility $\chi$ in QCD and ask how this compares with their dependence on $\lambda_{\scriptscriptstyle GB}$ and $\beta$ in the model theory that we are employing. Our goal is to set these couplings to values such that the magnitude of their effects on $\eta/s$ and $\chi$ in our model theory corresponds to the magnitude of the effect of the finiteness of the QCD gauge coupling on these observables, in the hope that doing so will allow us to use our model theory to get a qualitative impression of the effect of the finiteness of the gauge coupling on the dynamics of $\rho$ in collisions.

We start from a relatively well-established fact that when $\lambda_{\scriptscriptstyle GB}R^{2}$ terms in the gravitational dual are treated as enabling finite coupling corrections in the gauge theory, it is natural to take $\lambda_{\scriptscriptstyle GB}$ to have a negative sign. (See e.g. Refs. Grozdanov et al. (2016); Grozdanov and van der Schee (2017); Andrade et al. (2017); DiNunno et al. (2017); Atashi et al. (2017).) One way to motivate this sign is to look at the result for $\eta/s$ in (zero baryon number) Einstein-Gauss-Bonnet (4), where it is known nonperturbatively in $\lambda_{\scriptscriptstyle GB}$ and is given by Kats and Petrov (2009); Brigante et al. (2008)

[TABLE]

meaning that if we choose $\lambda_{\scriptscriptstyle GB}$ negative $\eta/s$ increases with increasing $|\lambda_{\scriptscriptstyle GB}|$, meaning with increasing inverse-gauge-coupling, or in other words, with decreasing gauge coupling. This is the natural sign, given that in a weakly coupled gauge theory the leading order dependence of $\eta/s$ on the gauge coupling $g$ at small $g$ is given by $\eta/s\sim{\rm const}/[g^{4}\ln(1/g)]$ Arnold et al. (2000), meaning that at weak coupling $\eta/s$ increases with decreasing gauge coupling and diverges as $g\to 0$. We will typically choose $\lambda_{\scriptscriptstyle GB}=-0.2$, which corresponds to an 80% increase in the value of $\eta/s$ relative to its value when $\lambda_{\scriptscriptstyle GB}=0$ and the gauge theory has infinite coupling. There is nothing sacred about this particular value, in particular given that $\eta/s$ in QCD will in reality have some temperature dependence, but this choice puts $\eta/s$ within the range estimated by comparing hydrodynamic calculations of the anisotropic expansion of the droplets of QGP produced in off-center heavy ion collisions with experimental data Romatschke and Romatschke (2007); Schenke et al. (2011); Bernhard et al. (2016); for reviews, see Refs. Heinz and Snellings (2013); Romatschke and Romatschke (2019); Busza et al. (2018).

At small but nonzero baryon number density, perturbative results in both the four-derivative couplings found in the Maxwell-Einstein-Gauss-Bonnet action (II) give rise to corrections to the ratio of $\eta/s$ that take the form

[TABLE]

the ellipsis denotes terms of $\sim{\cal O}(\lambda_{\scriptscriptstyle GB}^{3})$ and beyond. Sufficiently little is known about the $\mu$-dependence of $\eta/s$ in QCD that we will not attempt to use a calculation of the $\mu$-dependent contributions to (13) to estimate what range of $\beta$ we might investigate. Another reason why attempting this would be perilous is that the presently unknown contributions to $\eta/s$ that are of order $\lambda_{\scriptscriptstyle GB}^{2}$ are parametrically just as significant, given the scaling that we are using, as the known $\epsilon^{2}\beta$ and $\epsilon^{2}\lambda_{\scriptscriptstyle GB}$ terms, both of which can be computed from the field redefinitions discussed in Appendix A and the calculations of Refs. Myers et al. (2009); Cremonini et al. (2009).

We shall estimate a reasonable range for $\beta$ by looking at its effects on the baryon number susceptibility $\chi$. Note that in our model theory, $\chi\propto N_{c}^{2}$ because all degrees of freedom are in color-adjoint representations whereas in QCD, $\chi\propto N_{c}N_{f}$ because baryon number is carried by quarks, which come in $N_{c}$ colors and $N_{f}$ flavors. We can scale out this difference, however, by computing the ratio of $\chi$ in the theory with finite coupling to $\chi$ in the non-interacting limit in both QCD and our model theory, since the count of the number of degrees of freedom that carry baryon number cancels in this ratio. Beginning with our model theory, we note that we will show in Sec. III (see Eq. (60)) that to linear order in $\mu$ and to lowest order in $\lambda_{\scriptscriptstyle GB}$ and $\beta$ the baryon number density is given (for a choice of normalization, see Section III) by

[TABLE]

giving a susceptibility

[TABLE]

From Refs. Son and Starinets (2006); Teaney (2006); Casalderrey-Solana et al. (2014a) we know that in $\mathcal{N}=4$ SYM theory we have $\chi_{\infty}/\chi_{0}=\frac{1}{2}$ where $\chi_{\infty}$ is the susceptibility in the infinite coupling limit (corresponding to $\lambda_{\scriptscriptstyle GB}=\beta=0$ in the above), while $\chi_{0}$ is the susceptibility in the non-interacting limit. This gives $\chi_{0}=\pi^{2}T^{2}$ in our normalization scheme. Now we can calculate the ratio $\chi/\chi_{0}$, namely the ratio between $\chi$ at the intermediate coupling that corresponds to given values of $\lambda_{\scriptscriptstyle GB}$ and $\beta$ to $\chi$ in the noninteracting limit, obtaining

[TABLE]

This can be compared to lattice calculations of the same ratio in QCD Borsányi et al. (2012), which indicate that \frac{\chi}{\chi_{0}}\big{|}_{\rm QCD}\in[0.8,0.9] for temperatures from $250\text{\,}\mathrm{MeV}400\text{\,}\mathrm{MeV}$, corresponding to the temperatures of the strongly coupled QGP liquid produced in heavy ion collisions at RHIC and the LHC. Setting $\lambda_{\scriptscriptstyle GB}=-0.2$ and varying \frac{\chi}{\chi_{0}}\big{|}_{\rm QCD} in the given range we find $\beta\in[0.08,0.09]$. If we further vary $\lambda_{\scriptscriptstyle GB}$ in the broad range $[-\frac{1}{2},0]$, corresponding to values of $\frac{\eta}{s}$ between $\frac{1}{4\pi}$ and $\frac{3}{4\pi}$, we find that $\beta$ should lie in the broad range

[TABLE]

That said, we caution that working only to linear order in $\beta$ may not be a good approximation, in particular at the upper end of this range. Working only to this order we cannot reliably estimate the range of reliability of this linear approximation, but absent further information it is reasonable to guess from results like (16) that it starts to break down once $16\beta$ gets as large as $1$. Another way of making a reasonable guess is to note from (16) that values of $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0.1$ imply that $\rho$ increases by a factor 2, suggesting that the linear approximation is breaking down. And, if we express $\mu$ as a function of $\rho$ and $T$ we see that to the order we are working it is proportional to $1-3\lambda_{\scriptscriptstyle GB}-16\beta$, which for these particular values would be zero and hence clearly outside the linear regime. In the results that we shall present in Section IV, we will see that values of $\beta$ in the range (17) can have a large effect on the baryon number distribution, in some cases to a degree that supports our guess that working to linear order in $\beta$ may not be sufficient.

Summarizing these considerations, it seems reasonable to investigate the consequences of choosing $\lambda_{\scriptscriptstyle GB}$ small and negative and $\beta$ small and positive in the action (II) on the dynamics of $\rho$ in collisions. We shall present results in Section IV with this in mind, but first, in the next Section, we must describe how we analyze the dynamics governed by the action (II).

III Analysis of the model and details of the numerical setup

In this Section we proceed with the analysis of the action (II). We shall present the equations of motion, describe the initial conditions that we employ and the numerical methods that we use to solve for the time evolution, provide the dictionary that connects bulk quantities obtained via solving these holographic evolution equations to field theory observables, and close by remarking upon the hydrodynamics and thermodynamics of the model.

The action (II) is the most general 4+1-dimensional four-derivative theory containing an asymptotically AdS metric and a gauge field which is consistent with working to first order in $\lambda_{\scriptscriptstyle GB}$ and assuming the scaling $\lambda_{\scriptscriptstyle GB}\sim\beta\sim\epsilon^{2}\sim\mu^{2}/T^{2}$. The equations of motion that follow from variations of $g_{\mu\nu}$ and $A_{\mu}$ were derived in Appendix E of Ref. Grozdanov and Starinets (2017). The resulting Einstein equations are

[TABLE]

where the perturbatively small terms on the right-hand-side of the equation are

[TABLE]

The Maxwell equations, as modified by the four-derivative terms, are

[TABLE]

Both the equations of motion (18) and (22) contain at most second derivatives with respect to each of the coordinates. The fact that, with our scaling, $T^{\beta}_{\mu\nu}\sim\mathcal{O}(\epsilon^{2}\beta)\sim{\cal O}(\lambda_{\scriptscriptstyle GB}^{2})$ implies that to leading order in $\lambda_{\scriptscriptstyle GB}$ the $\beta$-term in the action only affects the Maxwell equations (22). We shall employ the convention

[TABLE]

This can be thought of as choosing our units, although looked at that way it looks unconventional. We make this choice for later convenience because with this choice the non-normalizable modes of the metric near the boundary are unaffected by our $\lambda_{\scriptscriptstyle GB}$ corrections.

Equations (18) and (22) have a solution of the following form:

[TABLE]

where $x_{\pm}=t\pm z$ are lightcone coordinates and where the functions $h_{i}(x_{+})$ and $a_{+(i)}(x_{+})$ can be chosen arbitrarily. This metric describes a gravitational wave moving at the speed of light in the $z$-direction in an asymptotically AdS spacetime. On the boundary it represents a sheet of baryon and energy density moving at the speed of light. A wave moving in the opposite direction is obtained by exchanging $x_{+}\leftrightarrow x_{-}$. The wave solution is planar (and translation-invariant) in the two-dimensional $\mathbf{x}_{\perp}$ plane. The functions $h_{i}(x_{+})$ and $a_{+(i)}(x_{+})$ determine the profile of the energy density and baryon density, respectively, along the “beam” direction. In our computation, we use a pair of well-separated left- and right-moving wave solutions, namely (24)–(25) and their counterparts with $x_{+}\leftrightarrow x_{-}$, as initial conditions. Specifying our initial conditions therefore requires choosing the eight functions $h_{i}(x_{\pm})$ and $a_{\pm(i)}(x_{\pm})$. We choose the four profile functions $h_{0}(x_{\pm})$ and $a_{+(0)}(x_{\pm})$ all to be Gaussians with the same width $w$:

[TABLE]

where with this choice $m^{3}$ is the energy density per unit transverse area of each incident sheet of energy. We then choose the four first order functions $h_{1}(x_{\pm})$ and $a_{\pm(1)}(x_{\pm})$ needed to complete the specification of our initial conditions such that the total ${\cal O}(\lambda_{\scriptscriptstyle GB})$ correction to the expectation value of the stress-energy tensor and the conserved $U(1)$ current vanishes in our initial conditions, so that our initial conditions are independent of $\beta$ and $\lambda_{\scriptscriptstyle GB}$. This will allow us to compare simulations with different $\beta$ and $\lambda_{\scriptscriptstyle GB}$ properly. The parameter $w$ sets the width of the sheets of energy and baryon number density. We shall show results for two different choices of $w$ which we shall refer to as thick and thin sheets of energy density. The two values of $w$ at which we shall quote results are $mw=0.1$ and $mw=1.5$, for thin and thick sheets respectively.

In our numerical solution of the equations of motion, we use the ingoing Eddington-Finkelstein coordinates Chesler and Yaffe (2011) in which the metric takes the form

[TABLE]

where the functions $C$, $F$, $S$, $B$ depend on the boundary time $t$, the radial holographic coordinate $r$, and the coordinate $z$ that lies along the direction of motion of the incident sheets of energy density, the “beam direction”. The two coordinates collectively denoted as $\mathbf{x}_{\perp}$ are perpendicular to $z$ in spatial $\mathbb{R}^{3}$. The metric (27) is based on assuming translational symmetry in this plane, appropriate for our analysis of the collision of sheets that are infinite in extent and translationally invariant in these transverse directions. The conformal boundary is at $r=\infty$. For the gauge field, we choose the radial gauge, meaning that its form reads

[TABLE]

In the numerical evolution, we work directly with the field strength $F_{\mu\nu}$. The coordinates used in Eq. (27) allow us to rewrite the full set of Einstein and Maxwell equations as a nested set of ordinary differential equations (ODEs) using the characteristic formulation of general relativity Chesler and Yaffe (2011); Bondi (1960). Details and examples of the numerical procedure in the (unperturbed) characteristic formulation can be found in a number of past publications (see e.g. Ref. Chesler and Yaffe (2014); van der Schee (2014); Attems et al. (2017); Waeber et al. (2019)). These include in particular the (tricky) coordinate transformation from (24) to (27), which luckily is straightforward to generalize to our $\lambda_{\scriptscriptstyle GB}$-corrected set-up (see also Grozdanov and van der Schee (2017)).

For our specific case of Gauss-Bonnet gravity with a Maxwell field, several aspects are noteworthy. To accommodate a perturbative expansion in $\lambda_{\scriptscriptstyle GB}$ (with our assumed scaling) we expand the metric functions, the gauge field, and later also various thermodynamic quantities as power series in $\lambda_{\scriptscriptstyle GB}$, namely $f=\sum_{n}f_{(n)}\lambda_{\scriptscriptstyle GB}^{n}$ (cf. Eq. (6)). For the metric this implies

[TABLE]

and similarly for all other functions. At every time step in the numerical evolution (solving Eqs. (18) and (22)), we then split the computation into two steps. First, we evaluate the unperturbed functions $\{C_{(0)},F_{(0)},\ldots\}$ and their derivatives according to the nested scheme involving a set of homogeneous differential equations. Second, we solve a nested set of nonhomogeneous differential equations for the first-order perturbations $\{C_{(1)},F_{(1)},\ldots\}$ and their derivatives. The same second-order differential operators act on the perturbations $\{C_{(1)},F_{(1)},\ldots\}$ as in the first step on $\{C_{(0)},F_{(0)},\ldots\}$. The nonhomogeneous source terms needed to solve for $\{C_{(1)},F_{(1)},\ldots\}$ are obtained using the zeroth-order solution computed in the first step. The only subtlety in the procedure involves the computation of the following time derivatives: $\partial_{t}^{2}B_{(0)}$, $\partial_{t}F_{(0)rt}$ and $\partial_{t}A_{0}$. These specific time derivatives do not arise in the differential operators of the homogeneous differential equations, meaning that they are not needed in the nested scheme for evolving the metric and gauge field, or in the computation of the constraints. They do, however, appear in the source terms in the nonhomogeneous equations for $\{C_{(1)},F_{(1)},\ldots\}$. This means that they affect the perturbed first-order solutions, which means that we need to keep track of $\partial_{t}B_{(0)}$, $F_{(0)rt}$ and $A_{0}$ at every step of the evolution and evaluate their time derivatives. Finally, we note that we need not compute the correction to the location of the apparent horizon: since we work perturbatively, the code will be stable if the horizon is located at $r=1+\mathcal{O}(\lambda_{\scriptscriptstyle GB})$. (Note, however, that these corrections were computed in Ref. Grozdanov and van der Schee (2017) in order to see the effects of $\lambda_{\scriptscriptstyle GB}$ corrections on the entropy density.) Of course, the equations generated by the source terms are much longer than the unperturbed equations, which means that the code runs roughly a factor ten times slower. The calculation of the time evolution of the collision of thin sheets of energy density then takes about two weeks to perform on an ordinary desktop computer.

To obtain the boundary field theory quantities that we are interested in, namely the field theory energy-momentum tensor and the baryon number current, from the five-dimensional bulk metric and gauge field functions, we shall require the near-boundary solutions of (18) and (22), which take the form

[TABLE]

The functions of $t$ and $z$ that appear here are the normalizable modes of the gauge field and metric near the AdS boundary; they depend upon the bulk dynamics and therefore must be found from the evolution of the full metric, as we have described above. In order to determine the field theory quantities of interest what we will need to extract from our evolution of the full metric are the functions of $t$ and $z$ appearing on the RHS of (30), (31), (32) and (33). To make this extraction easier, when we solve for the time evolution we recast the equations for $A_{\mu}$, $C$, $B$ and $F$ as equations for $r^{2}A_{\mu}$, $r^{2}(C-r^{2})$, $r^{4}B$ and $r^{2}F$, and solve them numerically in that form van der Schee (2014).

The main quantities in the field theory whose dynamics we wish to determine and study are the one-point functions of the boundary energy-momentum tensor $T^{\mu\nu}$ and the conserved $U(1)$ (baryon number) current $J^{\mu}$. We will work with the rescaled quantities

[TABLE]

The normalization factor of $\kappa_{5}^{2}/2L_{0}^{3}$ is given by $2\pi^{2}/N_{c}^{2}$ in $\mathcal{N}=4$ SYM theory at infinite coupling, i.e. when $\lambda_{\scriptscriptstyle GB}=\beta=0$. $L_{0}$ is the $\lambda_{\scriptscriptstyle GB}=0$ value of the AdS radius, which for convenience we have set equal to 1, see (23) Grozdanov and van der Schee (2017). A particularly straightforward way of ensuring that these normalization conventions are satisfied is to set $\kappa_{5}^{2}/8\pi=1/4\pi$ in all formulae. We will do so below.

The relationships that determine the energy-momentum tensor for the boundary gauge theory from the functions appearing in Eqs. (31)–(33) that are determined by the time evolution in the bulk gravitational theory can be obtained to ${\cal O}(\lambda_{\scriptscriptstyle GB})$ in the same way as they were previously obtained in Einstein-Gauss-Bonnet theory without a bulk gauge field, and take the form Brihaye and Radu (2008); Grozdanov and Starinets (2017)

[TABLE]

with

[TABLE]

We must in addition find the relationship that determines the conserved $U(1)$ (baryon number) current in the boundary field theory from the functions appearing in Eq. (30). The conserved boundary current can be identified as the boundary value of the canonical momentum with respect to $r$ of the bulk Maxwell field. We see this as follows. The gauge field equations of motion in the bulk imply that

[TABLE]

is the canonical momentum conjugate to the bulk gauge field $A_{\nu}$, satisfying $\Pi^{\mu\nu}=-\Pi^{\nu\mu}$, and with $\mathcal{L}$ being the Lagrangian corresponding to the action (II). The covariant divergence of an antisymmetric tensor satisfies $\nabla_{\mu}\Pi^{\nu\mu}=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}\Pi^{\nu\mu})$. Hence,

[TABLE]

where, to avoid confusion, in this paragraph alone we have found it necessary to use Latin indices $a,b,\ldots$ for the indices that run over the four-dimensional boundary coordinates only. From (48) we see that $\sqrt{-g}\,\Pi^{ra}$ is a conserved current within each constant-$r$ slice. We can therefore define the conserved boundary current as

[TABLE]

Inserting the near-boundary expansion of the gauge field (30) and the metric functions (31)–(33) into this expression, we find

[TABLE]

with

[TABLE]

which is the relationship that we need. We find no divergences in the $r\rightarrow\infty$ limit of $\sqrt{-g}\,\Pi^{ra}$, so no counterterms are needed. Note also that we have computed $\Pi^{\mu\nu}$ by taking the derivative with respect to $\epsilon A_{\mu}$, so that $\mathcal{J}\sim O(\epsilon)$, which is equivalent to first calculating $\Pi^{\mu\nu}$ and then recaling $A_{\mu}\rightarrow\epsilon A_{\mu}$, as described in the previous Section. We compute the baryon number density from $\rho=\mathcal{J}^{0}$. Finally, in (50) and (51) we can safely replace the Latin index $a$ by a Greek index $\mu$, returning to the notation that we have used elsewhere.

Having described the equations of motion, the initial conditions and the scheme for computing the time evolution, we can now perform collisions between two incident sheets of energy and baryon number and extract the resulting $\mathcal{T}^{\mu\nu}$ and $\mathcal{J}^{\mu}$. We shall present results in the next Section.

We anticipate that the matter produced in such collisions will rapidly hydrodynamize, becoming strongly coupled liquid plasma whose subsequent expansion and cooling dynamics is well-described by relativistic viscous hydrodynamics, including a conserved (baryon number) current. In order to check that the plasma does indeed hydrodynamize, after presenting our results we shall need to define the standard hydrodynamic variables (including the baryon number density) in terms of $\mathcal{T}^{\mu\nu}$ and $\mathcal{J}^{\mu}$ and test the validity of the hydrodynamic approximation (9)–(10) by comparing solutions to these hydrodynamic equations to the full results that we shall obtain by solving the holographic equations of motion that determine the functions defined in (30)–(33) and appearing in the energy-momentum tensor (36) and the current (50).

In order to make a comparison to hydrodynamics, we need to obtain the local fluid velocity $u^{\mu}$ and then the local energy density $\mathcal{E}_{\rm loc}$ and the local baryon number density $\rho_{\rm loc}$ defined in the rest frame of the fluid. We proceed as follows. At a boundary point $x$ of the bulk spacetime, we define $u_{\mu}$ as the timelike eigenvector of $-\mathcal{T}\indices{{}^{\mu}_{\nu}}(x)$, hence,

[TABLE]

The local baryon number density is then

[TABLE]

To evaluate (9) and (10) we also need the equilibrium expressions for $T$ and $\mu$ as functions of $\mathcal{E}$ and $\rho$. (We shall use these expressions to obtain the local $T$ and $\mu$ from $\mathcal{E}_{\rm loc}$ and $\rho_{\rm loc}$.) We obtain the relationships that we need from the black brane solution of the bulk theory gravitational equations of motion, as this is the dual of the equilibrium state in the boundary field theory. To linear order in $\lambda_{\scriptscriptstyle GB}$, the black brane solution is given by

[TABLE]

where $g_{0}$, $g_{1}$, $k_{0}$, $k_{1}$, $Q_{0}$ and $Q_{1}$ are free parameters. The Hawking temperature $T$, the chemical potential $\mu$, and the entropy density $s$ of the black brane solution are given by

[TABLE]

where the horizon radius $r_{+}\equiv r_{0}+\lambda_{\scriptscriptstyle GB}r_{1}+{\cal O}(\lambda_{\scriptscriptstyle GB}^{2})$ is defined as the largest solution to $C(r)=0$ and where $\Sigma$ is the black brane horizon area. Note that the entropy density of a black brane in Gauss-Bonnet gravity satisfies the same Bekenstein-Hawking formula as in pure Einstein theory Cai (2002), and, as we discuss in Appendix A we find that the same is true when we add the Maxwell and Maxwell-curvature coupling terms. Recall, however, that what we need is expressions that relate $T$ and $\mu$ to $\mathcal{E}$ and $\rho$. We can use (41), (46) and (51) to express the energy density $\mathcal{E}$, pressure $p$ and baryon number density $\rho$ as functions of $k_{i}$, $Q_{i}$ and $r_{i}$. Then, we can solve equations (56)–(58) together with $C(r_{+})=0$ for $k_{i}$, $Q_{i}$ and $r_{i}$ in terms of $T$ and $\bar{\mu}\equiv\mu/\epsilon$. Plugging this into our expressions for $\mathcal{E}$, $p$ and $\rho$ we find

[TABLE]

where we remind that $\bar{\rho}\equiv\rho/\epsilon$. Finally, we can solve these expressions for $\bar{\mu}$ and $T$ which yields

[TABLE]

The expressions (61) and (62) are the expressions that we shall use in order to obtain the local $T$ and $\mu$ from $\mathcal{E}_{\rm loc}$ and $\rho_{\rm loc}$, so that we can then use (52) and (53) in the evaluation of the expressions (9) and (10) for the two conserved tensors in the hydrodynamic description of the fluid.

In order to complete the evaluation of (9) and (10) we also need the transport coefficients $\eta/s$ and $\sigma/T$. We discussed the first of these already in the previous Section; to the order at which we are working, it is given by

[TABLE]

The conductivity can be obtained by using the relation $\sigma=D\chi$, where $\chi$ is the baryon susceptibility and $D$ the baryon diffusion constant. $D$ was calculated nonperturbatively in both couplings at $\mu=0$ in our theory in Ref. Grozdanov and Starinets (2017) and perturbatively at $\mu\neq 0$ in a theory related to ours by field redefinitions in Ref. Myers et al. (2009). Appendix A contains a discussion on how to convert the result of Ref. Myers et al. (2009) to our theory. The result, with the conventions given in Eqs. (23) and (35), is

[TABLE]

We shall take $\mathcal{E}_{\rm loc}$ and $\rho_{\rm loc}$ from our holographic simulation, use (61) and (62) to obtain $T$ and $\mu$, and then use (52) and (53) to evaluate (9) and (10). From these hydrodynamic solutions we can obtain hydrodynamic “predictions” for the time evolution of the longitudinal and transverse pressure as well as for the baryon number density and current, all of which we can compare to our full results for the evolution of these quantities, obtained from our holographic calculation. This comparison will allow us to assess the degree to which the matter produced in these collisions does in fact hydrodynamize.

As a final consistency check we have verified that thermodynamic relation

[TABLE]

is indeed satisfied. Note that the $\beta$-dependent correction to $\mathcal{J}^{\mu}$ in Eq. (51) contributes to (65) only at ${\cal O}(\lambda_{\scriptscriptstyle GB}^{2})$, since $\rho$ and $\mu$ are ${\cal O}(\epsilon)$, so that the leading $\beta$-independent contribution of $\mu\rho$ is already ${\cal O}(\lambda_{\scriptscriptstyle GB})$. We have computed a black brane solution and an expression for $\mathcal{T}^{\mu\nu}$ also to ${\cal O}(\lambda_{\scriptscriptstyle GB}^{2})$ and verified that (65) holds to second order as well.

IV Results

In this section we present the results of our numerical calculations of the collisions of thick and thin sheets of energy and baryon number density, as described in Section III.

IV.1 Colliding thick and thin sheets

We begin the discussion by considering the case with $\epsilon=0$. This means that we compute the stress-energy and baryon number density while neglecting the back-reaction of the baryon number density on the energy-momentum tensor. In equilibrium terms, this limit corresponds to the situation in which the ${\cal O}\left(\frac{\mu^{2}}{T^{2}}\right)$ corrections to $\bar{\rho}$ and $\mathcal{E}$ are neglected. On the bulk side of the duality, setting $\epsilon=0$ means that we neglect the backreaction of the gauge field dynamics on the Einstein equations (18). The evolution of energy and momentum in the $\epsilon=0$ system is therefore completely equivalent to what was studied in Ref. Grozdanov and van der Schee (2017) in Einstein-Gauss-Bonnet theory. Note that $\epsilon$ cancels from the Maxwell equation (22), which implies that the dynamics of the baryon number density whose time evolution we shall follow remains (highly) nontrivial even at $\epsilon=0$.

To begin with a view of our results that is both illustrative and instructive, we show in Figs. 1 and 2 the evolution of energy and baryon number densities with $\epsilon=0$ and $\lambda_{\scriptscriptstyle GB}=-0.2$ for thick and thin sheets, respectively. To the order that we are working, the energy density has no dependence on $\beta$. We plot the baryon number density for two values of $\beta$ in each case, which we choose to be $\beta\in\{0,0.1\}$ for thick sheets and $\beta\in\{0,0.025\}$ for thin sheets. We choose the (smaller) values of $\beta$ that we use for thin sheets in order to be safely within the range of values of $\beta$ where all the qualitative considerations that we discussed at the end of Section II indicate that working to linear order in $\beta$ is likely a safe approximation. And indeed, if we increase $\beta$ further to $\beta=0.05$ we find that after the collision of thin sheets the baryon number on the lightcone increases for a time, which seems unphysical and thus suggests that working to linear order in $\beta$ is inadequate. That said, in the case of collisions of thick sheets we find that $\beta$’s in this range have only modest effects and so for this case we have been more bold and extended our calculations up to $\beta=0.1$. It will take future higher-order calculations to reliably estimate up to what $\beta$ our linear approximation is under control.

Consistent with our general expectations about the phenomenology of heavy ion collisions, after dialing the coupling down from infinity towards intermediate values we find that the inverse coupling constant corrections that we have computed do indeed cause the collisions to become more transparent to both energy density and baryon number density in the sense that a larger amount of both energy and baryon number ends up moving close to the light-cone after the collision at reduced coupling. With $\epsilon=0$ our results for the energy density are the same as in Ref. Grozdanov and van der Schee (2017), so we refer the reader to that work for a full analysis of the results shown in the upper panels of Figs. 1 and 2. We shall focus on the baryon number dynamics.

IV.2 Effect of $\beta$ and $\lambda_{\scriptscriptstyle GB}$ on baryon number dynamics

In this work, we are primarily concerned with the analysis of the effects of $\beta$ and $\lambda_{\scriptscriptstyle GB}$ on the dynamics of baryon number density. In Figs. 3 and 4 we show snapshots that capture the evolution of the baryon number density profiles after the collisions (for $t\geq 0$) of thick and thin sheets, respectively. For comparison, at each time step we show a curve (solid) describing the baryon number dynamics in the full theory, with $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta$ nonzero (0.1 for the thick sheets; 0.025 for the thin sheets) as well as a thin-dashed curve for the theory with $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0$ (whose gravitational dual is Einstein-Maxwell-Gauss-Bonnet theory) and a thick-dashed curve for the theory with $\lambda_{\scriptscriptstyle GB}=\beta=0$ (whose gravitational dual is Einstein-Maxwell theory). The snapshots at different times are denoted by different colors.

In the case of the thick sheets, we see that turning on only $\lambda_{\scriptscriptstyle GB}=-0.2$ with $\beta=0$ slightly broadens the baryon number density profile at all times, indicating a somewhat increased transparency, similar to what is seen for the energy density Grozdanov and van der Schee (2017). Even at late times, as depicted in Fig. 5, the effect is small and the qualitative picture remains unchanged. We see that for both zero and nonzero $\lambda_{\scriptscriptstyle GB}$ when $\beta=0$, the thick sheets come to a full stop after which point the plasma expands.

When we turn on a nonzero $\beta$ (specifically $\beta=0.1$), we see in Figs. 3 and 5 that a small bump of baryon number appears on the lightcone, indicating that not all of the baryon number is fully stopped in the collision. In Fig. 5 we can see this bump develop as we compare the baryon number density at late times ($t=13/m$, the latest time plotted in Fig. 3) in collisions of thick sheets with increasing values of $\beta$ at a fixed value of $\lambda_{\scriptscriptstyle GB}$.

For the thin sheets at infinite coupling, we see in Fig. 4 that a significant amount of baryon number remains on the lightcone after the collision, and baryon number is gradually deposited in the plasma. The change in the baryon number density profile upon setting $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0$ is modest, while for the case of $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0.025$ (not $\beta=0.1$ as for thick sheets), a significant increase of baryon number is seen on the lightcone. This is clear from Fig. 4, where at the latest plotted times, the peak baryon number density near the lightcone is more than doubled compared to the infinitely coupled limit. In Fig. 6, we highlight the $\beta$-dependence of this increased transparency by focusing on the $\beta$-dependence of the baryon number density distribution in collisions of thin sheets at $t=5/m$, the latest time plotted in Fig. 4. Increasing $\beta$ to $\beta=0.05$ further increases the effect. We find that if we were to set $\beta\approx 0.3$, which is far too large for the linear approximation that we are using to be reliable, the baryon number density around $z=0$ is reduced to the point that it vanishes at the late time shown in Fig. 6. (In the case of the collision of thick sheets, as in Fig. 3, our calculation would indicate that the baryon number density around $z=0$ is reduced to near vanishing at late times after the collision for $\beta\approx 1$, which is even farther beyond the regime where the linear approximation that we are using can be trusted.)

To make this discussion more precise, we define a quantity $R_{\rm LC}(t)$ to be the fraction of the total baryon number at a given time that is found within $2w$ (with $w$ being the width of each incident sheet) of one lightcone or the other:

[TABLE]

where $I(t)\equiv[-t-2w,-t+2w]\cup[t-2w,t+2w]$. We find that at $t=13/m$ after the collision of thick sheets in Figs. 3 and 5, $R_{\rm LC}^{\rm(thick)}\in\{0.21,0.22,0.28\}$ for $(\lambda_{\scriptscriptstyle GB},\beta)\in\{(0,0),(-0.2,0),(-0.2,0.1)\}$. And, we find that at $t=5/m$ after the collisions of thin sheets in Figs. 4 and 6, $R_{\rm LC}^{\rm(thin)}\in\{0.16,0.20,0.36\}$ for $(\lambda_{\scriptscriptstyle GB},\beta)\in\{(0,0),(-0.2,0),(-0.2,0.025)\}$. We note that increasing $\beta$ from 0 to 0.1 at fixed $\lambda_{\scriptscriptstyle GB}=-0.2$ increases $R_{\rm LC}^{\rm thick}$ from 0.22 to 0.28 while increasing $\beta$ from 0 to 0.025 at fixed $\lambda_{\scriptscriptstyle GB}=-0.2$ increases $R_{\rm LC}^{\rm thin}$ from 0.20 to 0.36. In both cases, increasing $\beta$ increases the transparency of the collisions to baryon number density, with a higher fraction of the total baryon number ending up near the lightcones. And, we see that the relative importance of $\beta$ increases with decreasing sheet width: the increase in transparency becomes more sensitive to $\beta$ as the sheet width decreases.

In fact, the first-order perturbative equations that describe the time evolution and that we have solved to obtain these results are linear in both $\lambda_{\scriptscriptstyle GB}$ and $\beta$. This means that the baryon number density $\rho$, at fixed $(z,t)$, is a linear function of $\beta$ and $\lambda_{\scriptscriptstyle GB}$ (and of $\epsilon^{2}$, which has so far been set to zero). It follows that $R_{\rm LC}(t)$ is also linear in the same variables, and we can define the coefficients $R_{i}$ as

[TABLE]

In Fig. 7 we show the time evolution for the coefficients $R_{0}$, $R_{\beta}$ and $R_{\lambda_{\scriptscriptstyle GB}}$. Indeed, we observe that the $\beta$ term is responsible for the dominant contribution to the late time baryon number on the lightcone and that the $\beta$ term has a much greater effect for thin sheets. In particular, at the latest times plotted in Fig. 7 (which are the times at which Figs. 5 and 6 are plotted), we find that

[TABLE]

Comparing Eqs. (68) and (69) demonstrates that the effects of $\lambda_{\scriptscriptstyle GB}$ and $\beta$ depend on the width of the sheets. Their influence on the dynamics of baryon number transport and hence on its distribution in the final state is substantial in the collision of thin sheets, and becomes less for thicker sheets. For example, if we set $\lambda_{\scriptscriptstyle GB}=-0.2$ and compare $R_{\rm LC}$ with $\beta=0.05$ to $R_{\rm LC}$ with $\beta=0$, for thick sheets $R_{\rm LC}$ only increases from 22% to 25% while for thin sheets it increases from 20% to 51%.

From the point of view of the gravitational action (II), it is natural to ascribe this behavior to the fact that the two bulk couplings $\lambda_{\scriptscriptstyle GB}$ and $\beta$ which represent inverse-coupling corrections in the gauge theory multiply higher-derivative terms in the gravitational action. As such, when derivatives of the fields are larger (in the gauge theory and consequently in the bulk description) the consequences of the terms in the action multiplied by these couplings should become larger, larger relative to results obtained at infinite coupling in the gauge theory when the gravitational action is the two-derivative Einstein-Maxwell action. This is consistent with what we have found by direct calculation of the time evolution. What is particularly interesting is that the effect of $\beta$ on the dynamics is much more strongly dependent on the width of the sheets than that of $\lambda_{\scriptscriptstyle GB}$. The effect of $\beta$ is illustrated in Figs. 5 and 6. Since $\mathcal{E}$ and $\rho$ are linear in $\lambda_{\scriptscriptstyle GB}$, $\beta$ and $\epsilon^{2}$, the trends seen upon varying $\beta$ in these Figures are also obtained for any other values of $\lambda_{\scriptscriptstyle GB}$ and $\epsilon^{2}$, when working to linear order.

It is interesting to speculate on the qualitative implications of these results for QCD. The collision of thick sheets can be seen as a model for the collision of less ultrarelativistic, and therefore less Lorentz-contracted, heavy ions, with the thin sheets then modelling more ultrarelativistic heavy ion collisions, with higher beam energy. We now recall from our discussion in the Introduction that as the beam energy of a heavy ion collision is increased the baryon transparency increases: the baryon number density deposited in the QGP produced at mid-rapidity decreases and more of the baryon number ends up at higher rapidities, closer and closer to the lightcone. This means that as the beam energy of heavy ion collisions is increased, the baryon number density in these QCD collisions becomes more and more different from what is seen in holographic collisions at infinite coupling. In our calculation, increased $\beta$ corresponds to studying holographic collisions with a reduced gauge coupling. It is therefore tempting to speculate that the “reason” that a given value of $\beta$ makes a much bigger difference in our calculations (causing a much bigger increase in baryon transparency) for collisions of thin sheets than for collisions of thick sheets is that the infinite-coupling results for thin sheets are “more wrong” than those for thick sheets, and hence receive larger finite-coupling corrections.

We close our discussion of these results by noting that the values of $\beta$ that we have employed are around what we described as the reasonable range, (17), based upon consideration of the effect of $\beta$ on the baryon susceptibility in equilibrium. But, at the same time, these values of $\beta$ are in a regime where, as we expected based upon the discussion at the end of Section II, we find $\beta$-induced effects that are so large that the linear approximation is near or perhaps somewhat beyond its regime of validity. Given how different the baryon number distribution is in ultrarelativistic heavy ion collisions in QCD relative to that obtained in holographic calculations at infinite coupling, it is quite pleasing that reasonable values of $\beta$ have a substantial effect, in particular in the case of collisions of thin sheets. It is even more pleasing that these effects go in the right direction, increasing the baryon transparency and leading to more baryon number ending up near the lightcones after the collision.

A necessary consequence of these pleasing conclusions, however, is that further corrections that are higher order in $\beta$ (as well as all corrections that come in up to the same higher order in the inverse-gauge-coupling including those that are higher order in $\lambda_{\scriptscriptstyle GB}$) are likely important. Cranking up $\beta$ increases the baryon number on the lightcones and if we work only to linear order this increase is linear in $\beta$. This increase must be tamed at large enough $\beta$ by corrections that are higher order in $\beta$. We leave the calculation of higher-order-in-inverse-coupling corrections to future work.

IV.3 Spacetime rapidity distribution

Next, to further illustrate our results we transform from Minkowski coordinates $(t,z)$ to proper time $\tau=\sqrt{t^{2}-z^{2}}$ and spacetime rapidity $y={\rm arctanh}(z/t)$, compute the baryon number density in the local fluid rest frame $\rho_{\rm loc}$, and plot the spacetime rapidity distribution of $\rho_{\rm loc}$, shown in Fig. 8. For collisions of both thick and thin sheets, we see that reducing the gauge coupling reduces the baryon number density in the local fluid rest frame at mid-rapidity and broadens the rapidity distribution. This is similar to the effect of $\lambda_{\scriptscriptstyle GB}$ on the rapidity profile of the energy density analyzed in Ref. Grozdanov and van der Schee (2017). For the baryon number, however, there are in addition effects introduced by the $\beta$-dependent terms. And, as we saw before, the consequences of reasonable nonzero values of $\beta$ can be large. In particular, at $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0.025$, at a proper time $\tau=2.5/m$ after the collision of thin sheets the baryon number density at mid-rapidity is reduced by 56%, whereas with the same $\lambda_{\scriptscriptstyle GB}$ but $\beta=0$ this reduction is only 16%. The magnitude of the effects that we see with $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0.025$ further suggests that with these values of the couplings in the bulk gravitational action we are describing collisions in a gauge theory with an intermediate value of the gauge coupling. It is also interesting to note that the product $\tau\rho_{\rm loc}$ at later times is approximately constant, as would be expected if the dynamics becomes boost invariant at late times.

We note of course that baryon number is conserved, meaning that since the total baryon number present in the range of $y$ that we plot in Fig. 8 decreases as $\beta$ is increased, the baryon number density at (much) higher rapidities or on the lightcone increases with increasing $\beta$. These results thus further confirm the increasing transparency of the initial moments of these collisions to baryon number with increasing $\beta$. As we noted previously, the baryon density found near $y=0$ (namely near $z=0$) at late times in our calculations can be reduced to the point that it is close to vanishing in the collision of thin (thick) sheets if we increase $\beta$ to $\approx 0.3$ ($\approx 1$), but these values of $\beta$ are well beyond those for which our linear approximation can be trusted.

IV.4 Hydrodynamization

As promised in Sections II and III, we now check when and to what degree the matter produced in collisions like those in Figs. 1 and 2 hydrodynamizes. In Figs. 9 and 10 we plot selected components of $\mathcal{J}^{\mu}$ and $\mathcal{T}^{\mu\nu}$ at constant-$z$ slices as obtained directly from our holographic calculation together with the values obtained from the hydrodynamic constitutive relations. Recall that in Section III we described the procedure to start from $\mathcal{E}$ and $\rho$ as obtained from the full holographic calculation, obtain the local fluid velocity, and from that determine the local baryon number density in the local rest frame of the fluid, which we denote by $\rho_{\rm loc}$, and then use hydrodynamic constitutive relations to obtain the transverse and longitudinal pressure, denoted $\mathcal{P}_{T}$ and $\mathcal{P}_{L}$, and the baryon number current $J^{z}$. We then compare these three quantities obtained in this fashion via the laws of hydrodynamics to their values as obtained directly from the full holographic calculation. We show this comparison for thick (thin) sheets with $\beta=0.05$ ($\beta=0.025$). We see from Fig. 9 that the transverse and longitudinal pressures, $\mathcal{P}_{T}=T^{x_{\perp}}_{\ x_{\perp}}$ and $\mathcal{P}_{L}=T^{z}_{\ z}$, are well described by first-order hydrodynamics for nearly all times during and after a collision between thick sheets, even before the peak energy density has passed. And, hydrodynamics with a conserved current describes $J^{z}$ well starting around the time at which the energy and baryon number densities peak. From Fig. 10, we see that the transverse and longitudinal pressure are well described by hydrodynamics starting about $0.6/m$ after the time $2.0/m$ at which the energy density produced in a collision of thin sheets peaks, while $J^{z}$ is described well starting a little bit later, about $0.8/m$ after the time when the energy density peaks. Note that the temperature of the plasma at this time is $T_{\rm hyd}\approx 0.32\,m$, where we have used the relationship (59) between $T$ and $\mathcal{E}$, meaning that the hydrodynamization process after the collision of our thin sheets is completed about $0.26/T_{\rm hyd}$ after the energy density peaks, whereas the collision of our thick sheets is close to hydrodynamic throughout.

IV.5 Subleading $\mu/T$ effects for thick sheets

We now turn to the case where $\epsilon\neq 0$, where the lowest order ${\cal O}\left(\frac{\mu^{2}}{T^{2}}\right)$ corrections to $\rho$ and $\mathcal{E}$ are included in our calculation, meaning in particular that we include the backreaction of $\rho$ on $\mathcal{E}$. We are at the same time including what can be thought of as the backreaction of $\rho$ on $\rho$. In bulk terms, $\epsilon\neq 0$ means that we are including the backreaction of the Maxwell field on the bulk geometry, to linear order. In a collision of thick sheets we now find that $\frac{\bar{\mu}}{T}\approx 1.5+{\cal O}(\lambda_{\scriptscriptstyle GB})$ in the region where hydrodynamics is valid, so that $\left(\frac{\mu}{T}\right)^{2}=\epsilon^{2}\left(\frac{\bar{\mu}}{T}\right)^{2}\approx 2\epsilon^{2}+{\cal O}(\lambda_{\scriptscriptstyle GB}^{2})$. In Fig. 11, we show the effect of introducing a nonzero $\epsilon^{2}$, for values up to $\epsilon^{2}=0.75$, on the baryon number distribution after a collision of thick sheets. We see that the effects of choosing $\epsilon^{2}$ as large as 0.75 on the baryon number distribution are very small. If we extend our expression for the fraction of the total baryon number found near one lightcone or the other, $R_{\rm LC}(t)$ as written in (67), to include a term $R_{\epsilon^{2}}(t)\epsilon^{2}$, we find that $R_{\epsilon^{2}}(t=13/m)=0.009$. From these findings, and also from expectations based upon the results of Ref. Casalderrey-Solana et al. (2016), we expect the effect of $\epsilon^{2}$ to be small for thin sheets also.

IV.6 Asymmetric collisions of thin sheets

When colliding two incident sheets that both carry baryon number, it is impossible to determine from which sheet the baryon number deposited in the plasma originally came from. In a strongly coupled theory one could, for instance, imagine a complete ‘bounce’ where all left-moving baryon number on the lightcone after the collision came from the initially right-moving sheet and vice versa. In order to disentangle the contribution to the baryon number in the plasma from reflected and non-reflected baryon number, it is therefore interesting to look at a collision where only one of the incident sheets carries baryon number Casalderrey-Solana et al. (2016). Indeed, we need to do this check in order to confirm our interpretation of our results in terms of baryon transparency, which is to say in order to confirm that when we see increasing baryon number on the left-moving lightcone after the collision as we increase $\beta$, this baryon number does in fact come from the left-moving incident sheet.

In Ref. Casalderrey-Solana et al. (2016), it was found that in a collision of thin sheets — at infinite coupling — at the time $t=7/m$, roughly $40\%$ of the baryon number that arrived on an incident sheet coming from $z<0$ ended up at $z<0$, meaning that it was reflected. At finite coupling we expect the amount of reflected baryon number to be reduced, and indeed that is what we find, as shown in Fig. 12. At the last plotted time of $t=5.7/m$ we find that in the infinite coupling case $38\%$ of the baryon number that came in from $z<0$ has been reflected back, ending up at $z<0$, while $33\%$ of the baryon number is reflected for $\lambda_{\scriptscriptstyle GB}=-0.2$ and $\beta=0.025$. What is perhaps even more interesting is that if we compare the upper-left panel of Fig. 12 to the lower-right panel of Fig. 2, we see that the baryon number that ends up on the positive-going lightcone at late times is almost the same in both cases, meaning that it comes from baryon number going through the collision zone not from reflection. This confirms the interpretation of our results in terms of baryon transparency, and increasing baryon transparency with increasing $\beta$.

V Discussion and Outlook

We have achieved our goal of doing an initial analysis of the dynamics of baryon number, as well as energy-momentum, in a holographic model of heavy ion collisions at a finite value of the gauge coupling, to lowest order in the inverse coupling. In the dual gravitational theory this means that we have worked to linear order in the two new bulk couplings that appear in the most general four-derivative theory with a dynamical metric and gauge field in a 4+1-dimensional asymptotically anti-de Sitter spacetime. We have also assumed a small baryon number chemical potential compared to the equilibrium temperature, $\mu/T\ll 1$, consistent with heavy ion collisions at RHIC and at the LHC.

The two higher-derivative gravity couplings ($\lambda_{\scriptscriptstyle GB}$ and $\beta$) are both thought of as being proportional to the inverse ’t Hooft coupling in our model gauge theory, with proportionality coefficients which could only be uniquely fixed had we derived our bulk theory from string theory. Rather than doing so, we have used a bottom-up approach to developing the model studied here, working with the most general gravitational action to quartic order in derivatives of the bulk metric and gauge fields, given our choice of discrete symmetries and scalings. We have then motivated a qualitative sense for what the reasonable range of values for the bulk four-derivative couplings $\lambda_{\scriptscriptstyle GB}$ and $\beta$ may be. This cannot be done quantitatively as in an effective field theory because we do not know which deformation of ${\cal N}=4$ SYM theory our gravitational action with its inverse-coupling corrections is dual to, but at a qualitative level we have done so by computing the shear viscosity to entropy density ratio $\eta/s$ and the baryon susceptibility normalized by its value in free field theory $\chi/\chi_{0}$ in our model gauge theory and comparing the former to what we know from heavy ion collision phenomenology and the latter to what we know from lattice QCD. Based upon many previous studies Grozdanov and Starinets (2015, 2017); Grozdanov et al. (2016); Grozdanov and van der Schee (2017); Casalderrey-Solana et al. (2018) we expect qualitatively similar behavior to arise from higher-derivative gravitational actions whether they have been derived top-down from string theory or bottom-up, as here.

Our central qualitative result is that as we reduce the gauge coupling from infinity (by increasing $\lambda_{\scriptscriptstyle GB}$ and $\beta$ from zero) the holographic collisions whose dynamics we have calculated become increasingly transparent to baryon number. The baryon transparency turns out to be much more sensitive to $\beta$. As $\beta$ is increased from zero through its reasonable range of values, significantly more of the incident baryon number ends up near the lightcones after the collision — having passed through the collision zone — and less of it ends up stopped in the strongly coupled plasma produced at mid-rapidity. This is in stark contrast with previous results at infinite coupling Casalderrey-Solana et al. (2016). Our present findings show that introducing finite gauge coupling corrections to lowest order changes the dynamics of baryon number transport in these collisions in a direction that makes the final state look more similar to what is seen in heavy ion collisions in QCD. This strongly supports the hypothesis from Section I that this previously observed discrepancy between the dynamics of baryon number in heavy ion collisions and holographic models for collisions can be attributed to the strength of the coupling in the early moments of the collision. Furthermore, we find that for a given nonzero value of $\beta$ the increase in baryon transparency relative to $\beta=0$ is much more significant for holographic collisions of thin sheets than for collisions of thick sheets. This is plausibly consistent with the observation that in QCD as the beam energy is increased and the incident nuclei become more Lorentz contracted the distribution of the baryon number density after the collision becomes more and more different from that in holographic collisions with infinite coupling (with $\mu$ at mid-rapidity dropping and the rapidity at which the baryon number ends up increasing). This observation indicates that finite coupling corrections to the holographic collisions of thin sheets have “farther to go” if they are to take the infinite coupling results and make them more realistic, and indeed we find that the finite coupling corrections are larger for collisions of thin sheets.

We note that we have checked that holographic collisions of incident sheets of energy and baryon number hydrodynamize rapidly, even with the introduction of finite coupling corrections. And, we have checked that for values of $\mu/T$ as large as around 1 the backreaction of the Maxwell field in the bulk action on the gravitational geometry is small, indicating that the back reaction of the baryon number dynamics on the dynamics of the energy density is small. Finally, we have further confirmed the interpretation of our results in terms of increasing baryon transparency with increasing $\beta$ by analyzing the collision of two incident sheets only one of which carries baryon number.

Significant open questions remain. The most obvious next step would be to include corrections that are higher order in $\lambda_{\scriptscriptstyle GB}$ and $\beta$ in the gravitational action in the analysis of holographic collisions of incident sheets carrying baryon number as well as energy density, as well as other higher-than-four-derivative corrections whose contribution to the equilibrium thermal spectrum were first addressed in Ref. Grozdanov and Starinets (2019). For collisions of incident sheets carrying only energy density, without any baryon number, i.e. in the Einstein-Gauss-Bonnet theory with $\beta=0$, the first question was addressed in Ref. Grozdanov and van der Schee (2017) where it was shown that, at $\lambda_{\scriptscriptstyle GB}=-0.2$ (the value we have also used in this work), corrections at ${\cal O}(\lambda_{\scriptscriptstyle GB}^{2})$ are indeed small. The magnitude of sub-leading $\beta$-dependent terms remains to be determined in future work, but because the leading $\beta$-dependent terms that we have computed are so significant in their effects (in particular on the dynamics of baryon number transport in the collisions of thin sheets) we expect to see effects of the sub-leading $\beta$-dependent terms. Working to linear order in $\beta$ as we have done, the increase in the fraction of the total baryon number that ends up near the lightcones is linear in $\beta$. In collisions of thin sheets, with $\beta=0.05$ and $\lambda_{\scriptscriptstyle GB}=-0.2$ that fraction is already more than 50%, and with $\beta=0.1$ and $\lambda_{\scriptscriptstyle GB}=-0.2$ it would be more than 80% if we work to linear order in $\beta$. We therefore expect that in this range of values of $\beta$ the effect of contributions that are higher-order in $\beta$ should become significant.

A more challenging, and more important, future goal is to find a holographic model setting in which one can analyze collisions in which the coupling is weak (or at least finite and intermediate in value, not infinite) at very early times while the coupling at later times, in particular during the hydrodynamic expansion, is strong. Our interpretation of why the baryon number distribution in heavy ion collisions in QCD is so different from that in holographic collisions at infinite coupling is that in QCD the coupling is weak during the first moments of the collision, and if the nuclei are sufficiently highly Lorentz contracted then by the time QGP starts to form the incident baryon number has already passed through the collision zone. Later, the QGP that forms in QCD is indeed a strongly coupled liquid. So, to improve the holographic modelling of heavy ion collisions in QCD, instead of reducing the gauge coupling (by introducing inverse-coupling corrections) at all times as we have done it would be better to explore methods of doing so only at early times. This challenge of course reflects the asymptotic freedom of QCD, but addressing this challenge framed in these terms may be easier than implementing asymptotic freedom in a holographic gauge theory in full. It remains a challenge for future work, however.

Acknowledgements.

We are grateful to Peter Arnold, Jorge Casalderrey-Solana, Christian Ecker, Guy Moore, Andrei Starinets and Urs Wiedemann for helpful conversations. KR and WS gratefully acknowledge the hospitality of the CERN Theory group. We also thank Iain Stewart for making computational resources available. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics through grant DE-SC0011090. WS is supported by the Netherlands Organisation for Scientific Research (NWO) under the grant 680-47-458 (Veni). The work of ÅF was supported by an Aker Scholarship.

Appendix A General four-derivative action and field redefinitions

In this Appendix, we review the construction of the most general four-derivative action involving a metric tensor $g_{\mu\nu}$ and a vector gauge field $A_{\mu}$ in a five-dimensional asymptotically AdS spacetime, and show how to use field redefinitions and a few simplifying choices to cast this action in the form (II) that we use throughout this paper. The presentation largely follows Ref. Myers et al. (2009). The most general action can be written as

[TABLE]

where we split $\mathcal{L}$ into two-derivative and four-derivative terms $\mathcal{L}_{2}$ and $\mathcal{L}_{4}$, respectively. They are

[TABLE]

where $\tilde{L}$ is a length scale that we shall later relate to the AdS radius and

[TABLE]

where we treat all dimensionless coupling constants in $\mathcal{L}_{4}$ as perturbatively small. Other possible terms in the action can be removed by using standard metric identities and integration by parts. Because doing so leaves the dual field theory unchanged, we can now perform first-order field redefinitions $g_{\mu\nu}\to g_{\mu\nu}+\delta g_{\mu\nu}$ and $A_{\mu}\to A_{\mu}+\delta A_{\mu}$, where

[TABLE]

The $\mu_{i}$’s and the $\lambda_{i}$’s are constants that we may choose in order to specify a particular field redefinition. After this redefinition, the gravitational action transforms as $S\to S+\delta S$, where we write the transformations of the two- and four-derivative parts as

[TABLE]

The two-derivative part of the variation is

[TABLE]

while at the four-derivative order,

[TABLE]

with

[TABLE]

The coupling constants $\mathcal{C}_{n}$ in the action obtained after the field redefinition depend on the coupling constants in the original action and on the constants ($\mu_{i}$’s and $\lambda_{i}$’s) that specify the field redefinition. We can now choose these constants, specifying the field redefinition. First, though, since we are only interested in theories without any Chern-Simons-type terms (terms with the Levi-Civita symbol that are odd in $A_{\mu}$), we set

[TABLE]

Next, in order to keep the two-derivative part of the action in ‘canonical’ form after the field redefinition, we choose

[TABLE]

We then choose to fix $\mu_{5}$ and $\lambda_{1}$ as

[TABLE]

and define the AdS radius $L$ in terms of $\tilde{L}$ by

[TABLE]

It is important to notice that neither $\mathcal{C}_{3}=\alpha_{3}$ nor $\mathcal{C}_{6}=\beta_{3}$ can be changed by using the field redefinitions (73)–(74). These constants in the action after field redefinition are (related to) the Gauss-Bonnet coupling $\lambda_{\scriptscriptstyle GB}$ and the new coupling $\beta$:

[TABLE]

In the absence of any terms with the Levi-Civita symbol in the theory, i.e. with (84), we have seven remaining unfixed $\mathcal{C}_{n}$ with $n=\{1,2,4,5,7,8,9\}$ and the remaining freedom to choose five coefficients, $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, $\mu_{4}$ and $\lambda_{2}$. First, we choose $\mu_{1}$, $\mu_{2}$, $\mu_{3}$ and $\mu_{4}$ so as to satisfy the following set of linear equations, because doing so brings the action into the form that gives only two-derivative equations of motion:

[TABLE]

These equations imply that

[TABLE]

We are left with a single free coefficient, $\lambda_{2}$, and three $\mathcal{C}_{n}$’s. We specify $\lambda_{2}$ so as to make $\mathcal{C}_{9}$ to vanish:

[TABLE]

which gives

[TABLE]

There is then no remaining freedom, and the two remaining $\mathcal{C}_{n}$’s are given by

[TABLE]

where we have defined two new couplings $\zeta_{1}$ and $\zeta_{2}$. To first order in four-derivative couplings, this gives the action $S^{\prime}=S+\delta S$:

[TABLE]

Finally, the scaling of four-derivative couplings and the size of the gauge field that we employ in this work (see Section II) are $\lambda_{\scriptscriptstyle GB}\sim\beta\sim\zeta_{1}\sim\zeta_{2}\sim\epsilon^{2}$ and $|A_{\mu}|\sim\epsilon$, with small dimensionless $\epsilon\ll 1$. It is convenient to explicitly rescale $A_{\mu}\to\epsilon A_{\mu}$, so that to order $\epsilon^{4}\sim\lambda_{\scriptscriptstyle GB}^{2}\sim\ldots$,

[TABLE]

Both terms with four powers of $F_{\mu\nu}$ are proportional to $\zeta_{n}\epsilon^{4}\sim\epsilon^{6}$ and can therefore be neglected. Thus, we have arrived at the final form of the action (II) that we use in this paper, which results in purely second-order equations of motion.

Using the above field redefinition, we can connect the thermodynamic relations used in this work with the results discussed in Ref. Myers et al. (2009). The Lagrangian of Ref. Myers et al. (2009), up to a rescaling of $A_{\mu}\rightarrow A_{\mu}\tilde{L}$, is obtained from (71) and (A) by setting

[TABLE]

Since $c_{5}$ and $\kappa$ do not affect thermodynamics, we can set them to zero. Working with the assumption of a small chemical potential, i.e. that $A_{\mu}\sim O(c_{i})$, the results of Ref. Myers et al. (2009) for the baryon number and conductivity to the order we are working at are

[TABLE]

where we have adjusted for a relative factor of $\pi$ in our definition of the chemical potential. We use tildes to indicate the temperature and the chemical potential for the theory in Ref. Myers et al. (2009). Now, it turns out that at the boundary, the field redefinitions (73) and (74) with the choice of the coefficients in (107) simply reduce to global rescalings

[TABLE]

with

[TABLE]

In terms of the boundary metric, this is a constant Weyl transformation, under which the transformation properties of thermodynamic and hydrodynamic quantities are well known (see Refs. Baier et al. (2008); Loganayagam (2008); Grozdanov and Starinets (2015)). They follow from the condition that the generating functional of the conformal field theory must remain invariant under such transformations, giving

[TABLE]

An alternative way to state this, from the point of view of gauge-gravity duality, is that the generating functional (or the partition function) must remain invariant under field redefinitions. Following analogous logic, we impose that the generating functional must also remain invariant under the rescaling transformation of the gauge field, which gives

[TABLE]

The transformation of $\mu$ follows directly from the identification of the chemical potential with the gauge field on the boundary. The transformation of the local baryon number density can be found from the fact that $\mu\rho$ should transform the same way as $\mathcal{E}$, which does not change under the rescaling of the gauge field. The conductivity transformation follows from the Kubo formula, which schematically relates $\sigma\sim[\rho,\rho]$. Thus, we find the following relation between thermodynamic quantities and transport coefficients in our theory and the theory of Ref. Myers et al. (2009):

[TABLE]

Furthermore, we solve (90) to give

[TABLE]

With the above relations we can now use (108) to solve for $\sigma$ and $\rho$ in terms of $L$, $T$, and $\mu$, finding

[TABLE]

Plugging in $L=1+\lambda_{\scriptscriptstyle GB}/2$, $\mu=\epsilon\bar{\mu}$ and setting $\kappa_{5}^{2}=2$ to conform with the conventions used in this work (cf. Eq. (35)), we arrive at equations (60) and (64). Continuing further, we can also show agreement between the entropy, free energy and energy density from Ref. Myers et al. (2009) after the field redefinitions that we have specified and the same quantities in our theory. We note that the agreement between our entropy, calculated using the formula $A/4G_{N}$, and the entropy obtained from the result of Ref. Myers et al. (2009) after field redefinitions, which is calculated from the (higher-derivative) Wald formula Wald (1993); Iyer and Wald (1994), confirms that our theory (in the case of a black brane) has the property that the higher derivative corrections to $A/4G_{N}$ appearing in the Wald formula vanish, meaning that the entropy is given by the Bekenstein-Hawking result Wald (1993); Iyer and Wald (1994).