Far-from-equilibrium dynamics of a strongly coupled non-Abelian plasma with non-zero charge density or external magnetic field

Using holography, we study the evolution of a spatially homogeneous, far from equilibrium, strongly coupled N=4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N}=4 $$\end{document} supersymmetric Yang-Mills plasma with a non-zero charge density or a background magnetic field. This gauge theory problem corresponds, in the dual gravity description, to an initial value problem in Einstein-Maxwell theory with homogeneous but anisotropic initial conditions. We explore the dependence of the equilibration process on different aspects of the initial departure from equilibrium and, while controlling for these dependencies, examine how the equilibration dynamics are affected by the presence of a non-vanishing charge density or an external magnetic field. The equilibration dynamics are remarkably insensitive to the addition of even large chemical potentials or magnetic fields; the equilibration time is set primarily by the form of the initial departure from equilibrium. For initial deviations from equilibrium which are well localized in scale, we formulate a simple model for equilibration times which agrees quite well with our results.


Introduction
The discovery of gauge/gravity duality (or "holography") has enabled the study of previously intractable problems involving the dynamics of strongly coupled gauge theories. 1 In the limit of large gauge group rank N c , and large 't Hooft coupling λ, the strongly coupled quantum dynamics of certain gauge field theories may be mapped, precisely, into classical gravitational dynamics of higher dimensional asymptotically anti-de Sitter (AdS) spacetimes [4][5][6]. Numerical studies of the resulting gravitational dynamics can shed light on poorly understood aspects of the quantum dynamics of strongly coupled gauge theories.
Using the simplest example of gauge/gravity duality, applicable to maximally supersymmetric SU(N c ) Yang-Mills theory (N = 4 SYM), this approach has been applied to a succession of problems of increasing complexity involving far from equilibrium dynamics.
In this paper, we extend previous work on the dynamics of homogeneous but anisotropic N = 4 SYM plasma [7][8][9]. We examine the influence on the equilibration dynamics of a non-zero global U(1) charge density, or a background magnetic field. Inclusion of these effects is motivated by the physics of relativistic heavy ion collisions [19][20][21]. Hydrodynamic modeling of near-central events clearly indicates that the baryon chemical potential µ B in the mid-rapidity region is significantly smaller than the temperature, but not by an enormously large factor at RHIC energies. 2 Hence, it is desirable to understand the sensitivity of the plasma equilibration dynamics to the presence of a baryon chemical potential and associated non-zero baryon charge density. Similarly, it is clear that large, but transient, electromagnetic fields are generated in heavy ion collisions. A growing body of work [24][25][26][27][28] suggests that electromagnetic effects may play a significant role despite the small value of the fine structure constant. Electromagnetic effects on equilibrium QCD properties are also under study using lattice gauge theory [29][30][31][32].
The large N c , strongly coupled N = 4 SYM plasma we study is, of course, only a caricature of a real quark-gluon plasma. But it is a highly instructive caricature which correctly reproduces many qualitative features of QCD plasma (such as Debye screening, finite static correlation lengths, and long distance, low frequency dynamics described by neutral fluid hydrodynamics). Moreover, in the temperature range relevant for heavy ion collisions, quantitative comparisons of bulk thermodynamics, screening lengths, shear viscosity, and other observables show greater similarity between N = 4 SYM and QCD than one might reasonably have expected [33,34]. Since the composition of a plasma depends on the chemical potentials, or associated charge densities, of its constituents, studying the dependence of the equilibration dynamics on a conserved charge density provides a simple means to probe the sensitivity of the dynamics to the precise composition of the non-Abelian plasma. This, in small measure, may help one gauge the degree to which N = 4 SYM plasma properties can be extrapolated to real QCD plasma. At the very least, strongly coupled N = 4 SYM theory provides a highly instructive toy model in which one may explore, quantitatively, non-trivial aspects of non-equilibrium gauge field dynamics. 3 The remainder of the paper is organized as follows. Section 2 summarizes necessary background material. This includes the coupling of an Abelian background gauge field to a U(1) subgroup of the SU(4) R global symmetry group of N = 4 SYM. This U(1) symmetry may be regarded as analogous to either the baryon number U(1) B or electromagnetic U(1) EM flavor symmetries of QCD. Turning on a background magnetic field implies an enlargement of the theory under consideration from N = 4 SYM to N = 4 SYM coupled to electromagnetism (which we abbreviate as SYM+EM). The combined theory is no longer scale invariant; this has important implications which we discuss. This section describes JHEP07(2015)116 the 5D Einstein-Maxwell theory which provides the holographic description of the states of interest, presents our coordinate ansatz (based on a null slicing of the geometry), and summarizes relevant portions of the holographic dictionary relating gravitational and dual field theory quantities. This section also records the reduced field equations which emerge from our symmetry specializations, describes the relevant near-boundary asymptotic behavior, and summarizes properties of the static equilibrium geometries to which our time dependent solutions asymptote at late times.
The following section 3 briefly describes our numerical methods, which are based on the strategy presented in ref. [18]. When studying states with a non-zero charge density (but no background magnetic field) appropriate numerical methods for asymptotically AdS Einstein-Maxwell theory are immediate generalizations of methods which have previously been found to work well for pure gravity. However, the inclusion of a background magnetic field induces a trace anomaly in the dual quantum field theory which, in the gravitational description, manifests in the appearance of logarithmic terms in the near-boundary behavior of fields. Such non-analytic terms degrade the performance of spectral methods, on which we rely, and necessitate careful attention to numerical issues. Section 3 also describes the specifics of our chosen initial data.
Results are presented in section 4. We focus on the evolution of the expectation value of the stress-energy tensor. We first discuss the sensitivity of the equilibration dynamics to features in the initial data and, in particular, examine the extent to which the evolution shows nonlinear dependence on the initial departure from equilibrium. We find that only disturbances in the geometry originating deep in the bulk, very close to the horizon, generate significant nonlinearities. This is broadly consistent with earlier work [8,9]. However, for a very wide variety of initial disturbances, including ones which generate extremely large pressure anisotropies, we find remarkably little nonlinearity in the equilibration dynamics, often below the part-per-mille level.
We then present comparisons of the equilibration dynamics as a function of the charge density or background magnetic field. We focus on comparisons in which the form of the initial departure form equilibrium and the energy density, or the equilibrium temperature, is held fixed while either the charge density or magnetic field is varied. These comparisons reveal surprisingly little sensitivity to the charge density, or magnetic field, even at early times when the departure from equilibrium is large.
We verify the late time approach to the expected equilibrium states, and extract the leading quasinormal mode (QNM) frequency from the late time relaxation. Quasi-normal mode frequencies extracted from our full nonlinear dynamics are compared, where possible, with independent calculations of QNM frequencies based on a linearized analysis around the equilibrium geometry. This provides a useful check on our numerical accuracy.
We define an approximate equilibration time based on the relative deviation of the pressure anisotropy from its equilibrium value, and examine the dependence of this time on charge density or external magnetic field. Once again, changes in this quantity are largest for initial disturbances which originate very close to the horizon, but the overall sensitivity of the equilibration time to the charge density or magnetic field is remarkably modest.
The final section 5 discusses and attempts to synthesize the implications of our results. We present a simple model of equilibration times, for initial disturbances which are well -3 -

JHEP07(2015)116
localized in scale, which agrees rather well with our numerical results (but becomes less accurate for disturbances localized extremely close to the horizon). We end with a few concluding remarks. 4 2 Ingredients

N = 4 SYM in an external field
We study maximally supersymmetric SU(N c ) Yang-Mills theory (N = 4 SYM) on four dimensional Minkowski space when the conserved current for a U(1) subgroup of the SU(4) R global symmetry group either (a) has a non-vanishing charge density, or (b) is coupled to a background Abelian gauge field describing a uniform magnetic field. The embedding of the U(1) symmetry is chosen such that the U(1) commutes with an SU(3) subgroup of the SU(4) R global symmetry.
The coupling to the external field has the usual form 5 where j α (x) is the conserved U(1) current normalized such that the four Weyl fermions of N = 4 SYM have charges {+3, −1, −1, −1}/ √ 3 and the three complex scalars have charge +2/ √ 3. The overall factor of 1/ √ 3 in these charge assignments has no physical significance, but is chosen so that the trace anomaly and electromagnetic beta function (induced when this current is gauged) have convenient coefficients, as will be seen below. 6 The background U(1) gauge field A ext α (x) we take to have the form with µ the chemical potential which, in equilibrium, will be conjugate to the charge density j 0 , and B the amplitude of a constant magnetic field pointing in the x 3 direction. Although it should be straightforward to study dynamics when both the charge density j 0 and magnetic field B are non-zero, in this paper we focus for simplicity on the cases of either a non-zero charge density with vanishing magnetic field, j 0 = 0 and B = 0, or non-zero magnetic field with vanishing charge density, j 0 = 0 and B = 0. With a non-zero magnetic field B in the x 3 direction, changes in the background gauge field under a translation in the x 1 or x 2 directions, or a rotation in the x 1 -x 2 plane, can be compensated by a suitable U(1) gauge transformation. Hence, the theory retains full spatial translation invariance as well as rotation invariance in the x 1 -x 2 plane. 4 As this paper neared completion, we learned of the somewhat related work by A. Buchel, M. Heller, and R. Myers [39]. These authors examine quasinormal mode frequencies in N = 2 * SYM and argue that, in this non-conformal deformation of N = 4 SYM, the longest equilibration times are largely set by the temperature with little sensitivity to other scales. 5 We use a mostly-plus Minkowski space metric, ηµν ≡ diag(−1, +1, +1, +1). 6 These charge assignments are 1/ √ 3 times those used in ref. [40]. Overall rescaling of these charge assignments has implications for the holographic description which are noted below in footnote 11.

