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

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.

Previous authors [1] 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 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 λ ≡ g 2 N c , with g the gauge coupling. Although 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 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. [2][3][4][5].) The nonzero density in the 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. [1] can be thought of as a toy model for heavy ion collisions with nonzero baryon number density. The authors of [1] 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 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. [10]. These authors have performed a holographic study of coupling-constantdependence in holographic collisions in a hypothetical quantum field theory at large N c and 't Hooft coupling λ, including inverse-λ corrections to the large λ limit, albeit in collisions with no baryon number. The gauge theory here is typically assumed to be some deformation of 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. [1] to include finite coupling corrections. Equivalently, our goal is to extend the calculation of Ref. [10] 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 ρ. The baryon number chemical potential µ is thermodynamically conjugate to ρ. In this work, we begin from a model for the formation and evolution of such a state in a non-Abelian gauge theory (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 energymomentum and baryon number density of a four-dimensional state with an infinite 't Hooft coupling λ is described by an asymptotically Anti-de Sitter (AdS) geometry in five-dimensional Einstein-Maxwell theory. Its bulk action is where Λ ≡ −6/L 2 is the negative cosmological constant with L the AdS radius, and where the five-dimensional Newton's constant κ 5 is inversely proportional to the number of colors in the fourdimensional gauge theory: κ 5 ∝ 1/N c 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 µν is sourced by the five-dimensional metric perturbation and the four-dimensional conserved global current J µ in the boundary theory is sourced by the gauge field in the bulk, according to the standard holographic dictionary. (See, e.g., Ref. [11].) We will think of the dynamics of J µ as a model for the dynamics of baryon number current.
The holographic "experiment" of a collision of two sheets of energy -and baryon numberdensity 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 [12][13][14].) 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. [12][13][14][15][16][17][18][19][20][21][22] and the reviews in Refs. [5,11,23,24]. Collisions of sheets of energy density with nonzero µ and consequently with nonzero baryon density ρ were studied via their dual description in terms of the five-dimensional Einstein-Maxwell theory (2.1) in Ref. [1].
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 λ. 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 λ 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 (topdown) example is the leading-order 't Hooft coupling correction to the dynamics of energy and momentum in SU (N c ), N = 4 SYM theory. At zero baryon number density, the leading part of the corrected dimensionally reduced gravity action is where W is a contraction of Weyl tensors to the fourth power and ζ(z) is the Riemann zeta function. For details, see Refs. [25][26][27]. In type IIB string theory, the tension of the fundamental string sets the dual 't Hooft coupling of N = 4 SYM theory, i.e. α ∼ 1/ √ λ. Thus, a perturbative supergravity expansion in α , which starts at α 3 gives rise to field theory observables computed in a perturbative expansion to O(1/λ 3/2 ).
In a generic string theory, the corrections to the Einstein-Hilbert action start at O(α ) with curvature-squared terms ∼ R 2 instead of W ∼ R 4 . The most general such action is 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 α i would represent some combinations of corrections that are ∼ α and corrections that are ∼ 2 P /L 2 , where P is the Planck scale, with κ 2 5 = 3 P . While α corrections correspond to inverse coupling corrections in the field theory, 2 P /L 2 corrections correspond to 1/N c corrections. (See Refs. [10,28].) 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 higherderivative corrections as if they induce coupling constant corrections in the dual field theory. By performing a field redefinition of the metric g µν to O(α i ) (see e.g. Refs. [29,30]), the action (2.3) can be brought into the form where λ GB = 2α 3 . This is the Einstein-Gauss-Bonnet action, which gives rise to second-order equations of motion. In this theory, the coupling λ GB can in principle be treated non-perturbatively as long as it lies in the interval λ GB ∈ (−∞, 1/4] [29,30]. In this work, we will treat λ GB ∼ α 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 O(λ 2 GB ) in Ref. [10]. Now, in order to study coupling-dependent collisions at finite µ, we add all possible fourderivative terms involving the gauge field and its coupling to the metric with the exception of Chern-Simons-type terms odd in A µ , 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 λ GB . To find the relevant terms in the bulk action, we start by writing the most general four-derivative action involving g µν and A µ [28,[31][32][33]. Then, we perform the field redefinitions of both g µν and A µ 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: λ GB from Eq. (2.4) and a single new constant β that characterizes the strength of the higher-derivative coupling between g µν and A µ . A particularly convenient choice for the action that will be used in all of our analyses is [28] which, in analogy with the Einstein-Gauss-Bonnet theory in Eq. (2.4), gives rise to purely secondorder equations of motion for all components of g µν and A µ . Since both multiply four-derivative terms in the action, we shall assume throughout that β is of the same order as λ GB . We furthermore assume that µ 2 /T 2 is perturbatively small and of the order of the λ GB and β; we shall discuss this further below. We have introduced a control parameter as a device for later convenience; it could evidently be scaled away by redefining A µ .
Since µ/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 ρ that is much smaller than E 3/4 , with E their energy density. Consequently, µ/T is small in the hydrodynamic fluid produced after the collision. This assumption allows additional simplifications, as we explain below. Small ρ in the boundary theory means small A µ in the dual Einstein-Maxwell theory, and in particular means that the backreaction of the Maxwell field A µ on the five-dimensional geometry is small. We can implement this assumption by treating 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 µ ∝ and µ ∝ . To see this, let us for a moment consider the action with the normalization of the Maxwell field used in Eq. (2.1). Then ρ and µ are set by the size of the bulk gauge field |A µ | which has a mass dimension of 1/L. If we assume that ρ and µ are small and parametrically O( ), we have that |A µ | ∼ /L. It is then convenient to explicitly scale A µ → A µ and take A µ of order ∼ 1/L, which introduces powers of into the action as in Eq. (2.5). Note that the assumption of the smallness of A µ in Eq. (2.1) and consequently the smallness of in Eq. (2.5) is an assumption about the smallness of ρ in the incident sheets of energy and baryon number density. Thus, in Eq. (2.5) 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. (2.5) are λ GB and β.
Introducing as we do gives us a convenient way to take the smallness of ρ into account in a combined perturbative expansion: we shall assume that 2 is of order λ GB so that we can arrange our calculations in a single perturbation series in λ GB . We will then determine T µν and J µ / to 1 Note 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 [1,34,35]. 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. leading order, O(λ GB ). Our motivation here is simply that we wish to take λ GB , β and (and consequently ρ) to all be small and work to lowest nontrivial order in a combined expansion: the particular choice of scaling 2 ∼ |λ GB | is arbitrary; others could be investigated. However, if we had instead assumed ∼ |λ 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 O(λ GB ), and we would effectively be working at zero chemical potential rather than at small chemical potential. Since in our scheme 2 , λ GB and β scale in the same way we will frequently express 2 and β in terms of λ GB as where our choice of scaling means that both c and c β are O(1) when we count powers of .
With the scaling that we have described, the action Eq. (2.5) yields the most general equations of motion to linear order in λ GB . The Maxwell (F 2 ) and Gauss-Bonnet (R 2 ) terms are of the same order (since 2 ∼ λ GB ) while the term proportional to β 2 (RF 2 ) is subleading. However, the factor 2 cancels in the Maxwell equations, so the β term provides the leading correction to the Maxwell equations. Further details regarding the construction of the action Eq. (2.5) and the field redefinitions involved are discussed in Appendix A. We close our discussion here by noting that if we had not assumed that (and consequently ρ and µ/T ) is small then we would have had to include additional F 4 terms in the action (2.5).
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 energymomentum tensor T µν and a current J µ , respectively. Both conserved tensors need to be calculated to first order in our combined expansion in small λ GB , β and 2 . This means that T µν and J µ will take the form , (2.7) where we note that µ 2 /T 2 ∝ 2 and where all of T µν , J µ , and all of the δT µν and δJ µ 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 µν and J µ 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 where u µ is the velocity field normalized to u 2 = −1, µ and T are treated as hydrodynamic nearequilibrium quantities, and ∆ µν = u µ u ν + η µν is the projector in directions transverse to u µ . The thermodynamic quantities appearing in these expressions are the energy density E, pressure p and baryon number density ρ. Assuming local equilibrium, they satisfy the thermodynamic relation with s the entropy density, and all of them vary in space and time on wavelengths that are long compared to 1/T . Because N = 4 SYM theory and its deformations considered in this work are conformal, the equation of state is given by E = 3p. The two transport coefficients that enter into the hydrodynamic constitutive relations (2.7) and (2.8) are the shear viscosity η and the baryon number conductivity σ. (Recall that J µ is not the electric current, but the baryon current.) However, if J µ were to be weakly gauged, with the introduction of an additional U (1) gauge field, then σ would indeed become the electrical conductivity of the strongly coupled field theory. For further details on relativistic hydrodynamics, see Ref. [36]. We shall compute E and p to linear order in λ GB , meaning to O( 2 ), and in so doing will precisely reproduce the analysis of finite-coupling corrections in the absence of any baryon number density from Ref. [10]. Note that there are no O( 3 ) terms (∼ µ 3 or ∼ λ GB µ) in E and p, since these thermodynamic quantities are even in µ. We shall compute ρ, which at infinite coupling is proportional to for small ∼ µ, to order λ GB ∼ 3 . In terms of the gravity couplings λ GB and β, only λ GB will influence the energy density E and pressure p, while both λ GB and β will affect the baryon number density ρ. Working to O( λ GB ) in the analysis of ρ 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 O(1) quantitiesρ ≡ ρ/ ,μ ≡ µ/ andJ µ ≡ J µ / . Note that all these quantities are zeroth order in , meaning in particular that they are finite in the → 0 limit.
Before we turn to a detailed analysis of the action (2.5) 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 λ GB and β in our gravitational action. This certainly cannot be done in a quantitative way because although we do not know exactly what deformation of 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. (2.3), if we had derived the action (2.5) from a top-down stringy construction we would know precisely which gauge theory it is dual to and in such a construction both λ GB and β would have been fixed in terms of α and 2 P /L 2 . Instead, in our bottom-up approach reasonable values for these quantities can be estimated by asking that λ GB -and β-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 η/s and the baryon number susceptibility χ in QCD and ask how this compares with their dependence on λ GB and β 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 η/s and χ 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 ρ in collisions.
We start from a relatively well-established fact that when λ GB R 2 terms in the gravitational dual are treated as enabling finite coupling corrections in the gauge theory, it is natural to take λ GB to have a negative sign. (See e.g. Refs. [10,27,[37][38][39].) One way to motivate this sign is to look at the result for η/s in (zero baryon number) Einstein-Gauss-Bonnet (2.4), where it is known nonperturbatively in λ GB and is given by [29,40] meaning that if we choose λ GB negative η/s increases with increasing |λ 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 η/s on the gauge coupling g at small g is given by η/s ∼ const/[g 4 ln(1/g)] [41], meaning that at weak coupling η/s increases with decreasing gauge coupling and diverges as g → 0. We will typically choose λ GB = −0.2, which corresponds to an 80% increase in the value of η/s relative to its value when λ GB = 0 and the gauge theory has infinite coupling. There is nothing sacred about this particular value, in particular given that η/s in QCD will in reality have some temperature dependence, but this choice puts η/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 [42][43][44]; for reviews, see Refs. [5,45,46].
At small but nonzero baryon number density, perturbative results in both the four-derivative couplings found in the Maxwell-Einstein-Gauss-Bonnet action (2.5) give rise to corrections to the ratio of η/s that take the form the ellipsis denotes terms of ∼ O(λ 3 GB ) and beyond. Sufficiently little is known about the µdependence of η/s in QCD that we will not attempt to use a calculation of the µ-dependent contributions to (2.13) to estimate what range of β we might investigate. Another reason why attempting this would be perilous is that the presently unknown contributions to η/s that are of order λ 2 GB are parametrically just as significant, given the scaling that we are using, as the known 2 β and 2 λ GB terms, both of which can be computed from the field redefinitions discussed in Appendix A and the calculations of Refs. [31,32].
We shall estimate a reasonable range for β by looking at its effects on the baryon number susceptibility χ. Note that in our model theory, χ ∝ N 2 c because all degrees of freedom are in color-adjoint representations whereas in QCD, χ ∝ 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 χ in the theory with finite coupling to χ 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. (3.35)) that to linear order in µ and to lowest order in λ GB and β the baryon number density is given (for a choice of normalization, see Section III) by giving a susceptibility From Refs. [11,47,48] we know that in N = 4 SYM theory we have χ ∞ /χ 0 = 1 2 where χ ∞ is the susceptibility in the infinite coupling limit (corresponding to λ GB = β = 0 in the above), while χ 0 is the susceptibility in the non-interacting limit. This gives χ 0 = π 2 T 2 in our normalization scheme. Now we can calculate the ratio χ/χ 0 , namely the ratio between χ at the intermediate coupling that corresponds to given values of λ GB and β to χ in the noninteracting limit, obtaining This can be compared to lattice calculations of the same ratio in QCD [49], which indicate that (2.17) That said, we caution that working only to linear order in β 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 (2.16) that it starts to break down once 16β gets as large as 1. Another way of making a reasonable guess is to note from (2.16) that values of λ GB = −0.2 and β = 0.1 imply that ρ increases by a factor 2, suggesting that the linear approximation is breaking down. And, if we express µ as a function of ρ and T we see that to the order we are working it is proportional to 1 − 3λ GB − 16β, 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 β in the range (2.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 β may not be sufficient. Summarizing these considerations, it seems reasonable to investigate the consequences of choosing λ GB small and negative and β small and positive in the action (2.5) on the dynamics of ρ 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 (2.5).

III. ANALYSIS OF THE MODEL AND DETAILS OF THE NUMERICAL SETUP
In this Section we proceed with the analysis of the action (2.5). 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 (2.5) 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 λ GB and assuming the scaling λ GB ∼ β ∼ 2 ∼ µ 2 /T 2 . The equations of motion that follow from variations of g µν and A µ were derived in Appendix E of Ref. [28]. The resulting Einstein equations are where the perturbatively small terms on the right-hand-side of the equation are The Maxwell equations, as modified by the four-derivative terms, are Both the equations of motion (3.1) and (3.5) contain at most second derivatives with respect to each of the coordinates. The fact that, with our scaling, T β µν ∼ O( 2 β) ∼ O(λ 2 GB ) implies that to leading order in λ GB the β-term in the action only affects the Maxwell equations (3.5). We shall employ the convention 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 λ GB corrections.
Equations (3.1) and (3.5) have a solution of the following form: where x ± = t ± 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 + ↔ x − . The wave solution is planar (and translation-invariant) in the two-dimensional x ⊥ 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 (3.7)-(3.8) and their counterparts with x + ↔ x − , as initial conditions. Specifying our initial conditions therefore requires choosing the eight functions h i (x ± ) and a ±(i) (x ± ). We choose the four profile functions h 0 (x ± ) and a +(0) (x ± ) all to be Gaussians with the same width w: 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 ± ) and a ±(1) (x ± ) needed to complete the specification of our initial conditions such that the total O(λ 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 β and λ GB . This will allow us to compare simulations with different β and λ 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 [15] in which the metric takes the form 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 x ⊥ are perpendicular to z in spatial R 3 . The metric (3.10) 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 = ∞. For the gauge field, we choose the radial gauge, meaning that its form reads In the numerical evolution, we work directly with the field strength F µν . The coordinates used in Eq. (3.10) 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 [15,50]. 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. [19,22,51,52]). These include in particular the (tricky) coordinate transformation from (3.7) to (3.10), which luckily is straightforward to generalize to our λ GB -corrected set-up (see also [10]).
For our specific case of Gauss-Bonnet gravity with a Maxwell field, several aspects are noteworthy. To accommodate a perturbative expansion in λ GB (with our assumed scaling) we expand the metric functions, the gauge field, and later also various thermodynamic quantities as power series in λ GB , namely f = n f (n) λ n GB (cf. Eq. (2.6)). For the metric this implies and similarly for all other functions. At every time step in the numerical evolution (solving Eqs. , ∂ t F (0)rt and ∂ 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) , . . .}. This means that they affect the perturbed firstorder solutions, which means that we need to keep track of ∂ 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 + O(λ GB ). (Note, however, that these corrections were computed in Ref. [10] in order to see the effects of λ 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 (3.1) and (3.5), which take the form 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 (3.13), (3.14), (3.15) and (3.16). To make this extraction easier, when we solve for the time evolution we recast the equations for A µ , C, B and F as equations for r 2 A µ , r 2 (C − r 2 ), r 4 B and r 2 F , and solve them numerically in that form [51].
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 µν and the conserved U (1) (baryon number) current J µ . We will work with the rescaled quantities The normalization factor of κ 2 5 /2L 3 0 is given by 2π 2 /N 2 c in N = 4 SYM theory at infinite coupling, i.e. when λ GB = β = 0. L 0 is the λ GB = 0 value of the AdS radius, which for convenience we have set equal to 1, see (3.6) [10]. A particularly straightforward way of ensuring that these normalization conventions are satisfied is to set κ 2 5 /8π = 1/4π 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. (3.14)-(3.16) that are determined by the time evolution in the bulk gravitational theory can be obtained to O(λ GB ) in the same way as they were previously obtained in Einstein-Gauss-Bonnet theory without a bulk gauge field, and take the form [28,53] with 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. (3.13). 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 is the canonical momentum conjugate to the bulk gauge field A ν , satisfying Π µν = −Π νµ , and with L being the Lagrangian corresponding to the action (2.5). The covariant divergence of an antisymmetric tensor satisfies ∇ µ Π νµ = 1 where, to avoid confusion, in this paragraph alone we have found it necessary to use Latin indices a, b, . . . for the indices that run over the four-dimensional boundary coordinates only. From (3.23) we see that √ −g Π ra is a conserved current within each constant-r slice. We can therefore define the conserved boundary current as Inserting the near-boundary expansion of the gauge field (3.13) and the metric functions (3.14)-(3.16) into this expression, we find which is the relationship that we need. We find no divergences in the r → ∞ limit of √ −g Π ra , so no counterterms are needed. Note also that we have computed Π µν by taking the derivative with respect to A µ , so that J ∼ O( ), which is equivalent to first calculating Π µν and then recaling A µ → A µ , as described in the previous Section. We compute the baryon number density from ρ = J 0 . Finally, in (3.25) and (3.26) we can safely replace the Latin index a by a Greek index µ, 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 T µν and J µ . 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 T µν and J µ and test the validity of the hydrodynamic approximation (2.9)-(2.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 (3.13)-(3.16) and appearing in the energy-momentum tensor (3.19) and the current (3.25).
In order to make a comparison to hydrodynamics, we need to obtain the local fluid velocity u µ and then the local energy density E loc and the local baryon number density ρ 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 µ as the timelike eigenvector of −T µ ν (x), hence, The local baryon number density is then To evaluate (2.9) and (2.10) we also need the equilibrium expressions for T and µ as functions of E and ρ. (We shall use these expressions to obtain the local T and µ from E loc and ρ 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 λ GB , the black brane solution is given by where g 0 , g 1 , k 0 , k 1 , Q 0 and Q 1 are free parameters. The Hawking temperature T , the chemical potential µ, and the entropy density s of the black brane solution are given by where the horizon radius r + ≡ r 0 +λ GB r 1 +O(λ 2 GB ) is defined as the largest solution to C(r) = 0 and where Σ 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 [54], and, as we discuss in Appendix A we find that the same is true when we add the Maxwell and Maxwellcurvature coupling terms. Recall, however, that what we need is expressions that relate T and µ to E and ρ. We can use (3.20), (3.21) and (3.26) to express the energy density E, pressure p and baryon number density ρ as functions of k i , Q i and r i . Then, we can solve equations (3.31)-(3.33) together with C(r + ) = 0 for k i , Q i and r i in terms of T andμ ≡ µ/ . Plugging this into our expressions for E, p and ρ we find where we remind thatρ ≡ ρ/ . Finally, we can solve these expressions forμ and T which yields (3.37) The expressions (3.36) and (3.37) are the expressions that we shall use in order to obtain the local T and µ from E loc and ρ loc , so that we can then use (3.27) and (3.28) in the evaluation of the expressions (2.9) and (2.10) for the two conserved tensors in the hydrodynamic description of the fluid.
In order to complete the evaluation of (2.9) and (2.10) we also need the transport coefficients η/s and σ/T . We discussed the first of these already in the previous Section; to the order at which we are working, it is given by The conductivity can be obtained by using the relation σ = Dχ, where χ is the baryon susceptibility and D the baryon diffusion constant. D was calculated nonperturbatively in both couplings at µ = 0 in our theory in Ref. [28] and perturbatively at µ = 0 in a theory related to ours by field redefinitions in Ref. [31]. Appendix A contains a discussion on how to convert the result of Ref. [31] to our theory. The result, with the conventions given in Eqs. (3.6) and (3.18), is (3.39) We shall take E loc and ρ loc from our holographic simulation, use (3.36) and (3.37) to obtain T and µ, and then use (3.27) and (3.28) to evaluate (2.9) and (2.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 is indeed satisfied. Note that the β-dependent correction to J µ in Eq. (3.26) contributes to (3.40) only at O(λ 2 GB ), since ρ and µ are O( ), so that the leading β-independent contribution of µρ is already O(λ GB ). We have computed a black brane solution and an expression for T µν also to O(λ 2 GB ) and verified that (3.40) 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.

A. Colliding thick and thin sheets
We begin the discussion by considering the case with = 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 O µ 2 T 2 corrections toρ and E are neglected. On the bulk side of the duality, setting = 0 means that we neglect the backreaction of the gauge field dynamics on the Einstein equations (3.1). The evolution of energy and momentum in the = 0 system is therefore completely equivalent to what was studied in Ref. [10] in Einstein-Gauss-Bonnet theory. Note that cancels from the Maxwell equation (3.5), which implies that the dynamics of the baryon number density whose time evolution we shall follow remains (highly) nontrivial even at = 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 = 0 and λ GB = −0.2 for thick and thin sheets, respectively. To the order that we are working, the energy density has no dependence on β. We plot the baryon number density for two values of β in each case, which we choose to be β ∈ {0, 0.1} for thick sheets and β ∈ {0, 0.025} for thin sheets. We choose the (smaller) values of β that we use for thin sheets in order to be safely within the range of values of β where all the qualitative considerations that we discussed at the end of Section II indicate that working to linear order in β is likely a safe approximation. And indeed, if we increase β further to β = 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 β is inadequate. That said, in the case of collisions of thick sheets we find that β's in this range have only modest effects and so for this case we have been more bold and extended our calculations up to β = 0.1. It will take future higher-order calculations to reliably estimate up to what β 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 = 0 our results for the energy density are the same as in Ref. [10], 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. In this work, we are primarily concerned with the analysis of the effects of β and λ 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 ≥ 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 λ GB = −0.2 and β nonzero (0.1 for the thick sheets; 0.025 for the thin sheets) as well as a thin-dashed curve for the theory with λ GB = −0.2 and β = 0 (whose gravitational dual is Einstein-Maxwell-Gauss-Bonnet theory) and a thick-dashed curve for the theory with λ GB = β = 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 λ GB = −0.2 with β = 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 [10]. 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 λ GB when β = 0, the thick sheets come to a full stop after which point the plasma expands.
When we turn on a nonzero β (specifically β = 0.1), we see in Figs 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 β at a fixed value of λ 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 λ GB = −0.2 and β = 0 is modest, while for the case of λ GB = −0.2 and β = 0.025 (not β = 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 β-dependence of this increased transparency by focusing on the β-dependence of the baryon number density distribution in collisions of thin sheets at t = 5/m, the latest time plotted in Fig. 4. Increasing β to β = 0.05 further increases the effect. We find that if we were to set β ≈ 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 β ≈ 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 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: We find that at t = 13/m after the collision of thick sheets in Figs In both cases, increasing β 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 β increases with decreasing sheet width: the increase in transparency becomes more sensitive to β 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 λ GB and β. This means that the baryon number density ρ, at fixed (z, t), is a linear function of β and λ GB (and of 2 , which has so far been set to zero). It follows that R LC (t) is also linear in the same variables, and we can define the coefficients R i as In Fig. 7 we show the time evolution for the coefficients R 0 , R β and R λ GB . Indeed, we observe that the β term is responsible for the dominant contribution to the late time baryon number on the lightcone and that the β 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 Comparing Eqs. (4.3) and (4.4) demonstrates that the effects of λ GB and β 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 λ GB = −0.2 and compare R LC with β = 0.05 to R LC with β = 0, for thick sheets R 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 (2.5), it is natural to ascribe this behavior to the fact that the two bulk couplings λ GB and β 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 β on the dynamics is much more strongly dependent on the width of the sheets than that of λ GB . The effect of β is illustrated in Figs. 5 and 6. Since E and ρ are linear in λ GB , β and 2 , the trends seen upon varying β in these Figures are also obtained for any other values of λ GB and 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 β corresponds to studying holographic collisions with a reduced gauge coupling. It is therefore tempting to speculate that the "reason" that a given value of β 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 infinitecoupling 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 β that we have employed are around what we described as the reasonable range, (2.17), based upon consideration of the effect of β on the baryon susceptibility in equilibrium. But, at the same time, these values of β are in a regime where, as we expected based upon the discussion at the end of Section II, we find β-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 β 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 β (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 λ GB ) are likely important. Cranking up β increases the baryon number on the lightcones and if we work only to linear order this increase is linear in β. This increase must be tamed at large enough β by corrections that are higher order in β. We leave the calculation of higher-order-in-inverse-coupling corrections to future work.

C. Spacetime rapidity distribution
Next, to further illustrate our results we transform from Minkowski coordinates (t, z) to proper time τ = √ t 2 − z 2 and spacetime rapidity y = arctanh(z/t), compute the baryon number density in the local fluid rest frame ρ loc , and plot the spacetime rapidity distribution of ρ 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 λ GB on the rapidity profile of the energy density analyzed in Ref. [10]. For the baryon number, however, there are in addition effects introduced by the β-dependent terms. And, as we saw before, the consequences of reasonable nonzero values of β can be large. In particular, at λ GB = −0.2 and β = 0.025, at a proper time τ = 2.5/m after the collision of thin sheets the baryon number density at mid-rapidity is reduced by 56%, whereas with the same λ GB but β = 0 this reduction is only 16%. The magnitude of the effects that we see with λ GB = −0.2 and β = 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 τ ρ loc at later times is approximately . The presence of the λ GB corrections makes the rapidity distribution of the baryon number wider and lower, with more baryon number at higher rapidity. The effect of β is similar, but much stronger. Note also that at mid-rapidity the product τ ρ loc is approximately constant at late times, as would be expected if the dynamics approaches boost-invariance at late times. The wiggles in the results from our β = 0.025 calculations plotted in the right panel are artifacts of our numerical discretization; as we go to finer discretizations, the wiggles decrease while in other respects the results do not change.
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 β is increased, the baryon number density at (much) higher rapidities or on the lightcone increases with increasing β. These results thus further confirm the increasing transparency of the initial moments of these collisions to baryon number with increasing β. 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 β to ≈ 0.3 (≈ 1), but these values of β are well beyond those for which our linear approximation can be trusted.

D. 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 J µ and T µν 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 E and ρ 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 ρ loc , and then use hydrodynamic constitutive relations to obtain the transverse and longitudinal pressure, denoted P T and P L , and the baryon number current J z . We then compare these three quantities obtained in this fashion via the laws of The calculation is done with = 0, λ GB = −0.2 and β = 0.05. The calculation is done at z = 5/m, meaning that the lightcone is at t = 5/m. We also plot E and ρ, which we obtain from the full simulation and use as inputs to the hydrodynamic constitutive relations; the outputs of these relations are then the hydrodynamic "predictions" for P L , P T , and J z . We also compare the hydrodynamic prediction for ρ obtained from ρ loc via the hydrodynamic constitutive relation to that which we take from the full simulation. The hydrodynamic predictions for ρ and J z are related through u µ in that u µ J µ = u µ J µ hyd = −ρ loc . hydrodynamics to their values as obtained directly from the full holographic calculation. We show this comparison for thick (thin) sheets with β = 0.05 (β = 0.025). We see from Fig. 9 that the transverse and longitudinal pressures, P T = T x ⊥ x ⊥ and 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 corrections to ρ and E are included in our calculation, meaning in particular that we include the backreaction of ρ on E. We are at the same time including what can be thought of as the backreaction of ρ on ρ. In bulk terms, = 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μ T ≈ 1.5 + O(λ GB ) in the region where hydrodynamics is valid, so that µ Fig. 11, we show the effect of introducing a nonzero 2 , for values up to 2 = 0.75, on the baryon number distribution after a collision of thick sheets. We see that the effects of choosing 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 LC (t) as written in (4.2), to include a term R 2 (t) 2 , we find that R 2 (t = 13/m) = 0.009. From these findings, and also from expectations based upon the results of Ref. [1], we expect the effect of 2 to be small for thin sheets also.

F. 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 [1]. 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 β, this baryon number does in fact come from the left-moving incident sheet.
In Ref. [1], 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 λ GB = −0.2 and β = 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 β.

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, µ/T 1, consistent with heavy ion collisions at RHIC and at the LHC.
The two higher-derivative gravity couplings (λ GB and β) 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  Fig. 2 when both incident sheets carried baryon number, the peak in the baryon number density near z = t = 0 is much smaller here. The absence of the large peak near z = t = 0 is not surprising, since when both sheets carry baryon number this comes directly from summing the two original sheets, but the presence of the small peak here as well as the charge density on the negative-moving lightcone indicates that the sheet of energy density that comes in at the speed of light carrying no baryon number immediately obtains some baryon number from the other sheet. Initially, at mt = 1 this effect is stronger at infinite coupling (see lower panel), but at later times there is more charge on the lightcone for the collision done with nonzero λ GB and β, due to the weaker coupling. The amount of baryon number that ends up going through the collision zone and staying near the positive-going lightcone at late times is however almost the same as in the case of two colliding charged sheets. In the lower panel, we look at the time evolution of the baryon number density distribution at five time slices. In the upper-right panel, we compute the fraction of the baryon number that has been reflected (which is to say that is at z < 0 after t = 0) as a function of time for λ GB = −0.2, β = 0.025 compared to the infinite coupling case, λ GB = β = 0. sense for what the reasonable range of values for the bulk four-derivative couplings λ GB and β may be. This cannot be done quantitatively as in an effective field theory because we do not know which deformation of 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 η/s and the baryon susceptibility normalized by its value in free field theory χ/χ 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 [10,27,28,30,55] 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 λ GB and β 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 β. As β 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 [1]. 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 β the increase in baryon transparency relative to β = 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 µ 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 µ/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 β 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 λ GB and β 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. [56]. For collisions of incident sheets carrying only energy density, without any baryon number, i.e. in the Einstein-Gauss-Bonnet theory with β = 0, the first question was addressed in Ref. [10] where it was shown that, at λ GB = −0.2 (the value we have also used in this work), corrections at O(λ 2 GB ) are indeed small. The magnitude of sub-leading β-dependent terms remains to be determined in future work, but because the leading β-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 β-dependent terms. Working to linear order in β as we have done, the increase in the fraction of the total baryon number that ends up near the lightcones is linear in β. In collisions of thin sheets, with β = 0.05 and λ GB = −0.2 that fraction is already more than 50%, and with β = 0.1 and λ GB = −0.2 it would be more than 80% if we work to linear order in β. We therefore expect that in this range of values of β the effect of contributions that are higher-order in β 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. and show how to use field redefinitions and a few simplifying choices to cast this action in the form (2.5) that we use throughout this paper. The presentation largely follows Ref. [31]. The most general action can be written as where we split L into two-derivative and four-derivative terms L 2 and L 4 , respectively. They are whereL is a length scale that we shall later relate to the AdS radius and where we treat all dimensionless coupling constants in 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 µν → g µν + δg µν and A µ → A µ + δA µ , where δg µν ≡ µ 1L 2 R µν + µ 2L 4 F µλ F λ ν + µ 3L 2 Rg µν + µ 4L 4 F ρσ F ρσ g µν + µ 5 g µν , (A4) The µ i 's and the λ 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 → S + δS, where we write the transformations of the two-and four-derivative parts as The two-derivative part of the variation is L 2 + δL 2 = 1 + 6µ 1 + 30µ 3 + 3µ 5 2 R + 12 L 2 1 + 5µ 5 2 + while at the four-derivative order, (L 4 + δL 4 ) /L 2 = C 1 R 2 + C 2 R µν R µν + C 3 R µνρσ R µνρσ + C 4L 2 RF µν F µν + C 5L 2 R µν F µρ F ρ ν + C 6L 2 R µνρσ F µν F ρσ with The coupling constants C n in the action obtained after the field redefinition depend on the coupling constants in the original action and on the constants (µ i 's and λ 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 µ ), we set κ 1 = κ 2 = γ 1 = λ 3 = 0. (A15) Next, in order to keep the two-derivative part of the action in 'canonical' form after the field redefinition, we choose 12 We then choose to fix µ 5 and λ 1 as λ 1 = 6 (µ 1 + 2µ 2 + 5µ 3 + 10µ 4 ) , and define the AdS radius L in terms ofL by It is important to notice that neither C 3 = α 3 nor C 6 = β 3 can be changed by using the field redefinitions (A4)-(A5). These constants in the action after field redefinition are (related to) the Gauss-Bonnet coupling λ GB and the new coupling β: In the absence of any terms with the Levi-Civita symbol in the theory, i.e. with (A15), we have seven remaining unfixed C n with n = {1, 2, 4, 5, 7, 8, 9} and the remaining freedom to choose five coefficients, µ 1 , µ 2 , µ 3 , µ 4 and λ 2 . First, we choose µ 1 , µ 2 , µ 3 and µ 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: These equations imply that We are left with a single free coefficient, λ 2 , and three C n 's. We specify λ 2 so as to make C 9 to vanish: which gives There is then no remaining freedom, and the two remaining C n 's are given by where we have defined two new couplings ζ 1 and ζ 2 . To first order in four-derivative couplings, this gives the action S = S + δS: + βL 4 (RF µν F µν − 4R µν F µρ F ρ ν + R µνρσ F µν F ρσ ) + ζ 1 L 6 (F µν F µν ) 2 + ζ 2 L 6 F µν F νρ F ρσ F σµ .
Finally, the scaling of four-derivative couplings and the size of the gauge field that we employ in this work (see Section II) are λ GB ∼ β ∼ ζ 1 ∼ ζ 2 ∼ 2 and |A µ | ∼ , with small dimensionless The transformation of µ 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 µρ should transform the same way as E, which does not change under the rescaling of the gauge field. The conductivity transformation follows from the Kubo formula, which schematically relates σ ∼ [ρ, ρ]. Thus, we find the following relation between thermodynamic quantities and transport coefficients in our theory and the theory of Ref. [31]: T = e ω 1 T,μ = e ω 1 −2ω 2 µ,ρ = e 3ω 1 +2ω 2 ρ,σ = e ω 1 +4ω 2 σ, · · · . (A44) Furthermore, we solve (A21) to giveL With the above relations we can now use (A39) to solve for σ and ρ in terms of L, T , and µ, finding Plugging in L = 1 + λ GB /2, µ = μ and setting κ 2 5 = 2 to conform with the conventions used in this work (cf. Eq. (3.18)), we arrive at equations (3.35) and (3.39). Continuing further, we can also show agreement between the entropy, free energy and energy density from Ref. [31] 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. [31] after field redefinitions, which is calculated from the (higher-derivative) Wald formula [59,60], 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 [59,60].