Nonlinear conductivity and the ringdown of currents in metallic holography

We study the electric and heat current response resulting from an electric field quench in a holographic model of momentum relaxation at nonzero charge density. After turning the electric field off, currents return to equilibrium as governed by the vector quasi-normal modes of the dual black brane, whose spectrum depends qualitatively on a parameter controlling the strength of inhomogeneity. We explore the dynamical phase diagram as a function of this parameter, showing that signatures of incoherent transport become identifiable as an oscillatory ringdown of the heat current. We also study nonlinear conductivity by holding the electric field constant. For small electric fields a balance is reached between the driving electric field and the momentum sink — a steady state described by DC linear response. For large electric fields Joule heating becomes important and the black branes exhibit significant time dependence. In a regime where the rate of temperature increase is small, the nonlinear electrical conductivity is well approximated by the DC linear response calculation at an appropriate effective temperature.

In this paper we turn to an investigation of the nonlinear response of such holographic models, following the quench of an applied electric field. From the gravitational point of view, this amounts to an investigation of the nonlinear dynamics associated to black branes in AdS with inhomogeneous Dirichlet boundary conditions. One should expect qualitatively different dynamical behaviour for bulk evolutions with such boundary conditions, as compared with those with homogeneous boundary conditions. At the very least the boundary CFT dynamics are distinct, with the former conserving momentum either JHEP10(2016)008 only approximately or not at all. 1,2 Thus at late times we do not expect to find hydrodynamic behaviour, but rather only exponentially fast decay described by the quasi-normal modes (QNMs) associated to momentum relaxation. The QNM spectrum itself is known to change qualitatively with the strength of inhomogeneity, resulting in transitions between coherent and incoherent metals [28][29][30], offering up further interesting features for the bulk dynamics.
From the field theory point of view, an electric field is a natural tool to quench an inhomogeneous system with finite charge density. 3 The non-conservation of momentum allows the possibility of approaching a finite steady state, at least in the linear response regime. In particular, as a special case this includes the relaxation of currents to zero in the laboratory frame after the electric field is turned off. Here we can study the imprint of the QNM spectrum on the current relaxation, in particular how it is affected by the aforementioned coherent to incoherent shifts in the QNM spectrum. 4 For constant electric fields we compute nonlinear electric and thermoelectric conductivity. Here the ideal scenario would involve a steady state solution, perhaps utilising an external heat bath. This is precisely the kind of scenario that is argued to arise in some probe brane constructions [34][35][36][37][38][39][40][41]. In the present backreacted context a steady state will not be reached since energy is continually added and the system exhibits Joule heating without bound. It would be interesting to consider models which realise a coupling to an external heat bath incorporating backreaction -we consider some of results presented here as a precursor to this problem, which we shall leave to future work.
Remarkably however, even in the absence of a steady state, nonlinear electrical conductivity has been computed exactly in the absence of charge [42], explored using a Vaidya spacetime. Similarly in this paper we will be able to make clean statements about the nonlinear electric and thermoelectric conductivities at finite charge density in the absence of a steady state, at least in a limit where the rate of heating is low.
To motivate our choice of bulk model, we note that at finite charge density, the electric field will drive both electric and heat currents. It is therefore desirable that the model used incorporates some form of momentum relaxation. The minimal ingredients required are provided by an Einstein-Maxwell-axion model [19]. Each axion field, φ I , is dual to an operator O I sourced by the value of φ I at the AdS boundary, φ (0) I . To break translational invariance the sources are taken to be linear in the boundary spatial coordinates φ (0) If two axions are used, the equilibrium state can be made isotropic, although driving the system with an electric field will break isotropy in general. The axion model brings a clear computational benefit since the bulk geometry remains independent of x I . 5 1 Previous work on the dynamics associated to momentum relaxation was discussed in a complementary class of holographic experiments in [7].
2 See [26,27] for approaches to a hydrodynamic-like description when momentum relaxation is weak, in the context of axion models. 3 For studies of electric field quenching in a probe brane context see for example [31][32][33]. 4 We thank Richard Davison for encouraging pursuit of these features at nonzero charge density. 5 The only x I dependence appears in φI and the gauge field, A, where it is known in advance.