JHEP07(2015)116
We will be interested in initial states which: (i ) have non-trivial expectation values T αβ (x) and j α (x) for the stress-energy tensor and U(1) current density, respectively; (ii ) are invariant under spatial translations as well as O(2) rotations in the x 1 −x 2 plane; and (iii ) are invariant under the SU(3) R subgroup of the SU(4) R global symmetry which commutes with our chosen U(1).
Since all N = 4 SYM fields transform in the adjoint representation of the SU(N c ) gauge group, the stress-energy and U(1) current expectation values both scale as O(N 2 c ) in the large N c limit. For later convenience, we define a rescaled energy density ε and charge density ρ, via SYM is a conformal field theory with a traceless stress-energy tensor. Adding a chemical potential µ introduces a physical scale, but does not modify the microscopic dynamics of the theory and hence does not affect the tracelessness of the stress-energy tensor. In contrast, introducing an external magnetic field does affect the microscopic dynamics and, in particular, generates a non-zero trace anomaly, 7 The trace anomaly generated by the external magnetic field implies that the theory is no longer scale invariant. For example, the ground state energy density, as a function of magnetic field, need not have the simple form of some pure number times B 2 . This will be seen explicitly below. The trace anomaly implies that there must be logarithmic dependence on a renormalization point. To interpret this dependence, it is appropriate to adopt the perspective that adding an external magnetic field means that the theory under consideration has been enlarged -it is now N = 4 SYM coupled to U(1) electromagnetism (SYM+EM). The complete action of the theory is the SYM action, minimally coupled to the U(1) gauge field, plus the Maxwell action for U(1) gauge field, The electromagnetic coupling e 2 (having been scaled out of covariant derivatives) appears as an inverse prefactor of the Maxwell action. We regard the electromagnetic coupling e 2 7 We define the external gauge field such that no factor of an electromagnetic gauge coupling appears in the interaction (2.1), in our U(1) covariant derivatives, or in the trace anomaly (2.5). The coefficient of − 1 4 F 2 µν in the trace anomaly (2.5) equals the EM beta function coefficient b0, given below in eq. (2.9). (Note that the sign of the trace anomaly depends on the metric convention in use.) as arbitrarily weak. Hence, quantum fluctuations in the U(1) gauge field are negligibly small, allowing us to view the EM gauge field as a classical background field. 8 However, just as in QED, fluctuations in the SYM fields which are electromagnetically charged will cause the electromagnetic coupling e 2 to run with scale. The associated renormalization group (RG) equation for the inverse coupling has the usual form, with the one-loop beta function coefficient 9 Here, q α f = (3, −1, −1, −1)/ √ 3 and q a s = (2, 2, 2)/ √ 3 are the charge assignments of the four Weyl fermions and three complex scalars, respectively. Integrating this renormalization group equation leads, as usual, to with the RG invariant scale Λ EM denoting the Landau pole scale where the (one loop approximation to the) electromagnetic coupling diverges. The total stress-energy tensor derived from the combined action (2.6) will equal the N = 4 SYM stress-energy tensor, augmented with minimal coupling terms to the EM gauge field, plus the classical Maxwell stress-energy. An essential point, however, is that while the total stress-energy tensor is well-defined, partitioning the stress-energy tensor into separate SYM and EM contributions is inherently ambiguous, as the individual pieces depend on the renormalization point. We define In an arbitrary background SU(4) gauge field, the divergence of the SU(4)R current acquires an anoma- This anomaly, when specialized to our chosen U(1) subgroup, is proportional to the sum of the cubes of our fermion charges and is non-zero, α (q α To make the combined SYM+EM theory well defined, one could add to the theory additional fermions, charged under the U(1) but with no SYM interactions, which would cancel this U(1) anomaly. As we are not concerned with quantum fluctuations in the U(1) gauge field, the presence of this U(1) anomaly (in the absence of compensating spectators) is irrelevant for our purposes. 9 A non-renormalization theorem in supersymmetric N = 4 SYM implies that the short distance behavior of the current-current correlation cannot depend on the 't Hooft coupling λ [41]. This implies that the leading EM beta function coefficient b0 does not depend on λ, and hence may easily be evaluated in the λ → 0 limit. 10 Note that T αβ EM (µ) is not the metric variation of some renormalized EM action (whose separation from the total action would not be well-defined). Rather, eq. (2.12) is simply defining T αβ EM (µ) as the classical EM stress-energy tensor multiplied by the scale-dependent inverse EM coupling. The partitioning (2.11) of the stress-energy tensor puts all quantum corrections other than the running of the EM coupling into the SYM contribution ∆T αβ SYM (µ). The scale dependence must, of course, cancel between the two terms because the total stress-energy tensor is a physical quantity. Therefore, the scale dependence in the SYM contribution to the stress-energy must simply compensate the known running of the inverse electromagnetic coupling (2.8) in the Maxwell stress-energy tensor (2.12), (2.14) Specializing to zero temperature states in a constant static magnetic field B, the scale dependence (2.14) plus dimensional analysis implies that the SYM contribution to the ground state energy density is a non-analytic function of magnetic field, with c 0 some pure number. (Here and henceforth, when considering physics in a nonzero magnetic field ε(µ) ≡ ∆T 00 SYM (µ)/κ denotes the SYM portion of the rescaled energy density.) In the second form of eq. (2.15), the analytic term has been absorbed by defining a scale dependent "fiducial" magnetic field amplitude, (2.16) Note that the ground state energy acquires a simple quadratic form when the renormalization point is chosen to scale with the magnetic field, ε(|B| 1/2 ) = c 0 B 2 . Our numerically determined value for the coefficient c 0 is given below in eq. (2.68). When considering low temperature physics in a background magnetic field, T 2 |B|, it is natural to choose a renormalization point µ = O(|B| 1/2 ), as this is the relevant scale which cuts off long range fluctuations in the charged SYM fields. We will employ two choices for the renormalization point. One choice is µ = 1/L, with L the AdS curvature scale (discussed below); this choice is computationally convenient but not physically significant. We will also report and discuss results with µ = |B| 1/2 . For later convenience, we define abbreviations for the (rescaled) energy density evaluated at these two renormalization points,

Holographic description
The holographic description of SYM states, within our sector of interest, in the limit of large N c and large 't Hooft coupling λ, is given by classical Einstein-Maxwell theory on 5-dimensional spacetimes which are asymptotically AdS 5 [42]. The 5D bulk action is with G 5 ≡ π 2 L 3 /N 2 c the 5D Newton gravitational constant, Λ ≡ −6/L 2 the cosmological constant, and L the AdS curvature scale. 11 Setting to zero the variation of the action with respect to the metric gives the Einstein equation,  [40,42].) However, as stated above, in this paper we consider solutions with non-zero chemical potential µ or non-zero magnetic field B, but not both µ and B nonzero. For such solutions, the Chern-Simons term makes no contribution to the dynamics and hence may be neglected.
As usual in holography, the expectation value T αβ (x) of the stress-energy tensor is determined by the subleading near-boundary behavior of the 5D metric G M N . The leading near-boundary behavior of the bulk gauge field A M will be fixed by our chosen external U(1) gauge field (2.2), while the expectation value j α (x) of the U(1) current density is determined by the subleading near-boundary behavior of the bulk gauge field. The precise relations will be shown below.
Following ref. [18], we choose a coordinate ansatz, based on generalized Eddington-Finklestein (EF) coordinates, which is natural for gravitational infall problems. The metric has the general form where r is the bulk radial coordinate and x ≡ {x α }, α = 0, · · ·, 3, denotes the four remaining spacetime coordinates. The spacetime boundary lies at r = ∞; the {x α } may be regarded as coordinates on the spacetime boundary where the dual field theory "lives". Curves of varying r, with x held fixed, are radially infalling null geodesics, affinely parameterized by r. The one-form w ≡ w α dx α (which is assumed to be timelike) depends only on x, not on r. These infalling coordinates remain regular across future null horizons. The form of the ansatz (2.20) remains invariant under r-independent diffeomorphisms, The coefficient of the Maxwell action may, of course, be set to an arbitrary value by suitably rescaling the bulk gauge field AM . However, as the on-shell variation of the gravitational action with respect to the boundary value of the gauge field defines the associated current, such rescaling changes the normalization of the U(1) current in the holographic description. It will be seen below that the coefficient of the Maxwell term in our action (2.18) is correctly chosen so that the U(1) current normalization is consistent with our previous charge assignments. If charge assignments are chosen, for example, to be larger by a factor of √ 3, then either the Maxwell term in the action (2.18) must be multiplied by a factor of 3, or else one must regard the boundary value of the bulk gauge field as equaling √ 3 times the QFT gauge field (and the charge density in the bulk theory as equal to the QFT charge density divided by √ 3), as was done in ref. [40]. as well as radial shifts (with arbitrary x dependence), We use the diffeomorphism freedom (2.21) to transform the timelike one-form w to the standard form −dx 0 (or w α = −δ 0 α ). Our procedure for dealing with the radial shift invariance (2.22) is discussed below in subsection 2.6.
We are interested in geometries which, at large r, asymptotically approach (the Poincaré patch of) AdS 5 . This will be the case if g αβ (x, r) approaches η αβ as r → ∞, with η αβ ≡ diag(−1, 1, 1, 1) the usual Minkowski metric tensor. Demanding that the metric and bulk gauge field satisfy the Einstein-Maxwell equations, one may derive the near-boundary asymptotic behavior of the fields. Using radial gauge, A r = 0, for the bulk gauge field, and a suitable choice of the radial shift (2.22) (which eliminates O(1/r) terms in g αβ ), one finds that for solutions of interest, the metric and gauge field have asymptotic expansions of the form The coefficient h (4) αβ of the logarithmic term in the metric is only non-zero when there is an external EM field, For a constant magnetic field in the x 3 direction, h . The subleading asymptotic coefficients g (4) αβ (x) and A (2) α (x) cannot be determined solely from a near-boundary analysis of the field equations, and depend on the form of the solution throughout the bulk. However, asymptotic analysis does show that 3 i=1 g The subleading metric coefficients g (4) αβ (x) and h (4) αβ (x) encode the expectation value of the SYM stress-energy tensor [43,44]. The appropriate holographic relation is T µν = κ g (4) µν − η µν tr ( g (4) ) + [ln(µL) , and C is an arbitrary renormalization-scheme dependent constant. 13 We adopt a specific value, In Fefferman-Graham (FG) coordinates, for which ds 2 ≡ (L 2 /ρ 2 ) g αβ (x, ρ) dx α dx β + dρ 2 , one has g αβ (x, ρ) ∼ η αβ + g (4) αβ (x) + h (4) αβ ln L ρ ρ 4 + O(ρ 6 ln ρ) as ρ → 0. Eq. (2.26) gives the relation between the subleading asymptotic metric coefficients in our infalling EF coordinates and FG coordinates. 13 To perform the required holographic renormalization one must add a counterterm depending logarithmically on the UV cutoff. (See, for example, refs. [44][45][46].) As always, such a logarithmic counterterm comes with an inevitable finite ambiguity. which will make the subsequent explicit expression (2.42a) for the energy density as simple as possible.

Symmetry specialization
As noted earlier, we are interested in studying solutions of Einstein-Maxwell theory which are spatially homogeneous. This implies that all metric functions depend only on x 0 and r. The arbitrary function λ in the residual radial shift diffeomorphism (2.22) will depend only on x 0 . Henceforth, for convenience, we will use v as a synonym for x 0 ; v is a null time coordinate. (In other words, v = const. surfaces are null slices of the geometry.) At the boundary, v coincides with the time t of the dual field theory. We also impose invariance under O(2) rotations in the x 1 -x 2 plane. This implies that only the g 00 , g 03 , g 33 , and g 11 = g 22 components of g αβ are non-zero. Our Einstein-Maxwell theory (without a Chern-Simons term) is also invariant under spatial parity, or x 3 → −x 3 reflections, and for simplicity we will also impose parity invariance. This requires the vanishing of g 03 .
For the bulk gauge field, the choice of radial gauge, A r = 0, plus our imposed symmetries imply that The corresponding bulk field strength, which is what appears in the field equations, can have a constant (x and r independent) magnetic field plus a radial electric field, with all other components vanishing.
14 As in ref. [40], one can also use a comparison of holographic and QFT evaluations of the U(1) anomaly to confirm that the U(1) current normalizations are consistent.

JHEP07(2015)116
As in ref. [18], it is convenient to rename the non-vanishing metric components as where A, B, and Σ are functions of v and r. The resulting line element is Henceforth, A will always denote the metric function multiplying dv 2 (times −1/2), not the bulk gauge field. The function Σ is the spatial scale factor (with Σ 3 dx dy dz the spatial volume element), while B characterizes the spatial anisotropy (which should not be confused with the magnetic field amplitude B).
The radial derivative ∂ r is a directional derivative along infalling radial null geodesics. It proves convenient to define a corresponding directional derivative along outward radial null geodesics, The field equations which result from varying the action (2.18), inserting the above symmetry specializations, and re-expressing v-derivatives in terms of the d + modified time derivative (2.35), take a remarkably compact form. The Einstein equations are: where primes denote radial derivatives, h ≡ ∂ r h. As discussed in ref. [18], the anisotropy function B encodes the essential propagating degrees of freedom. The functions Σ and A may be regarded as auxiliary fields, determined by solving eqs. (2.36a) and (2.36b) using data on a single time slice. Information about the time evolution of B is contained in equation (2.36c). Equations (2.36d) and (2.36e) may be viewed as boundary value constraints -if they hold at one value of r, then the other equations ensure that these equations hold at all values of r.
Maxwell's equations reduce to the statements that neither the magnetic field B, nor the radial electric flux density E Σ 3 , have any radial or temporal variation. In other words, B = const., as already indicated in (2.32), and The form (2.37) of the radial electric field simply reflects Gauss' law in 4+1 dimensions, combined with charge conservation and spatial translation invariance, which imply that ρ cannot have any temporal or spatial variation. The bulk gauge field A M does not appear in the field equations (except via the field strength), but one may choose to regard A M as satisfying the radial gauge condition, A r = 0, plus the condition that the time component A v vanish at the horizon. This fixes the residual r-independent gauge freedom which remains after imposing radial gauge. With these choices, the chemical potential µ is the boundary value of A v in the late time (v → ∞) equilibrium limit. Equivalently (in radial gauge), the chemical potential µ equals the difference between the boundary and horizon values of A v , in the equilibrium geometry. This coincides with the line integral of the radial electric field from horizon to boundary, which gives the work needed to move a unit charge from the boundary to the horizon. As usual, the charge density and the chemical potential are thermodynamically conjugate. One may consider the chemical potential µ to be a function of the (rescaled) charge density ρ, or vice-versa. The choice of perspective ("canonical" vs. "grand canonical") has no bearing on the dynamics.

Asymptotic analysis
Asymptotic analysis of these equations is straightforward. We impose a flat boundary geometry with the requirement that lim r→∞ g αβ (x, r) = η αβ , implying for our renamed metric functions. Solutions to Einstein's equations (2.36) with this leading behavior may be systematically expanded in integer powers of 1/r and (for non-zero magnetic field) logarithms of r. One finds: The constant a 4 and the function b 4 (v) cannot be determined just using asymptotic analysis, and the radial shift λ(v) is completely arbitrary. The coefficient a 4 encodes the energy density which, due to homogeneity, cannot vary in time, while b 4 (v) encodes the anisotropy in the spatial stress. Using the holographic relation (2.25) and our convention (2.27) for defining the stress-energy tensor, one finds

Scaling relations
Consider independent rescaling of the boundary and radial coordinates, also satisfy the Einstein equations (and our asymptotic conditions) with rescaled parameters The subleading asymptotic coefficients a 4 and b 4 become Using the holographic relation (2.42) for the stress-energy expectation, one finds that with a rescaled renormalization point µ ≡ α µ. If α = γ, then these transformations are just a trivial rescaling of all quantities according to their dimension. But transformations with α = γ are non-trivial and scale bulk and boundary quantities by different amounts. In particular, transformations with α = 1 but γ = 1 rescale the AdS curvature scale L without affecting the boundary coordinates or boundary parameters (B, ρ, or µ), showing that the value of L has no physical significance (in the large N c , large λ limit for which classical gravity provides the dual description). This illustrates, explicitly, the independence of the boundary field theory on the AdS curvature scale L.

Apparent horizon
With a non-zero homogeneous energy density, the dual geometries of interest will have an apparent horizon at some radial position, r = r h (v) [18]. Since we are investigating nonequilibrium dynamics, one might expect the horizon position to change significantly before ultimately settling down as the system equilibrates. However, as illustrated in figure 1 it is possible, and very convenient, to use the residual radial shift diffeomorphism freedom (2.22) to place the apparent horizon at a fixed radial position, A short exercise [18] shows that the condition for an apparent horizon to be present at r =r h is that this location be a zero of the modified time derivative of the spatial scale factor, This condition serves to fix the radial shift λ(v). It is convenient to regard this condition as a combination of a constraint on initial data (implemented by finding the radial shift λ(v 0 ) which is needed to satisfy (2.49) at some initial time v 0 ) together with the requirement that the horizon be time-independent, ∂r h /∂v = 0, which requires that the time derivative of d + Σ vanish at the apparent horizon. Evaluating this condition, and using the Einstein equations (2.36d), (2.36e) to simplify the result, determines the value of the metric function A at the horizon, For the metric to be non-singular on and outside the apparent horizon, the spatial scale factor Σ must be non-vanishing for r ≥r h .

Equilibrium solutions
Given some initial non-equilibrium state of the system, the dynamical evolution should asymptotically approach a thermal equilibrium state. In the gravitational description, this implies that the geometry should, at late times, approach some static black brane solution. The specific black brane solution will depend on the values of the conserved energy and charge densities in the chosen initial state, and on the value of the background magnetic field.

JHEP07(2015)116
Schwarzschild. For initial states with vanishing charge density and magnetic field, the bulk geometry will equilibrate to the 5D AdS-Schwarzschild black brane solution. A standard form of this metric is The radial coordinater should not be confused with our Eddington-Finklestein coordinate r. The zero of U (r) determines the horizon location, 53) and the horizon temperature [given by (2π) −1 times the surface gravity at the horizon] is proportional to the horizon radius, In our infalling EF coordinates (2.34), this AdS-Schwarzschild solution is described by Using the holographic relation (2.42), one sees that the parameter m is related to the (rescaled) equilibrium energy density ε ≡ T 00 /κ via Reissner-Nordstrom. If the initial state has a non-zero charge density but vanishing magnetic field, then the bulk geometry will equilibrate to a 5D Reissner-Nordstrom (RN) black brane [42]. This metric may be written in the form (2.51), with The charge density ρ of the Reissner-Nordstrom brane is bounded from above by the extremal charge density ρ max , given by The relation (2.56) between the energy density and the mass parameter m is unchanged. Hence, the extremal charge density ρ max = 4 3 ε 3/4 . It is convenient to express ρ in terms of the fraction x of the extremal charge density, . Also shown are the high and low temperature asymptotic forms (dashed lines). In the high temperature regime, πT ρ 1/3 , the curve approaches the Schwarzschild result, ε = 3 4 (πT ) 4 . In the low temperature (or near extremal) regime, πT The horizon radiusr h is given by the outermost positive root of U (r); explicitly, The horizon radius (divided by L) varies from m 1/4 down to ( 1 3 m) 1/4 as x varies from 0 to 1. The horizon temperature T h is given by The horizon temperature decreases with increasing charge density, and vanishes as the charge density approaches the extremal value (or x → 1). From the perspective of the dual field theory, for any given value of the charge density there is a lower bound on the energy density, ε min (ρ), which must be a monotonically increasing (and convex) function of ρ. This implies that for any given value of the energy density, there will be a maximum charge density, corresponding to a ground state with vanishing temperature. The equilibrium chemical potential µ, thermodynamically conjugate to the charge density ρ, is given by Physically distinct non-extremal solutions may be labeled by the value of one dimensionless ratio such as ε/|ρ| 4/3 [or (πT ) 4 /|ρ| 4/3 or ε/(πT ) 4 , etc.]. Figure 2 shows a log-log plot of the curve representing these solutions in the plane of ε/|ρ| 4/3 and (πT ) 4 /|ρ| 4/3 . In our infalling EF coordinates, the RN black-brane solution is described by and B(r) = 0.

JHEP07(2015)116
Magnetic branes. When the magnetic field is non-zero, the bulk geometry will equilibrate to a stationary magnetic black brane solution. These solutions are not known analytically, but have been studied numerically [40,43]. In our infalling coordinates, the solutions satisfy the static specialization of eqs. (2.36). 15 The near-boundary behavior of asymptotically AdS 5 solutions is given by the expansions (2.41) (but with no time dependence). The extremal, zero-temperature magnetic brane solution interpolates smoothly between AdS 3 × R 2 near the horizon and AdS 5 near the boundary. In our infalling coordinates, a series in fractional powers of δr ≡ r−r h describes deviations from the AdS 3 × R 2 geometry near the horizon, The constant η cannot be determined from a near-horizon analysis and must be suitably adjusted after integrating eqs. (2.64) to obtain the desired boundary geometry. There is a single extremal magnetic brane solution, modulo the rescaling transformations (2.43)-(2.45) (which relate solutions with any non-zero values of the magnetic field B and curvature scale L), and radial shift diffeomorphisms (2.22).
For non-extremal solutions (with non-zero B but vanishing ρ), metric functions near the horizon have power series expansions in δr ≡ r−r h of the form The coefficient a 0 is proportional to the horizon temperature, The other two undetermined constants, s 0 and β, which control the horizon values of the spatial scale factor and the anisotropy function, must be suitably adjusted after integrating 15 The resulting equations may be written explicitly as: eqs. (2.64) to select solutions which have the desired near-boundary behavior (with an isotropic boundary metric). If B 2 T 4 , then the resulting magnetic brane geometry is a small perturbation away from the Schwarzschild black brane (2.55), while if B 2 T 4 then the geometry may be regarded as interpolating between the BTZ black brane (× R 2 ) near the horizon and AdS 5 near the boundary [40].
There is a one parameter family of non-extremal solutions, modulo the rescaling transformations (2.43)-(2.45) (and radial shift diffeomorphisms). Physically distinct solutions may be labeled by the value of the dimensionless ratio ε B /B 2 [or (πT ) 4 /B 2 or ε B /(πT ) 4 , etc.]. The left panel of figure 3 shows a log-log plot of our numerically determined curve representing these solutions in the plane of ε B /B 2 and (πT ) 4 /B 2 . Extrapolating our lowest temperature numerical results to zero temperature, we find estimates of for the coefficient c 0 or the equivalent fiducial scale B * defined by eqs. (2.15) and (2.16).
If one chooses to measure energy density and magnetic field in units set by the curvature scale L, then one may traverse the one-parameter family of magnetic brane solutions by varying |B|L 2 for a fixed value of ε L L 4 (or vice versa). The holographic relation (2.42) shows that these curvature scale dependent quantities are related to the intrinsic dimensionless In all cases, the radial shift λ has been suitably adjusted to fix the horizon radius at u = 1. Note that the horizon value of the anisotropy function is not a monotonic function of magnetic field at fixed ε L , but is monotonic when ε B is held fixed.
For these non-extremal magnetic brane solutions, the anisotropy function B(r) increases (and the scale factor Σ decreases) smoothly as one moves inward from the boundary toward the horizon. The right panel of figure 4 shows a similar set of solutions with increasing magnetic field, but now with the energy density at the scale |B| 1/2 held fixed, ε B = 8 L −4 . With the physical parameter ε B held fixed, the horizon value of the anisotropy function is now monotonically increasing with magnetic field. The temperatures of these solutions (in order of increasing B) are given by πT L = 1.807, 1.797, 1.738, 1.620, and 1.433.

Computational strategy
We apply the computational strategy presented in ref. [18] to our case of homogeneous isotropization in Einstein-Maxwell theory. For convenience, we choose units in which the AdS curvature scale L = 1.
Required initial data, on some v = v 0 time slice, consists of an initial choice for the anisotropy function B(v 0 , r) and the radial shift λ(v 0 ), along with chosen (time independent) values of the energy density ε, charge density ρ, and magnetic field B. As noted above, the holographic relation (2.42a) shows that fixing the energy density ε at the scale -19 -JHEP07(2015)116 µ = 1/L is equivalent to fixing the asymptotic coefficient a 4 . Our choices for the initial anisotropy function will be detailed below in subsection 3.3.
Given a set of initial data, the linear second order radial ordinary differential equation (ODE) (2.36a) may be integrated to find the spatial scale factor Σ(v 0 , r). The two leading terms in the asymptotic behavior (2.41a) provide the integration constants needed to specify uniquely the desired solution. Next, one solves eq. (2.36d), a linear first order radial ODE for d + Σ. The near-boundary asymptotic behavior of this function is and homogeneous solutions to eq. (2.36d) behave as r −2 near the boundary. Hence, the chosen value of the energy density ε uniquely specifies the desired solution. With B, Σ, and d + Σ determined on the v 0 time slice, one next solves eq. (2.36c), a linear first order radial ODE for d + B. The desired asymptotic behavior of this function is d + B ∼ −2b 4 r −4 + O(r −5 ) while homogeneous solutions to eq. (2.36c) behave as r −3/2 near the boundary. So the needed integration constant corresponds to requiring the absence of any such homogeneous solution. Finally, one solves the second order linear ODE (2.36b) to determine A(v 0 , r). Homogeneous solutions are linear or constant functions of r. From the asymptotic behavior (2.41b), one sees that the value of the radial shift λ(v 0 ) fixes the coefficient of the homogeneous solution linear in r and provides one of the two needed boundary conditions. The second boundary condition, needed to fix the constant homogeneous solution, is provided by the horizon stationarity condition (2.50), which determines the value of A on the apparent horizon.
Having solved for the modified time derivative d + B(v 0 , r), and A(v 0 , r), one reconstructs the ordinary time derivative of the anisotropy function via (3.1) The time derivative of the radial shift, ∂ v λ(v 0 ), is extracted from the asymptotic behavior (2.41b) of A by evaluating the r → ∞ limit of A − 1 2 (r+λ) 2 . These time derivatives of B and λ provide the information needed to advance in time. Using a standard numerical integration scheme, one takes a small step forward in time, advancing v to v 0 +∆v for some timestep ∆v. Repeating this process, one progressively determines the metric functions on a sequence of equally spaced time slices, v = v k ≡ v 0 + k ∆v. On each time slice, the asymptotic coefficient b 4 (v), needed to determined the stress tensor (2.42), is extracted from the large r behavior of the anisotropy function B. (In the presence of a non-zero magnetic field, one extracts b 4 from the large r limit of a subtracted, rescaled version of the anisotropy function which removes the leading logarithmic piece in eq. (2.41c). This is detailed in the next subsection.)

Numerical methods
We use an inverted radial coordinate u = 1/r, and arbitrarily choosē as our apparent horizon location. This makes our computational domain a fixed radial interval, 0 ≤ u ≤ u h ≡ 1. We use a 4th order Runge-Kutta method (described in ref. [18]) -20 -

JHEP07(2015)116
for time integration. This requires four integrations of our radial ODEs per time step, but yields much better accuracy, for a given timestep ∆v, than a lower order method.
To integrate the radial ODEs (2.36a)-(2.36d), we have used both traditional shortrange finite difference approximations, and spectral methods [47]. In the latter approach, one implicitly represents the radial dependence of functions as a (truncated) series of Chebyshev polynomials. Explicitly, functions are represented by the list of their values on a discrete, finite collocation grid consisting of the points and derivatives are represented by (dense) M × M matrices acting on the finite list of function values. The (truncated) spectral expansion converts each ODE into a straightforward linear algebra problem. Boundary conditions are simply encoded into the first or last rows of the resulting matrix [47].
Although there are subtleties (described momentarily) in applying spectral methods to our problem, we have found the use of spectral methods to be clearly superior to finite difference approximations, yielding both faster computation and more accurate results. Using an M point discretization of the computation domain, short-range finite difference methods have errors which scale as an inverse power of M while spectral methods, in favorable cases, produce errors which fall exponentially with increasing M .
Spectral methods presume that one is approximating functions which are regular and well-behaved on the computational domain. However, our radial ODEs have regular singular points at u = 0 or r = ∞ (due to the r 2 growth of the scale factor near the boundary), and our functions Σ, d + Σ, and A all diverge at the u = 0 endpoint. Therefore, we define subtracted functions in which the known singular near-boundary behavior is removed. To minimize loss of numerical precision in spectral approximations of derivatives, it is also helpful to rescale the subtracted functions so that they do not vanish faster than linearly as u → 0. If the magnetic field is zero, so no logarithmic terms are present in the nearboundary behavior (2.41), then our subtracted functions are analytic in a neighborhood of the u ∈ [0, 1] radial interval and spectral methods converge exponentially. With a nonzero magnetic field, logarithmic terms are unavoidably present, showing that u = 0 is a branch point of the metric functions. This degrades convergence of the spectral series, leading to power-law convergence at a rate which depends on the behavior of the leading non-analyticity. Consequently, it is desirable to subtract logarithmic terms to as high an order as is practicable. We chose to introduce subtracted/rescaled functions (denoted with a subscript 's') via:

JHEP07(2015)116
All our numerical work is performed using these subtracted/rescaled functions; we directly solve for Σ s , (d + Σ) s , (d + B) s , and A s . 16 The expressions (3.4) are used to reconstruct the original functions when needed. We also use a subtracted/rescaled anisotropy function B s (r), introduced via This removes leading logarithmic terms and introduces a convenient rescaling. Henceforth, B s will be referred to as the subtracted anisotropy function.
In the above subtractions, the series in 1/r multiplying each logarithm are just truncated expansions of (r+λ) −k for k = 2, 3 or 4. The choice to truncate these binomial series was arbitrary, but we found that our numerics were sufficiently well-behaved with the above subtractions. At higher orders in 1/r, additional terms appear which involve the asymptotic coefficient b 4 and its (a-priori unknown) time derivatives, as well as higher powers of logarithms.

Initial data
As indicated above, one must select the value of the energy density (or asymptotic coefficient a 4 ) and the initial value of the radial profile of the anisotropy function, B(v 0 , r). And one must choose the value of the magnetic field B or charge density ρ. For the charged case, with vanishing magnetic field, physics can only depend on the dimensionless combination of charge and energy densities ρ/ε 3/4 , so a (positive) value of ε may be chosen arbitrarily without loss of generality. Given a choice of ε, possible values of the charge density ρ are limited by the extremality bound |ρ| ≤ ρ max = 4 3 ε 3/4 . For the initial anisotropy function, we chose to focus on Gaussian profiles. In the λ = 0 frame, We investigate the dependence of results on the parameters of the Gaussian (amplitude A, width σ, and mean position r 0 ) in the first part of section 4. 17 Motivated by the fact that in our coordinates lines of varying r, at fixed v, are radially infalling geodesics, we will refer to this initial Gaussian as a "pulse" of initial anisotropy.
For the magnetic case, as discussed above, the breaking of scale invariance implies the presence of logarithmic terms in the asymptotic behavior of the anisotropy function. We simply add the log terms shown in eq. (3.5) to the Gaussian (3.6). With vanishing radial shift, λ = 0, our chosen initial anisotropy function takes the form Arguably, a more natural choice for the magnetic case initial data might be to add a Gaussian to the full equilibrium solution for the anisotropy function in the chosen magnetic

JHEP07(2015)116
field. This could be seen as nicely paralleling the charged case (in which the equilibrium solution has vanishing anisotropy). Nevertheless, we will stick with our somewhat arbitrary choice (3.7), which is an adjustable initial pulse added to the correct asymptotics. As will be seen, the Gaussian pulse will nearly always be the dominant portion of the deviation from equilibrium and the driving force of the resulting anisotropy in the boundary stress. We doubt that differing choices in the precise form of the slowly varying function to which the Gaussian pulse is added would impact, in any significant way, the characteristic equilibration times or other significant features of the results presented below. After choosing the initial anisotropy function (in the λ = 0 frame) the initial value of the radial shift, λ(v 0 ), is adjusted to fix the location of the apparent horizon, as discussed in section 2.6.
It should be noted that, in all cases (charged, uncharged, magnetized) it is quite possible to select physically inconsistent initial data. This happens when the initial anisotropy, for a given energy density, is so large that no apparent horizon shields a coordinate singularity produced by a vanishing scale factor Σ. This sets a natural limit on the amplitude of the initial pulse which is meaningful to study.

Neutral plasma
Before presenting results for equilibration in charged or magnetized plasmas, we first discuss general features of the time evolution in the uncharged, unmagnetized case and examine the sensitivity of results to specific features in the initial data. As noted above, we choose a Gaussian profile (3.6) for the initial anisotropy function, with an adjustable amplitude, width, and mean position. Typical evolution of our subtracted/rescaled anisotropy function, B s (v, u) ≡ u −3 B(v, u), is shown on the left in figure 5. One sees the initial pulse profile on the back side of the plot at v = 0. The figure clearly shows the influence of the pulse propagating outward and reflecting off the boundary at u = 0. The outgoing portion of the pulse essentially propagates along an outgoing radial null geodesic which, in our coordinates, near the boundary are 45 • lines at constant values of u + v. The influence of the anisotropy pulse, after the reflection, largely falls through the horizon along an ingoing radial null geodesic which is instantaneous in our null time coordinate v. The asymptotic coefficient b 4 , which equals the slope of B s at u = 0, controls the anisotropy in the stress tensor, with ∆P/κ = 3 b 4 . Hence, the reflection of the pulse off the boundary directly produces the pressure anisotropy ∆P in the boundary theory. The time dependence of the relative pressure anisotropy, 18 defined as ∆P/(κ ), is plotted on the right in figure 5. . The (rescaled) energy density ε is used to the set the scale for time. One sees that the effect of the initial Gaussian pulse propagates outward, essentially along an outgoing radial null geodesic, reflects off the boundary, and then largely falls through the horizon (along an ingoing radial null geodesic which is instantaneous in v). After one "bounce", the anisotropy rapidly approaches zero. Right: the corresponding relative pressure anisotropy, ∆P/κε ≡ 1 2 (T 11 + T 11 − 2T 33 )/κε, induced in the boundary field theory, as a function of time. Note how the peaks of the pressure anisotropy correspond directly to the reflection of the anisotropy pulse off the boundary. late times, when the departure from equilibrium is small, the evolution should be well described by a linearized approximation to the full nonlinear dynamics. The linearized dynamics of infinitesimal perturbations away from equilibrium may be represented as a sum of quasinormal modes (QNM), which are eigenfunctions of the linearized dynamics with complex frequencies, φ(t) = Re(A e −iωt ) with Im ω < 0. The lowest quasinormal mode (for which −Im ω is minimal) dominates the late time approach to equilibrium. For the Schwarzschild black brane, quasinormal mode frequencies have been previously evaluated by Starinets [48]. From the late time behavior of our full nonlinear evolution, it is straightforward to extract an estimate of the lowest quasinormal mode frequency. Comparing with the independent results of ref. [48] provides a useful test of the accuracy of our numerics. Fitting the late time (4 v ε 1/4 25) portion of our calculated pressure anisotropy to a decaying, oscillating exponential, |A|e (Im ω)v cos[(Re ω)v + φ], yields an estimate of the lowest QNM frequency ω which agrees with ref. [48] to five digits, ω/(πT ) ≈ 3.11946 − 2.74663 i.
We will see the same vanishing of the anisotropy function at late times for the case of charged plasmas, whose gravitational duals equilibrate to an isotropic Reissner-Nordstrom black brane solution. For the magnetic case, however, at late times there is a non-zero profile for the anisotropy function, reflecting the spatial anisotropy of equilibrium magnetic brane solutions.
We now turn to an examination of the dependence of the pressure anisotropy on the parameters of the initial Gaussian (3.6). Of particular interest will be the dependence  Figure 6. Left: the pressure anisotropy ∆P/κε produced by the same initial pulse shown in figure 5 (blue curve), overlaid with the pressure anisotropy produced by an initial pulse with half the amplitude (purple curve). Halving the initial amplitude roughly halves the induced pressure anisotropy. Right: the "nonlinearity" (NL), defined as the difference between the pressure anisotropy produced by the larger pulse and twice the pressure anisotropy produced by the half amplitude pulse. of the response on the amplitude and position of the initial pulse. Less interesting is the dependence on the width of the pulse, which affects the duration of the reflection off the boundary (and also produces changes more naturally associated with the position of the pulse).
In figure 6, we compare the pressure anisotropies created by the pulse shown in figure 5, and an otherwise identical pulse with half the amplitude. Visually, one sees that the smaller amplitude pulse produces roughly half the pressure anisotropy as does the larger pulse, but with a virtually identical time course. The peak pressure anisotropy (divided by energy density), for both pulses is over 4, significantly larger than unity. Hence both initial pulses represent large departures from equilibrium. Given the highly nonlinear nature of the Einstein equations, one might have expected to see clear signs of nonlinearity in the dependence of the pressure anisotropy on the initial pulses. However, even for these pulses producing large departures from equilibrium, the amplitude of the peaks in the induced pressure anisotropy are nearly linear in the amplitude of the initial Gaussian pulse.
The right hand panel of figure 6 makes this comparison quantitative. This shows the nonlinearity (NL) defined as the difference between the pressure anisotropy ∆P/κε of the larger initial pulse and twice the pressure anisotropy produced by the halved initial pulse. Compared to the pressure anisotropies themselves, the relative size of the nonlinearity is roughly one part in 10 7 . This suggests that the dynamics, as probed by these initial pulses, are surprisingly close to a linear dynamical system.
In asymptotically AdS gravitational solutions, deviations of the geometry from that of pure AdS space necessarily vanish as one approaches the boundary. Hence, one might expect nearly linear dynamics to be evident for initial pulses which are localized sufficiently close to the boundary, while anticipating much larger nonlinearities for initial pulses localized closer to the horizon. To test this expectation, we used a very large Gaussian -the same initial Gaussian profile which generated figure 5 but with the amplitude increased 19 19 For the chosen values of position and width of the Gaussian, plus energy density ε = 3 4 L −4 , this amplitude is close to the upper limit set by demanding the existence of an apparent horizon.  2. × 10 -7 NL Figure 7. Comparisons of initial anisotropy functions (left column), induced pressure anisotropy (middle column), and nonlinearity (right column) defined as in figure 6, for a series of five Gaussian initial anisotropy functions differing only in their depth in the bulk. From top to bottom, the mean position of the initial pulse, in the λ = 0 frame, is r 0 = {4, 2, 1, 1 2 , 1 4 }. In all cases the energy density is ε = 3 4 L −4 . The plots in the first row come from the same initial data as in figure 5, but with the amplitude increased by a factor of 40 (A = 0.02, r 0 = 4, σ = 1 2 ). The left column shows the initial anisotropy function scaled by u −3 and plotted as a function of the inverse radial coordinate u, after adjusting the radial shift λ to fix the apparent horizon at u = 1. by a factor of 40 -and then examined the resulting evolution when the mean position of the initial Gaussian was progressively shifted deeper into the bulk. Each row of figure 7 displays results for one of these five cases. In each row, the left hand panel shows the initial anisotropy function, but plotted as a function of the inverse radial coordinate u, after adjusting the radial shift λ to fix the apparent horizon at u = 1, as discussed in section 2.6. The middle panels show the resulting pressure anisotropy as a function of time, and the rightmost panels display the nonlinearity (NL), again defined as the difference between the pressure anisotropy of the given initial pulse and twice the pressure anisotropy after halving the initial amplitude.
From the middle column of plots, one sees that the magnitude of the pressure anisotropy decreases significantly as the initial pulse is moved deeper into the bulk. Moreover, both the time it takes for the effect of the pulse to reach the boundary, and the width of the resulting peaks in the pressure anisotropy, grow with increasing depth of the initial pulse. This reflects the usual holographic mapping between bulk and boundary: phenomena deeper in the bulk correspond to lower energy or longer time scales in the boundary field theory.
Turning to the nonlinearity plots in the right hand column, one sees that the magnitude of the nonlinearity also decreases as the pulse moves deeper into the bulk. Dividing the peak nonlinearity by the peak pressure anisotropy gives a relative measure of nonlinearity. This is about 1 × 10 −4 for the top row, 3 × 10 −5 for the middle row, and 3 × 10 −6 for the bottom row. So in this comparison, as the initial pulse is pushed deeper into the bulk, the relative nonlinearity decreases systematically.
This comparison does not, however, imply that nonlinearities are never significant. The amplitude of the Gaussian pulse in the initial anisotropy function was kept fixed in figure 7, resulting in a decreasing size of the induced pressure anisotropy as the pulse moves deeper into the bulk. While the first two rows of the figure show pressure anisotropies which are large departures from equilibrium, ∆P/κε 1, the final rows with ∆P/κε 1 represent small departures from equilibrium. It is natural to ask what happens if one instead increases the amplitude as the pulse is moved into the bulk, so as to keep fixed the size of the induced pressure anisotropy.
Such a comparison is shown in figure 8, The upper row of the figure shows the time dependence of the pressure anisotropy and the nonlinearity for the same pulse which generated figures 5 and 6 (A = 5 × 10 −4 , r 0 = 4, σ = 1 2 ), while the lower row shows the pressure anisotropy and nonlinearity of a pulse with the same shape as the last row of figure 7, but with larger amplitude (A = 2.5, r 0 = 1 4 , σ = 1 2 ). The energy density remains fixed, ε = 3 4 L −4 . The peak pressure anisotropy is similar in the two cases, and corresponds to a large departure from equilibrium. For the latter case of a pulse deep in the bulk, large enough to induce a far from equilibrium pressure anisotropy, the nonlinearity is significant, much larger than the previous examples. However, even for this case, the size of the nonlinearity relative to the peak pressure anisotropy is only about 5%, NL/(∆P/κε) ≈ 0.05. 20 20 We have also examined the level of nonlinearity by comparing the pressure anisotropy resulting from a sum of two different Gaussians to the sum of anisotropies induced by the individual Gaussians. The results were comparable to those discussed above and do not warrant separate discussion. 2 ). Bottom row: corresponding plots of pressure anisotropy and nonlinearity for a pulse located deeper in the bulk (A = 2.5, r 0 = 1 4 , σ = 1 2 ) with amplitude adjusted to produce a similar peak pressure anisotropy. In both cases the energy density is ε = 3 4 L −4 . Substantial nonlinearity is present for this case, where an initial pulse deep in the bulk has sufficient amplitude to produce a large departure from equilibrium.
The data shown in figure 7 inspire several further questions. Looking down the middle column of the figure, one sees that the onset of the response (i.e., the time of the first peak in the pressure anisotropy) increases as the initial pulse moves deeper into the bulk but seems, perhaps, to be approaching a maximum value. Is this really true, or can one craft initial data for which the onset of the response is much greater? As shown on the left panels of the lower rows of the figure, when the average position r 0 of the Gaussian pulse is moved into the bulk, an increasingly large portion of the Gaussian ends up lying behind the apparent horizon. And when plotted as a function of our computational coordinate u = (r−λ) −1 (withr the λ = 0 frame radial coordinate), initial pulses with small values of r 0 are only moderately localized near the horizon -even though these pulses had constant widths when viewed as functions of r. As may be seen by comparing the left and middle columns of figure 7, it is the leading edge of the anisotropy pulse (the region of near-maximal slope) which produces the first large response in the boundary anisotropy. Looking at the last two rows of the figure, one may question whether we are doing an adequate job exploring the response from initial disturbances which are localized close to the horizon. Will initial pulses which are more strongly localized near the horizon show significantly greater nonlinearity? Figures 9 and 10 show results of an effort to explore these questions. Figure 9 shows the initial anisotropy function, along with the resulting pressure anisotropy and nonlinearity,  for a significantly narrower "deep pulse" (A = 1 10 , r 0 = 10 11 , σ = 1 20 ). And figure 10 shows a 3D plot of the time dependent anisotropy function, plus the induced pressure anisotropy, for an extremely narrow deep pulse (A = 1 10 , r 0 = 1, σ = 1 200 ). 21 The energy density in both cases remains fixed, ε = 3 4 L −4 . For the narrowest pulse, the amplitude A = 1 10 is near the upper limit which can be studied without destabilizing the horizon. In both figures 9 and 10, the peak pressure anisotropy ∆P/κε is large compared to unity, showing that these pulses are producing far from equilibrium initial states.
As seen in these figures, pulses which are more sharply localized very near the horizon do lead to a delayed onset in the resulting pressure anisotropy, occurring at vε 1/4 ≈ 1.75 for the case of figure 9 and vε 1/4 ≈ 3 for our narrowest pulse in figure 10. But, within the range of pulse widths we have studied, the onset of the pressure anisotropy response is only delayed by a factor of 2-3 compared to the case of figure 8. From the left panel of figure 10, showing the time dependence of the (rescaled) anisotropy function B s , one sees that the outward movement of the pulse toward the boundary resembles the behavior of outgoing null geodesics originating very close to an event horizon, which "hug" the horizon for JHEP07(2015)116 extended periods before eventually escaping. One may wonder if the onset of the response in the pressure anisotropy could be delayed indefinitely by going to narrower and narrower initial pulses localized at the apparent horizon. We do not have a firm analytic argument, but doubt that this is possible if one simultaneously demands that the amplitude of the response remain bounded away from zero. 22 The relative nonlinearity for the case of a narrow deep pulse shown in figure 9 is quite small, about half a percent. But for the extremely narrow pulse with near maximal amplitude of figure 10, the nonlinearity, relative to the pressure anisotropy, reaches the 10% level. Linearization of the dynamics about equilibrium must provide an accurate approximation to the full nonlinear dynamics when the deviation of the geometry from the equilibrium Schwarzschild black brane solution is sufficiently small, as will be true at sufficiently late times. For asymptotically anti-de Sitter geometries, where metric functions have the asymptotic forms (2.41), this will also be the case for initial data involving perturbations localized sufficiently close to the boundary. (See refs. [8,9,49] for related discussion.) It should be noted that reasonably good agreement between linearized dynamics and full nonlinear evolution was previously reported in ref. [8] and further explored in ref. [9]. In these works, the authors found agreement at a 20% level between the linearized and full dynamics for a variety of initial anisotropy profiles. Our results examining, more systematically, the dependence of the relative nonlinearity on the parameters of our initial Gaussian anisotropy function complement and extend this earlier work. Overall, despite prior knowledge of refs. [8,9], we are still surprised by the small, often extremely small, levels of nonlinearity which we find even at early times when the induced pressure anisotropy is maximal, for initial perturbations localized deep in the bulk and producing large departures from equilibrium.