JHEP10(2016)008
Some basic aspects of the role of the electric field in the axion model can be understood in reference to the Ward identities, where T µν is the QFT stress tensor, J µ the U(1) current and F (0) = dA (0) is a classical, boundary field strength. Throughout this paper the sources are taken to be, where I = 1, 2 labels the two spatial boundary directions, x = x 1 , y = x 2 . With the exception of φ (0) I and A (0) all quantities here are independent of x I . Inserting these expressions into (1.1) we find, where = T tt is the energy density, ρ = J t is the charge density, J = J x is the electrical current, J E = T tx is the energy current. 6 The first equation (1.4) shows that energy is injected by the electric field, which results in Joule heating. The second, (1.5), accounts for the momentum, and because of the axion sink term there is the possibility of a steady state where the terms on the hand side cancel, balancing momentum relaxation with the driving electric field. Finally we note that (1.2) gives ∂ t ρ = 0. Since we shall be subject to the nonlinearities of the full evolution of the Einstein-Maxwell-axion system, we resort to a numerical integration of the bulk equations of motion. Details of the numerics are given in appendices A and B. We will study two different cases for the profiles of E(t): • Top hat E(t) -current relaxation. In the first case we smoothly turn on E(t) hold it at some constant value, and then turn it off again. This experiment will confirm crucial aspects of the model, first and foremost, its stability. We confirm that the system returns to a member of the family of equilibrium solutions at a rate consistent with its longest lived QNMs. In particular we will examine the QNM spectrum for different values of the momentum relaxation parameter k, and show its imprint on current relaxation. We will see a qualitative change in the relaxation of the heat current following a quench to an incoherent regime. Additionally, the return to equilibrium will allow us to get a complete picture of the gravitational situation, allowing the identification of the black hole event horizon. • Step E(t) -nonlinear conductivity. For the second case, we smoothly turn on E(t) and hold it at a constant value E f . At sufficiently late times, transient behaviour associated to the quenched electric field dies down and we enter a regime JHEP10(2016)008 where the system is effectively responding only to a constant E = E f . When E is constant we choose to characterise the response of the system using electrical (σ) and thermoelectric (ᾱ) DC conductivities, defined as follows, where Q is the heat current in the absence of thermal gradients [8]. Here we have restricted to the xx-entries in each conductivity matrix, with the yx-entries vanishing.
The thermal conductivity will not play a role since there is no temperature gradient. Note that σ andᾱ will be time dependent in general due to the effect of heating, and they may also depend explicitly on E f . Such examples of perpetually driven CFTs are less well studied in favour of a system which eventually returns to equilibrium. An interesting set of examples are provided in [43].
In the context of the second class of experiments we can gain analytic control over σ,ᾱ in two limits of parameter space. For sufficiently small E f , linear response can be used, finding [8,19] In this regime O 1 = −E f ρ/k and the right hand sides of (1.4) and (1.5) vanish at linear order in E f , corresponding to a steady state. For ρ = 0 one can go beyond linear response analytically for any E f , finding σ = 1,ᾱ = 0. (1.9) In fact, these expressions apply for any choice E(t), showing that J(t) responds instantaneously to E(t). These conductivities are extracted from a Vaidya-like spacetime which can be constructed with axions [44]. This generalises the conductivity results of [36,42] to include a momentum relaxation parameter. 7 In the absence of net charge, the added feature of momentum relaxation may be somewhat redundant, nevertheless these solutions serve as a useful reference point. Indeed we will show that an instantaneous electric current is the first response of the system in the charged examples that follow. Away from these two limits the responses of J and Q are subject to the nonlinearities of the system. Remarkably, we find that the expression for σ in (1.8) provides an excellent accounting of the evolution of J far from linear response, provided some local notion of temperature during the evolution. This applies even when Joule heating introduces significant time dependence of the black brane. On the other handᾱ does not appear amenable to the same treatment, exhibiting nonlinear dependence on E f . The paper is organised as follows. In section 2, we present the equilibrium Einstein-Maxwell-axion black branes. In section 3 we review Vaidya-like solutions to the Einsteinaxion system, and for the case of an electric field compute the nonlinear σ andᾱ at ρ = 0, as given in (1.9). We then turn to the numerical solutions for a top hat electric field profile JHEP10(2016)008 in section 4 including a discussion of QNMs and the transition to an incoherent regime. We examine nonlinear conductivity in section 5. We conclude in section 6. In appendix A we detail the numerical setup used, with convergence tests presented in appendix B.

Gravitational model and equilibrium black branes
Let us briefly recap the model and present the black brane solutions of [19,44], presented in ingoing Eddington-Finkelstein coordinates and in a generalised coordinate frame on the boundary. These give the equilibrium, finite temperature, charged metallic state which we will later force out of equilibrium using E(t). We employ the action specialised to D = 4, where Λ = −3 −2 and F = dA. The units employed throughout set 2κ 2 4 = 1 and we set = 1. 8 This theory admits a family of equilibrium black branes which can be written in the following form, where x µ = (v, x, y) with constant 3-vectors u and n I , I = 1, 2. The functions appearing are given as follows, which is a solution to the equations of motion provided we satisfy the following orthogonality condition, which is equivalent to the 6 relations, u 2 = −1, η µν n µ I n ν J = δ IJ , n I µ u µ = 0. The stress tensor and current are given by, where = 2m. For convenience we work in a fixed laboratory reference frame, setting u = ∂ v and n I = ∂ x I , which satisfies all of the above conditions for a solution.
The thermodynamical quantities can be written most succinctly in a form parameterised using the event horizon coordinate location r 0 (where f = 0),

Vaidya-like solutions with axions
It is also possible to construct Vaidya-like spacetimes in the presence of the axion sources [44].
Here we review these solutions and compute their nonlinear electric and thermoelectric conductivities in the vein of [42]. Specifically, in the D = 4 Einstein-axion system, we can add an additional stress tensor source,T ab , leading to the equations of motion For a null stress tensor of the form,T vv = r 2 S(v) with other entries zero, the following electrically neutral solution is obtained, with the sourced energy conservation equation, As pointed out in [44] an electric field can be used to generate a stress tensor of this type. Specifically for a Maxwell field in the bulk, F = E(v)dx ∧ dv gives rise to a stress tensor of the required form, with S(v) = E(v) 2 . For Einstein-Maxwell in the absence of axions such solutions were given in [36,42]. For this solution the electric and heat currents are J = E and Q = 0 -specialising to a constant electric field E(v) = E f we see that the nonlinear electric and thermoelectric conductivities are given by (1.9). Note that this does not straightforwardly extend to D = 4; the conductivity is no longer dimensionless and is influenced by heating, as was explored in the Einstein-Maxwell context in [42]. Similarly we have not been able to extend these results to ρ = 0; momentum is nonzero and the solutions do not fit into the form above. The remainder of this paper is focussed on a numerical study of the ρ = 0 case. The analytical results presented here have been used as a test case for the numerical evolution described in the remainder of this paper. Finally we note that the time dependent solution in this section can be used to quench the system from an initial equilibrium state with < 0 to a final equilibrium state with exactly = 0, which occurs at k = √ 2/r 0 . Despite the remarkable behaviour of the perturbations of this state [28,30] there appears to be nothing remarkable about a quench so designed.

JHEP10(2016)008 4 Current relaxation
In this section we drive the system out of equilibrium and then let the energy and electrical currents return to their equilibrium value -zero, in the laboratory frame defined by the axion sources. We take the electric field profile where Θ(t) ≡ 1 2 tanh t w + 1 where w is taken to be some short timescale compared to the intrinsic dynamical response time of the system.
The return to equilibrium of the currents is governed by the QNM spectrum in the vector sector of perturbations, which changes qualitatively with varying k/ √ ρ and T / √ ρ.