Charged plasma
We now turn to the equilibration of charged plasmas (by which we mean SYM plasmas with a non-zero density of the global U(1) conserved charge -not a plasma in which electromagnetic forces are included in the dynamics and Coulomb repulsion plays a significant role). As noted earlier in section 2.7, the bulk geometry should equilibrate to a non-singular Reissner-Nordstrom black brane solution provided the charge and energy densities satisfy the extremality bound (2.58), ρ < ρ max = 4 3 ε 3/4 . However, for values of the charge density near ρ max , we find that sufficiently large initial metric perturbations can destroy the apparent horizon. We expect that such initial data are unphysical, not representing SYM initial states which could be produced by an operational procedure such as turning on time dependent external fields (which correspond to time dependent boundary data in the holographic description). In any case, we limit our attention to initial perturbations for which an apparent horizon is present at all times. We find that if one suitably decreases the amplitude of the initial departure from equilibrium while increasing the charge density, one can approach ρ max from below while maintaining the existence of an apparent horizon.  Figure 11 (left) compares the time dependence of the pressure anisotropy which results from initial data consisting of precisely the same Gaussian initial anisotropy function B(v 0 , r) (in the λ = 0 frame) and energy density as in figure 5, and a charge density ρ equal to 0, 20, 40, 60, or 80% of the extremal density ρ max . The immediately obvious qualitative result is that the five different curves are so close together than they cannot be visually distinguished! Varying the charge density (at fixed energy density and fixed initial anisotropy function) has stunningly little impact on the subsequent time evolution. This is quantified in the right panel of figure 11 which plots the difference in the pressure anisotropy ∆P/κε between the cases of ρ = 0.8 ρ max and ρ = 0. Comparing the scales on the right and left hand plots, one sees that for this initial anisotropy function the sensitivity to the charge density is less than one part in 10 4 .
In the plots of figure 11, we used the fourth root of the (rescaled) energy density, ε 1/4 , to set the scale for time. Since the energy density was held constant in the comparisons of figure 11, this was a simple and convenient choice. For the five cases shown in the figure, ε/(πT ) 4 = 0.75, 0.76 , 0.79, 0.86, and 1.03 for ρ = 0, 20, 40, 60 and 80% of ρ max , respectively. And, for comparison, the values of the equilibrium chemical potentials corresponding to these charge densities are given by µ/T = 0, 0.34, 0.73, 1.26, and 2.21, respectively.
If the initial anisotropy pulse begins deeper in the bulk, then the sensitivity to the charge density is larger. Figure 12 shows a comparison of pressure anisotropies for differing charge densities, now using the deep pulse initial anisotropy function whose radial profile has the shape shown in figure 9. 23 As seen on the left panel, the amplitude of the response increases significantly as the charge density varies from 0 to 80% of ρ max . However, the time course of the equilibration (e.g., the times of the first or second peaks in the response, 23 One might guess that the pulse used in the bottom row of figure 8 would exhibit greater sensitivity to charge since it had a larger nonlinearity than the deep pulse of figure 9. However, the latter pulse exhibits much greater sensitivity to charge. or the zero-crossing between these peaks) is only modestly affected, with changes of 3% percent or less.
The fact that the sensitivity to the charge density is greatest for pulses close to the horizon is to be expected. In the equilibrium geometry (2.57), one sees that as the radius r increases from the horizon, the influence of the charge density decreases rapidly relative to the other terms in the metric. Only near the horizon, and close to extremality, does the charge density produce an O(1) effect on the equilibrium geometry.
At the beginning of this work, we expected that one interesting outcome would be information on the change in equilibration time produced by varying the plasma charge density. By "equilibration time", we mean some rough but useful measure of when the departure from equilibrium is no longer substantial. To make this a bit more quantitative we adopt, somewhat arbitrarily, the criterion for all times t > t eq , as indicating that the system is near equilibrium at time t eq . Looking at the left panels of figures 11 and 12, it is clear that the effect of the charge density on any reasonable measure of equilibration time can be summarized easily: there is very little effect! Even in figure 12, where the sensitivity to charge density is the largest we have found with horizon-preserving initial data, 24 the time t peak of the initial response peak and the approximate equilibration time t eq both vary by less than 5%.
Since figures 11 and 12 plot time in units set by the energy density, the high degree of insensitivity of the relaxation time course to the charge density seen in these figures might lead one to think that the total energy density is playing a special role in setting the time scale of relaxation. But it should be borne in mind that these figures show comparisons in 24 Achieving good numerical accuracy becomes increasingly difficult as one pushes toward extremality, where the equilibrium solution bifurcates. Investigating very near extremal behavior more carefully is an interesting topic we leave to future work. which, by design, both the initial anisotropy function (in the λ = 0 frame) and the total energy density have been held fixed. Because the ratio of energy density to temperature (to the fourth power) varies significantly with increasing charge density, it is natural to ask whether the degree of (in)sensitivity of the relaxation dynamics to the charge density is substantially different if one holds fixed the equilibrium temperature instead of the energy density. Figure 13 shows such a comparison. Plotted are the pressure anisotropies resulting from the same initial anisotropy function used in figures 5 and 11, and charge densities of either 0 or 80% of extremality, but now with the energy density in either case suitably adjusted to fix the equilibrium temperature, πT = 1/L. One again sees a significant change in the amplitude of the response with increasing charge density, but now with the temperature held fixed, increasing charge density decreases the amplitude of the response. Nevertheless, with time now plotted in units of (πT ) −1 , one again sees negligible (≈ 0.01%) change in the time course of the equilibration dynamics as the charge density varies from zero and 80% of ρ max .
Performing the same constant temperature response comparison using the deep pulse initial anisotropy function (whose radial profile is shown in figure 9), we find a larger -but still quite small -variation in the time course, approximately 2%, as the charge density varies from zero to 80% of ρ max .
We have also examined the degree of nonlinearity in the above examples of equilibrating charged plasmas. The results for the relative size of the nonlinearity are quite similar to our earlier results for equilibrating uncharged plasmas. Because of this, we will refrain from presenting explicit nonlinearity plots for charged plasmas.
Finally, as noted earlier, at sufficiently late times the relaxation must be accurately described by a superposition of quasinormal modes (eigenfunctions of the linearized dynamics about the equilibrium solution). Extracting the leading quasinormal mode frequency by fitting the late time ( a decaying, oscillating exponential, as described in the previous section, yields the results shown in table 1. The second and third columns (with uncertainties) show our estimates for the real and imaginary parts of the leading QNM frequency in units of ε 1/4 , while the fourth column (without uncertainties) shows our estimates converted to units of πT . The rightmost column shows independent results of Janiszewski and Kaminski [50] obtained by analyzing the linearized small fluctuation equations about the RN black brane solution. The agreement is a satisfying confirmation of our numerical accuracy. Interestingly, the imaginary part of the lowest QNM frequency varies by over 15% between ρ = 0 and ρ = 0.8 ρ max . This is enormously larger than the part in 10 4 sensitivity seen in figure 11, and substantially bigger than the largest (≈ 5%) sensitivity we found in the evolution time course with deep pulses, figure 12. The implications of these very differing sensitivities will be discussed further in section 5.

Magnetized plasma
We now present results of an analogous investigation of equilibration in plasmas (with vanishing charge density) in a homogeneous magnetic field B. Our discussion will parallel, as much as possible, the previous treatment of charged plasmas. But the breaking of scale invariance by the magnetic field, discussed in section 2.1, produces some notable differences. One change, seen in section 2.4, is that the anisotropy function B(v, r) must now contain logarithmic terms in its near boundary behavior. We choose our initial anisotropy function to have the form (3.7) in which an adjustable Gaussian is added to the required leading logarithmic term. As in the previous discussion of charged plasmas, we will keep fixed the parameters of the Gaussian part of the initial anisotropy function B(v 0 , r) as we dial up the external magnetic field. We will also hold fixed the energy density ε L defined at a renormalization point µ = 1/L. As shown by the holographic relation (2.42), this is the same as holding fixed the asymptotic coefficient a 4 . In the following plots, axis labels involving energy density ε, without any explicit indication of scale, will denote the energy density evaluated at the curvature scale, ε(1/L) = ε L . Similarly, the pressure anisotropy ∆P should be understood as ∆P(1/L) unless otherwise indicated explicitly.
Instead of keeping ε L fixed as the magnetic field B is varied, one could choose to hold fixed the energy density ε B defined at the scale set by the magnetic field, µ = |B| 1/2 . Since the AdS curvature radius L is not a physical scale present in the dual QFT, fixing ε B instead of ε L is arguably more natural. However, for computational reasons it is easier to hold fixed ε L as the magnetic field is increased. The issue is that accurate numerical calculations become progressively more difficult the deeper one penetrates into the high-field/low-temperature regime, T 2 /|B| 1. (This is analogous to the difficulty of approaching extremality in the charged case, where the horizon temperature also vanishes.) By holding fixed ε L instead of ε B , we are able to perform scans in B which avoid dipping too deeply into the very low temperature region.
Another difference concerns the definition of pressure (or stress) anisotropy. With our choice (2.27) for fixing the scheme dependent ambiguity in the stress-energy tensor, the resulting holographic relation (2.42), when evaluated at µ = 1/L, puts the trace anomaly solely in the transverse components of the stress, T 11 and T 22 . So the pressure anisotropy (4.1), defined as the difference between transverse and longitudinal stress, when evaluated at µ = 1/L has a "kinematic" contribution of − 1 4 κB 2 plus a "dynamic" contribution of 3κ b 4 (v). In presenting results below, we will omit the uninteresting static kinematic contribution, and just plot the dynamic contribution (relative to the energy density), evaluated at the renormalization point µ = 1/L. 25 A further difference comes from the fact that equilibrium magnetic brane solutions are intrinsically anisotropic. The anisotropy function B does not vanish at late times, but rather settles down to the profile of the equilibrium magnetic brane solution discussed in section 2.7. This is illustrated in figure 14, which shows the (subtracted/rescaled) anisotropy function B s as a function of time v and inverse radial depth u. One sees similar features as in figure 5: the initial pulse propagates outward, reflects off the boundary, and largely disappears into the horizon. But in addition one also sees that the anisotropy function is approaching the non-trivial profile of a static magnetic brane solution.
Correspondingly, the pressure anisotropy in the dual field theory asymptotically approaches a non-zero constant. To examine equilibration, it is the difference between the pressure anisotropy and its asymptotic value which is of interest. As a measure of (near) equilibration, the condition (4.2) is naturally replaced by for all times t > t eq . The time dependence of the resulting (dynamical contribution to the) pressure anisotropy is shown in figure 15 for a series of magnetic fields, BL 2 = 0, 1, 2, 3, and 25 Examining the holographic relation (2.42), one sees that simply shifting the renormalization point to µ = e 1/4 /L would accomplish the same removal of this uninteresting kinematic contribution to the pressure anisotropy.  figure 15. 26 Small temporal variations are also evident after v ≈ 1.3. These are produced by the relaxation of the initial non-Gaussian background profile of the anisotropy function to the correct equilibrium form. These small variations at relatively late times would be absent if we had constructed our initial anisotropy function by adding a Gaussian perturbation to the equilibrium solution (instead of merely adding a Gaussian to the leading asymptotic term). Given our choice of initial data, there are two distinct time scales in the equilibration shown in figures 14 and 15, the first associated with the boundary reflection and subsequent infall of the Gaussian pulse, and the second with the time it takes for the background anisotropy profile to reach its equilibrium form. It is the latter which is responsible for the late time variations; fortunately, there is relatively little ambiguity in separating the two contributions to the dynamics.
Once again, a notable feature of in the comparison of figure 15 is the similarity in the time dependence of the pressure anisotropy during the pulse-driven period of significant departure from equilibrium (0.2 vε 1/4 L 1.0), when the initial perturbation is reflecting from the boundary. Any reasonably defined measure of equilibration time t eq clearly does not vary much with magnetic field, and neither does t peak for these relatively near boundary pulses. This insensitivity result relies, of course, on the constancy of the initial Gaussian perturbation in the anisotropy function, and also on our choice to hold fixed the energy density at the scale 1/L. Figure 16 compares the response, for different values of magnetic field, when one holds fixed the equilibrium temperature T (as well as the initial Gaussian perturbation), instead of fixing the energy density ε L . With time now plotted in units of (πT ) −1 , one also sees remarkable similarity in the time dependence of the response, with the times of the first, second, or third peaks in the pressure anisotropy varying by less than 0.3% as B/T 2 varies from 0 to 26.7.
Sensitivity to the magnetic field is significantly larger when the initial pulse is placed very close to the horizon. This is shown in figure 17, which plots the time dependence of the dynamical pressure anisotropy for magnetic field values BL 2 = 0, 1, 1.5, 2, and 2.5, using the same "deep pulse" Gaussian parameters (A = 1 10 , r 0 = 10 11 , σ = 1 20 ), and fixed energy density ε L = 3 4 L −4 , which generated figure 9. For this series of solutions we have 26 The tiny positive late time pressure anisotropy barely visible in the BL 2 = 1 curve is a consequence of our removal of the static kinematic part of the anisotropy; the final value of the total pressure anisotropy, at the scale 1/L, monotonically decreases with increasing BL 2 in this series of solutions.   to increasingly large and negative final values for the anisotropy. For the lowest curve with BL 2 = 2.5, the contribution of the Gaussian pulse is completely swamped by the contribution from the relaxation of the background profile of the anisotropy function to the correct equilibrium form. Excluding this curve, the time t peak of the first peak, as well the rough equilibration time t eq , characterizing the portion of the evolution arising from the Gaussian pulse, vary at most by 20% as the magnetic field ranges from 0 to 2/L 2 . This is the largest difference in the relaxation time course we have seen in our exploration of magnetized plasmas with Gaussian initial perturbations.
A constant temperature comparison (not plotted), analogous to figure 16 but using the same deep pulse as in figures 9 and 17, shows variations in the time course of up to 15% as |B|/T 2 ranges from 0 to 22 -beyond which the response of the pulse cannot be clearly distinguished from the background evolution.
To conclude this section, we report in table 2 equilibrium properties, plus our estimates for the lowest quasinormal mode frequency, for values of magnetic field which extend from small fields well into the strong field regime. The ratio B/T 2 ranges from zero to just over 30. The equilibrium energy density, pressure, and pressure anisotropy are given in units of T 4 , and have been converted from the µ = 1/L renormalization point used in our calculations (and above presentation) to the intrinsic scale µ = |B| 1/2 . Results for the lowest quasinormal mode frequency are given both in units of ε

Discussion
The above results show that to a good (often extremely good) level of accuracy: 1. the pressure anisotropy response is a linear functional of the initial anisotropy pulse profile; 2. the time course of the response, measured in units set by the energy density, 27 is insensitive to the charge density or background magnetic field when the pulse profile and the energy density are held fixed; 3. the time course of the response, measured in units set by the equilibrium temperature, is insensitive to the charge density or background magnetic field when the pulse profile and equilibrium temperature are held fixed.
How can one synthesize these observations? Consider some feature in the time course of the response, such as the time of the first (or second) peak in the pressure anisotropy, or the approximate equilibration time discussed above. For simplicity, we will focus on the time t peak of the first pressure anisotropy peak. This time must be some function of the equilibrium state parameters (energy density plus charge density or magnetic field), as well as the Gaussian pulse parameters (depth, width, and amplitude). Consider first the charged case. The response time t peak is a function of the energy density ε, the fraction x of the extremal charge density, and the Gaussian pulse parameters r 0 , σ, and A. But since the temperature decreases monotonically with increasing charge density (for fixed energy density) one may equally well regard the equilibrium state as labeled by the energy density ε and temperature T , and write t peak /L = f (εL 4 , T L, r 0 /L, σ/L, A) , (5.1) for some function f of the indicated arguments. We have written all quantities in dimensionless form using, in effect, our computational units. Our results on the degree of nonlinearity imply that the function f is nearly independent of the last argument, the Gaussian amplitude A. Only for our narrowest pulse, located right at the horizon, did the relative nonlinearity reach 10%. Away from this corner of parameter space, the nonlinearity was substantially smaller, rapidly falling to much less than a percent as the initial pulse becomes less localized at the horizon. So, to a good approximation, one can regard the function f as being independent of A. For narrow pulses, σ r 0 , the dependence of the response time t peak on the pulse width is negligible; there is a smooth σ → 0 limit. To the extent that linearity is a good approximation, one may regard the response from wider pulses as superpositions of the response from narrower pulses (suitably arranged so that their sum reconstructs the wider pulse). As noted in the discussion of figure 7, the first peak in the response is clearly associated with the propagation of the leading edge of the anisotropy pulse -the region of near-maximal slope on the side of the pulse closest to the boundary. Therefore, for pulses of non-negligible width σ centered at some depth r 0 , one should expect that the time of 27 Specifically, εL in the magnetic case. the first peak in the response will be most similar to the corresponding response time for a narrower pulse located not at the depth r 0 , but rather at a depth of r 0 + nσ for some positive O(1) multiplier n. At the same level of accuracy determined by the degree of nonlinearity, one should be able to merge the dependence of the response time t peak on the depth r 0 and width σ of the initial pulse into a single effective depth r eff given by r 0 + nσ. (The accuracy of this simplification will be discussed below.) Hence, the above functional dependence for t peak can be replaced by a simpler form, t peak /L ≈ g(εL 4 , T L, r eff /L) , (5.2) for some function g. Now, the results of section 4.2 (including figures 11 and 12 and associated discussion) show that the time course of the pressure anisotropy response has remarkably little dependence on the charge density when comparisons are made holding fixed the initial pulse parameters and the energy density. Since varying the charge density at fixed energy density is, as noted above, equivalent to varying the equilibrium temperature, this implies that the function g describing the response time (5.2) is nearly independent of the second argument, T L. At the same time, the comparisons at fixed temperature also discussed in section 4.2 ( figure 13 and associated discussion) show that the time course of the pressure anisotropy response also has remarkably little dependence on the charge density when the initial pulse parameters and the equilibrium temperature are held fixed. This implies that the function g describing the response time (5.2) is nearly independent of the first argument, εL 4 . Hence, at a level of accuracy determined by the minimal level of nonlinearity, and the minimal dependence on charge density in comparisons at both constant energy density and constant temperature, the response time t peak must be a function of only the effective depth of the initial pulse, t peak /L ≈ h(r eff /L) , for some function h. Finally, this function must be consistent with the scaling relations discussed in section 2.5, which imply that L 2 /r scales in the same fashion as a distance (or time) in the boundary theory. Consequently, the dependence of the response time on the effective depth must have the form for some dimensionless constant C. A similar line of reasoning is applicable to the magnetic case. Since we used ε L ≡ ε(L) and BL 2 as parameters labeling the equilibrium magnetic brane geometry in our comparisons of magnetic plasma response, it is convenient to view the response time t peak as depending on these parameters plus the Gaussian pulse parameters, t peak /L = f (ε L L 4 , BL 2 , r 0 /L, σ/L, A) , (5.5) for some function f . Once again, the observed near-linearity of the response allows us to simplify this to the form t peak /L ≈ g(ε L L 4 , BL 2 , r eff /L) ,  Figure 18. Time of the first peak in the pressure anisotropy response as a function of the inverse effective depth of the pulse, u eff ≡ 1/r eff . The left panel shows results for narrow pulses with σ = 1/20, while the right panel shows results from wide pulses with σ = 1/2. Blue circles represent data points from neutral (uncharged, unmagnetized) plasmas, maroon squares represent data points from charged plasmas at 80% of extremality, and gold diamonds represent data points from magnetized plasmas with BL 2 = 1.5. The green straight line shows the prediction of our simple model (5.7) with n = 2.5 and C = 2.04.
for some function g. The near-independence of the time course of the response on the magnetic field B, for fixed ε L and a fixed initial pulse, implies that the function g is nearly independent of its second argument. Because ε L does not transform homogeneously under the scaling relations (2.43)-(2.45) (due to the use of the curvature scale L instead of a physical scale in the dual QFT for setting the renormalization point) consistency with the scaling relations requires that the function g be independent of its first argument and depend inversely on the third. So, just as for the charged case, for some dimensionless constant C. Figure 18 compares this simple model with a sample of our data for neutral, charged, and magnetized plasmas. The left panel shows data for relatively narrow pulses of width σ = 1/20, while the right panel shows data from rather wide pulses with width σ = 1/2. The abscissa for both panels is the inverse effective depth u eff ≡ 1/r eff . (Data points at the largest values of u eff shown in these plots come from pulses centered very near the horizon.) In both panels, rather good agreement with the simple model is found when the multiplier n in the effective depth r eff ≡ r 0 + nσ is chosen to be 2.5, and the coefficient C ≈ 2.04 . (5.8) For both sets of data, the model is least accurate for pulses located very close to the horizon. Accuracy for the case of magnetized plasma is a bit worse than for charged or neutral plasma. But for all cases, even in the near-horizon regime, this simple model works at about the 20% level or better. As pulses move away from the horizon, the accuracy rapidly improves.

JHEP07(2015)116 6 Conclusion
In weakly coupled plasmas, adding a conserved charge density to the system (e.g., flavor charge in a QCD plasma) significantly changes screening lengths. When the associated chemical potential is of order πT , the relative change in the Debye screening length is O(1) [51]. Such changes in screening lengths significantly affect transport coefficients like the viscosity [52]. Lattice studies [53] of the effect of a baryon chemical potential in the deconfined phase, when the system is not asymptotically weakly coupled, claim to find measurable sensitivity in the Debye screening length, comparable to perturbative estimates, but with quite large error bars. Lattice QCD studies of magnetoresponse [29][30][31][32] also find substantial changes in thermodynamics when the magnetic field energy density becomes large compared to T 4 . Consequently, when this work on strongly coupled N = 4 SYM plasma was initiated, we expected to find significant changes in equilibration dynamics when a conserved charge density is added to the plasma, or when the system is placed in a background magnetic field. The most notable result we have found is that this expectation was wrong. At least within the range of charge densities we studied, up to 80% of extremality, the equilibration dynamics is remarkably insensitive to the presence of a conserved charge density. Additionally, magnetic fields which are well into the strong field region, B/T 2 1, induce almost no change in the equilibration time course.
Efforts to use results of holographic calculations in strongly coupled N = 4 SYM as the basis for predictions about real heavy ion collisions [21] are inevitably hampered by our limited understanding of the effect on the relevant dynamics of changing the theory from real QCD to a supersymmetric Yang-Mills model theory. From this perspective, the insensitivity of the equilibration dynamics to the charge density is reassuring, as this provides an example where changes in the plasma constituents have very little impact on the overall dynamics.
A further notable feature in our results is the remarkably small degree of nonlinearity in the dynamics governing the pressure anisotropy. Despite the fact that one is solving the highly nonlinear Einstein equations, the dependence of the induced pressure anisotropy on the initial anisotropy function is surprisingly close to linear. We find deviations from linearity of at most ≈ 30%, and this only for initial disturbances which are crafted to reside very deep in the bulk. (This is consistent with earlier work in ref. [9].) For disturbances localized even modestly above the horizon, the degree of nonlinearity quickly drops to subpercent levels. This near-linearity holds even when the system is far from equilibrium, with pressure anisotropies which are large compared to the energy density.
In the charged case, the sensitivity of the lowest quasinormal mode frequency to the charge density is significantly larger than the sensitivity we find in, and shortly after, the far-from-equilibrium period of the equilibration dynamics. The lowest quasinormal mode dominates the equilibration process at sufficiently late times when all higher modes have decayed away and become negligible relative to the lowest mode.
The near-linearity which we find in the dynamics implies that the deviation from equilibrium, even during the far-from-equilibrium portion of the process, can be represented -43 -

JHEP07(2015)116
quite accurately as a sum of quasinormal modes obeying linearized small fluctuation equations. The minimal sensitivity of the equilibration dynamics to the charge density, substantially less than the sensitivity of the lowest quasinormal mode frequency, suggests that during most of the equilibration process many quasinormal modes above the lowest one are contributing, with the sensitivity to the charge density quickly falling with increasing mode number. This is something which could be tested directly in a linearized analysis but, to our knowledge, has not yet been done.
In the magnetic case, we find only modest (few percent) changes induced by the magnetic field in the time course to equilibration, and in the lowest quasinormal mode frequency characterizing very late time behavior. It would be nice to have independent calculations of quasinormal mode frequencies for magnetic branes, including studies of the field dependence of higher mode frequencies.
Although we no longer have reason to expect large effects, it would be natural to generalize the study of homogeneous equilibration to the case of plasmas with both non-zero charge density and a background magnetic field. Extensions to more complicated inhomogeneous settings, such as colliding shock waves, are also of interest. It would be desirable to gain a clearer understanding of the connection between our choices of gravitational initial data on null slices and boundary observables such as multi-point stress-energy correlators. Can more operational procedures, such as time-dependent background fields [7], produce far-from-equilibrium states which resemble those produced by our "deep pulse" initial data? We hope future work can shed light on some of these topics.