Vector QNMs at ρ = 0
One universal statement is that for sufficiently small values of k 2 there is a long lived QNM with a purely imaginary frequency, well separated from the other QNMs. In this regime the metal is coherent, with a Drude-peak appearing in the conductivity at low frequency. One can then associate a momentum relaxation rate Γ rel with the decay rate of this QNM, which has been computed using a matching calculation [26], For larger values of k 2 one must examine the QNM spectrum in detail.
The vector QNMs are given by ingoing, normalisable solutions to a pair of decoupled equations for gauge invariant master fields, labelled by (±), for details see [19,29]. At ρ = 0 the QNM spectrum indicates both coherent and incoherent behaviour is possible [28][29][30]. In [28] the transition between the two behaviours was demonstrated to coincide roughly with the location of a pole collision, producing a damped-oscillatory pole as a dominant contribution to QQ . We find that this transition persists at ρ = 0 in the (−)-sector of QNMs, whilst the (+)-sector is dominated by a longer lived purely decaying mode. Since in the general case two-point functions of J and Q receive contributions from both (±)-sectors of QNMs, the relaxation of J and Q are dominated by the longer-lived (+)-sector mode.
However, the currents can be decoupled -see [29] for details. To each of the (±) sectors one can assign currents, J ± ≡ a ± (J E + γ ± kJ), which satisfy J ± J ∓ = 0, where a ± is an overall normalisation and This relation allows one to relate the black hole QNM spectrum in each sector to the thermoelectric response. However, as can be seen in the structure of (4.3), the precise combination of physical currents which achieves this decoupling depends on the details of this particular holographic model. Universal decoupling is achieved at ρ = 0 where simply J + = J and J − = J E = Q. But, there is a second limit, k → ∞ where remarkably J + = J and J − = Q [29]. In other words, if we take large enough k, the amplitude of the (−)-sector QNMs become parametrically enhanced over the (+)-sector QNMs in the phenomenology of Q, whilst the (+)-sector QNMs remain dominant for J.
JHEP10(2016)008 We have not performed an exhaustive analysis of this vector QNM spectrum. In figure 1 we show the longest lived QNMs in each sector (±) as a function of k/ √ ρ at T / √ ρ = 7/10. We reiterate that even though the (−)-sector modes are shorter lived, they can still contribute significantly to the response of Q provided k is large enough. From figure 1 we see that the (+)-sector appears to show no significant features, but the (−)sector becomes dominated by oscillatory-decaying modes at large enough k/ √ ρ, just as in the case ρ = 0. Fortunately, this behaviour is also the regime where the (−)-sector can contribute to the response of Q. Thus the relaxation of Q will undergo a qualitative change as k/ √ ρ is dialled.
At this temperature one can see this oscillatory mode results from a pole collision, which closely resembles the neutral case. As one takes higher T / √ ρ the case ρ = 0 is approached an the red and black curves in figure 1 connect in the vicinity of their closest approach, giving the behaviour observed in [28,30] where the Drude pole continuously connects to the pole collision. At lower T / √ ρ the oscillatory mode remains but the dominant Re(ω) = 0 piece of the (−)-sector disappears, along with the pole collision.
In the next subsection we examine the far from equilibrium dynamics of J and Q for quenches which return to equilibrium in the two qualitatively different regimes of the QNM spectrum shown in figure 1.

Quenches
In this section we study electric field quenches which eventually settle down to equilibrium states whose QNMs are given by the spectrum shown in figure 1. As argued, by dialling k/ √ ρ we expect to see a qualitative change in the relaxation of Q, due to the pole collision and off-axis mode at large k/ √ ρ. The parameters considered here are labelled by the dashed grey lines in figure 1.
We first examine the general features of a quench at k/ √ ρ = 2. For this example we choose E c /ρ = 3 with √ ρt * 5, with initial temperature T i / √ ρ = 3/10. In figure 2 we show the behaviour of J, J E and the black brane horizon. During the quench in which E is turned on, and shortly after, J(t) E(t) as in the neutral theory. After E returns to zero the system returns to an equilibrium black brane with T / √ ρ 0.690.
The relevant QNMs for this example are purely decaying, with the relaxation of J and Q governed at late times by the (+)-mode, which for the specific T reached here is has JHEP10(2016)008  73i. For Q we show a fit of a linear combination of these two QNMs. For J we fit only the (+)-sector mode. As we previously argued, it is expected that the (−)-sector make a dominant contribution to Q despite being much shorter lived, since in the limit k → ∞, Q is sensitive only to this sector. Indeed in figure 4 we see the oscillatorydecaying behaviour of Q (or ringdown, appropriating the dual black hole terminology), but not of J, as expected. Eventually the contribution of the longer lived (+)-sector mode appears in Q due to 1/k effects; we have verified that increasing k extends the time over which the (−)-sector mode governs Q.

Nonlinear conductivity
In order to investigate nonlinear thermoelectric response we quench the system away from equilibrium with a step profile, where Θ(t) is defined in section 4.

Defining an effective temperature
In the numerical sections that follow it will be useful to have a notion of effective temperature, T E (t), during time evolution. We define T E (t) as the temperature of the equilibrium state to which the system eventually settles down if we stop driving the system at time t. This definition has the benefit of carrying a clear physical interpretation at any time, even far from equilibrium, and it is precisely the temperature when the system is in a stable  equilibrium. It also extends beyond the use of just the local energy density, (t), accounting for the possibility of other time dependent charges. In the presence of an instability, determining T E (t) clearly can become an involved task, and whether or not it would be a useful quantity in such cases remains to be seen. In all the examples given here, no instabilities are seen and the late time equilibrium solutions are known, given in section 2. Moreover at fixed k and ρ (which are both constants for the evolution) equilibrium can be labelled by the energy density, , thus we can easily compute T E (t) by computing the temperature of the equilibrium state with energy (t).

Linear regime
Here we take the example of E f /ρ = 10 −3 showing J(t)/E f and Q(t)/E f /T in figure 5 together with the expected values assuming DC linear response, (1.8). As anticipated the immediate response is for J to track E during the quench, corresponding to the kink from J/E f = 0 to J/E f = 1 around t = 0. Following this, momentum grows until it is balanced by the momentum sink, and the charged, linear response steady state is achieved. In figure 6 we verify the approach to the steady state is governed by the longest lived vector QNM, by plotting the scalar vev, which in the steady state takes the value O 1 = −E f ρ/k.
In this example the electric field is small but finite, and so there is a small Joule heating effect, with a rate of temperature increase, Γ Joule . Near equilibrium the Ward Identity (1.4) can be used to give, where c ρ is the specific heat, to be evaluated for the equilibrium black brane. Deviations will occur for strong electric fields, or after significant heating. We show T E (t) in the right panel of figure 6, together with the timescale Γ Joule .

Examples with significant Joule heating
For a steady state we may reasonably expect J to depend on T / √ ρ -as indeed it does in linear response -but we may also expect J to depend nonlinearly on E f . In the absence of a steady state there is the additional complication of time dependence. However, if the rate of heating, Γ Joule , is sufficiently low, we may approximate the time dependence of J through its temperature dependence by promoting T → T E (t), Higher order corrections in the heating rate Γ Joule could also be considered systematically in a derivative expansion, but for now we focus on this leading effect. We shall demonstrate that the T E dependence is remarkably simple -good agreement is achieved by taking the DC linear response result for σ, which may be written in terms JHEP10(2016)008 of the temperature at equilibrium, T , and promoting T → T E (t). In other words, σ does not appear to depend on E f explicitly. Let us refer to σ computed in this way as σ lin . A practical short-cut for this is to first eliminate µ from σ (1.8) using the thermodynamic relations (2.9) and (2.10) to obtain, Since ρ is conserved and k is fixed, it is only the energy density, , which we update during the evolution. Note that the linear response result forᾱ (1.8) is constant since ρ is conserved. In figure 7 we show J/E f as measured during the evolution as a function of T E at E f /ρ = 1. After some transient period σ approaches the attractor governed by the linear response approximation, independently of the initial temperature. We have verified that in the regime where the linear response result applies, subleading corrections due to higher derivatives of T E are small. 9 We have also verified similar behaviour for the momentum relaxation values, k/ √ ρ = 1/2 and k/ √ ρ = 8. We note that for large k/ The plateau coincides with a period of rapid temperature increase; one possible explanation for it is then that the growth in temperature is much faster than the response time due to the presence of charge and the system appears neutral, hence the instantaneous J = E behaviour just as in section 3.
Finally we illustrate the dependence on E f in figure 8, up to E f /ρ = 10. The linear response curve, σ lin is eventually reached for each case. Γ Joule is higher for larger E f , therefore the system is at a higher T E by the time any transients have settled down and σ lin is reached. For this reason, whether or not the σ lin approximation applies for arbitrarily large E f /T 2 E becomes difficult to assess via a quench. Here we have at least good agreement with σ lin in cases where E f /T 2 E ∼ 1. We note that the linear response result forᾱ does not give good agreement over the same timescales indicatingᾱ has an explicit dependence on E f . This is demonstrated in figure 9, where at low enough T i / √ ρ at fixed E f /ρ = 1 the effective temperature dependence ofᾱ deviates from the linear response result,ᾱ lin .

Final comments
We have presented the time evolution of a holographic metal under the influence of an applied electric field at finite temperature. At ρ = 0 the system responds instantaneously 9 Note similar agreement for the σ approximation can be observed in example of the top hat electric field profile presented in section 4, corresponding to the dashed blue curve in figure 2.
JHEP10(2016)008 to the applied electric field, encoded by a Vaidya-like geometry. There is Joule heating but no momentum and no heat current. Such solutions may be useful reference points for the construction of steady states, or as we have investigated here the addition of charge, where they govern the initial response of the system. At ρ = 0 we analysed the response of the system using a nonlinear, numerical evolution of the bulk Einstein-Maxwell-axion system of equations. First we studied finite-time quenches after which which the system is allowed to return to equilibrium. In particular no instabilities were encountered and the system behaved as expected, returning to equilibrium with the approach governed by the vector QNMs associated to momentum relaxation.
As previously noted, the vector QNM spectrum can exhibit qualitatively different behaviour depending on k. We showed the imprint of these features on the relaxation of J and Q. Most notably in a large k incoherent regime, Q receives parametrically enhanced contributions from a sector of QNMs which oscillate and decay, readily identifiable in the relaxation of Q following an electric field quench. We note that at fixed k/ √ ρ a transition from exponential decay to oscillatory decay of Q following a quench can also be obtained by varying T / √ ρ. These oscillatory QNMs generalise those studied at ρ = 0 where they are important for the thermal conductivity [28] in the vicinity of a coherent to incoherent JHEP10(2016)008 25 After some time the attractor behaviour is reached, marked by the black dashed line, which is the DC linear response conductivity after the equilibrium temperature is promoted to T E , as described in the text. The amount of heating before the linear response regime is reached increases with E f , leading to a higher T E value when they eventually agree. Here k/ √ ρ = 2.
transition. It would be interesting to understand whether the off-axis mode contributing to the oscillations in Q in the incoherent regime are in any way generic. We note that the situation described here resembles that of [46], where the exchange of dominance of purely imaginary and off-axis QNMs impacted the order parameter equilibration after a quench in a holographic model of superfluidity.
Next we studied a step quench to a constant electric field. For sufficiently weak, constant electric fields, an approximate steady state is reached described by linear response. Since the electric field was small but finite, Joule heating still occurs but is only relevant over much longer timescales. Going beyond linear response the effect of Joule heating becomes significant, and the energy density grows. Nevertheless after some initial response time, the electrical conductivity is well approximated by the equilibrium DC linear response result once the increasing effective temperature (as defined in section 5.1) is taken into account, σ lin .
The agreement of the nonlinear evolution with σ lin is somewhat surprising; in general we should have expected a dependence on E f . This is not unprecedented however, since at ρ = 0 we showed that σ = 1 at any E(t), just as in the case ρ = k = 0. It would be interesting to consider corrections to this approximation by performing an expansion in time derivatives of T E . It would also be interesting to perturbatively add charge to the ρ = 0 solutions presented in section 3. The problem considered here is inherently time dependent due to Joule heating. A natural next step is to prevent the energy in the system from growing without bound, so that we can reach a steady state at fixed temperature at late times. This may resemble a set up where the system is coupled to an external heat bath. Indeed this is the interpretation drawn in [37] for a steady state in the context of the D3/D5 probe brane system. A first glance at the right hand side of (1.4) suggests the possibility of preventing energy growth by turning on a time dependent source for an axion. Note that one example of a stationary solution with a time dependent scalar source is given already in section 2, by moving away from the present laboratory frame. Another possibility is considering an electric field applied to a strip or to a finite region on the boundary, leaving an infinite undriven heat bath.
Finally we note that the model considered here is a electrical conductor with an AdS 2 ×R 2 infrared scaling behaviour at T = 0. It would be interesting to contrast nonlinear transport for holographic metals and insulators across a variety of infrared behaviours. 10 Another interesting direction includes different spacetime dimensions, where σ is not di-10 See [47,48] for discussions of transport near quantum critical points.

JHEP10(2016)008
mensionless and can lead to qualitatively different behaviour [42]. Additionally it may be worthwhile to consider the nonlinear response of other classes of momentum-relaxing black branes, such as those which include a spatially modulated charge density. JHEP10(2016)008 using (A.5), getting F from (A.6), making sure to average it between times v and v + ∆v to ensure second order behaviour. The remaining equations are constraints; we ensure that they are solved for our initial data and consequently they are preserved for the evolution, provided we enforce energy conservation (1.4). As a note of caution, we have found that using a different choice for (A.6) in order to perform the evolution can lead to a numerical instability, which can seen by monitoring the constraints for sufficiently long times on an equilibrium solution. With the choice made in (A.6) the scheme is numerically stable and can be evolved for long times, and with good convergence as we show in appendix B.
We use Chebyschev collocation in the r-direction. It is convenient to factor out the leading terms at the AdS boundary, and work with Y ∈ {h vv , h vx , h s , h b , h, h v , h x }, which satisfy Y = 0 at the AdS boundary, imposed as a Dirichlet boundary condition. The function λ(t) is a gauge parameter, arising from the invariance of the form of the ansatz (A.1) under r → r/(1 + λ(t)r). For the majority of evolutions in this paper we use λ to adjust the radial coordinate so that the apparent horizon during evolution sits at r = 1. The only exception is the example shown in figures 2 and 3, where we fix λ = 0. For initial data we use the analytical black brane solutions presented in (2). Correspondingly we arrange our equilibrium initial data such that its horizon is situated at r = 1, in practise this is achieved by adjusting ρ for some fixed initial temperature, T i / √ ρ, and momentum relaxation parameter k/ √ ρ. This will allow us to perform large injections of energy without the need for excision. We do not impose any boundary conditions at r = 1. The undetermined data near the boundary associated to the holographic stress tensor and scalar vev is accessible via one radial derivative of the functions Y . In particular, after performing holographic renormalisation (see for example [50] and in the context of the axion model, [19]) T xx = −h (1) vv + 6h  with other entries zero. In addition these satisfy (1.1) and (1.2) courtesy of first order equations amongst the Y (1) ≡ ∂ r Y | bdy , as well as T µ µ = 0.

B Convergence testing
In this section we present convergence tests over the full dynamical range of one of the evolutions presented. A successful convergence test will show second order convergence with ∆v (for Crank-Nicolson) and exponential convergence with N , the number of Chebyschev points in the r-direction. Note that due to the exponential convergence with N , the residual becomes dominated by ∆v 2 -errors even for relatively small values of N , and inevitably it becomes impractical to reduce ∆v in order to overcome this. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.