Holographic pump probe spectroscopy

We study the non-linear response of a 2+1 dimensional holographic model with weak momentum relaxation and finite charge density to an oscillatory electric field pump pulse. Following the time evolution of one point functions after the pumping has ended, we find that deviations from thermality are well captured within the linear response theory. For electric pulses with a negligible zero frequency component the response approaches the instantaneously thermalizing form typical of holographic Vaidya models. We link this to the suppression of the amplitude of the quasinormal mode that governs the approach to equilibrium. In the large frequency limit, we are also able to show analytically that the holographic geometry takes the Vaidya form. A simple toy model captures these features of our holographic setup. Computing the out-of-equilibrium probe optical conductivity after the pump pulse, we similarly find that for high-frequency pulses the optical conductivity reaches its final equilibrium value effectively instantaneously. Pulses with significant DC components show exponential relaxation governed by twice the frequency of the vector quasinormal mode that governs the approach to equilibrium for the background solution. We explain this numerical factor in terms of a simple symmetry argument.


Introduction
The gauge/gravity duality applied within the context of strongly correlated many-body quantum systems started out as an interesting, yet limited, source of intuition on some properties of quantum critical matter. In the past decade, it has evolved into a powerful framework capable of taking into account a number of phenomenological aspects that should not be neglected when dealing with realistic models, such as crystal lattices, disorder, non-relativistic dispersion relations, etc. [1,2].
An important advantage of the holographic approach is its capacity of describing within a unique framework both equilibrium and out-of-equilibrium quantum systems by mapping them to tractable problems in general relativity, which can be systematically analyzed in real time without any need for conceptually new approaches. In the past, most of the attention towards far-from-equilibrium situations in this framework has gone to the JHEP07(2018)065 formation of quark gluon plasmas in heavy ion collisions. Holographic models relevant for this process incorporating a number of realistic features have been suggested and explored, with interesting results in relation with experiments [3]. Studies of far from equilibrium situations directly relevant to condensed matter systems have been, on the other hand, relatively scarce and mostly limited to toy models (see e.g. [4][5][6][7][8][9][10][11][12][13]). At the same time, recent advances in ultrafast experimental techniques in condensed matter physics have put a demand for a theoretical framework capable of explaining and predicting observed phenomena [14][15][16][17]. One is therefore led to ask whether, given the current state-of-the-art, time-dependent holographic models can make contact with experiments. In this paper we make a step in this direction by proposing a model for pump-probe experiments in which one follows the optical response of a holographic strange metal after it has been taken into a highly excited state by an electromagnetic pulse.
Our starting point is the minimal model considered in [18], which describes a 2+1 dimensional strange metal at finite temperature and density, in presence of a weak momentum relaxation mechanism obtained through axion fields linear in the boundary spatial coordinates. This efficiently reproduces the effects of explicit translational symmetry breaking [19] while preserving a homogeneous and isotropic bulk geometry (see also [20][21][22] for related holographic models). To mimic a pump pulse, we quench the holographic system by applying for a finite amount of time an oscillatory electric field. For simplicity we take it to be in the form of a modulated Gaussian wave packet of mean frequency ω P . In this way the system is driven into a highly excited out-of-equilibrium state, which then relaxes towards a new equilibrium state at a higher temperature, but equal charge density.
In contrast with the zero density case where, both with [23] and without [24] momentum relaxation the bulk dynamics results into a simple Vaidya geometry, at finite density the response of the system to the external electric field becomes more complicated. In fact, although the electric field always sets the charges into motion, explicitly breaking spatial isotropy and inducing on the boundary non-trivial currents, at finite density this also causes non-zero momentum densities. Holographically this corresponds to having additional metric components and field excitations, which generically make the problem not treatable analytically.
We study the resulting non-linear bulk dynamics with numerical methods and follow the evolution of the boundary one point functions as ω P is varied.
Although we work in the non-linear regime, we find that deviations from thermality after the pump pulse ends are surprisingly well captured within the linear response theory and their decay controlled by quasinormal modes (QNMs). In particular, at zero frequency the purely imaginary longest lived QNM of the vector sector governs the decay toward the new equilibrium configuration. 1 As the pump frequency is dialed up, we find that the response of the bulk geometry is increasingly well approximated by a bulk solution of the Vaidya form. That is, we observe that as soon as the pump electric field is turned off the boundary one point functions almost instantaneously approach their final equilibrium JHEP07(2018)065 configurations. In fact, in the limit ω P → ∞ we are also able to show analytically that the bulk solution takes precisely the Vaidya form, and one point functions thermalize instantaneously.
From the bulk point of view the origin of this dynamics can be understood from an analysis of QNM amplitudes. As we show explicitly, the amplitude of each QNM contribution is determined by the Fourier transform of the electric pump pulse evaluated at the frequency of the mode in question. The almost instantaneous approach to thermality at large frequencies is then explained by the absence of overlap between the pump spectrum and the frequency of long lived modes. A very simple toy model realized in terms of a driven harmonic oscillator effectively captures the main features of the bulk solution.
With the numerical background in hand we then proceed to compute the main observable of interest for a pump-probe experiment, the probe optical conductivity after the quench. In the same way as in a pump-probe experiment, we consider the optical response of the holographic strange metal to a probe pulse that is applied only after the pump has ended. This is incorporated in the definition of out-of-equilibrium conductivity we adopt.
Similarly to what happens for the background solution, the conductivity thermalizes almost instantaneously whenever the pump pulse has a negligible DC component. This behavior, although surprising from the boundary point of view, is completely natural with the insight provided by the analysis of the bulk background solution: if the geometry is described by a Vaidya solution, by causality in the bulk, the response to any perturbation applied after the light-like Vaidya shell will be insensitive to any detail of the quench other than the final equilibrium configuration. On the other hand, we find that for pump pulses with a DC component the optical conductivity relaxes with a rate set by twice the lowest vector QNM frequency. The appearance of this QNM, which governs momentum relaxation, can be understood from the fact that the zero frequency component in the pulse corresponds to a static electric field, which accelerates the finite density system. When the pulse is over, the resulting finite momentum has to relax in order to reach equilibrium. To explain the factor of two, which is less intuitive from a boundary point of view, we provide a careful but general analysis of linearized bulk fluctuations relying on the symmetries of the final equilibrium configuration.
A brief summary of our main results appeared before in [25]. There we proposed this model as an idealized setup for realistic pump-probe experiments, and the almost instantaneous thermalization as an extreme limit of fast thermalization that might manifest itself experimentally in certain regimes, similarly to what has been observed in the creation of quark gluon plasma in heavy ion collisions. In this paper, we present the computations behind them, as well as a number of new results, including the surprisingly good estimate of the size of non-thermality from a linear response analysis, and the explanation based on symmetry of the appearance of twice the lowest vector QNM frequency.
The rest of the paper is organized as follows. In the next section, we define the bulk model of a strange metal with momentum dissipation and briefly review its equilibrium properties. In section 3, details of the used numerical techniques are outlined. In section 4, we provide the non-equilibrium background solution computed numerically and discuss the behavior of the corresponding boundary one point functions. In section 5, we introduce JHEP07(2018)065 a toy model of a rapidly driven oscillator, which captures some of the important features of our holographic model and makes the phenomenological picture more transparent. Section 6 contains the main physical result of the paper, the time-dependent AC conductivity. Finally, in section 7 we conclude with a general discussion of our results, including prospects and challenges for comparison with experiment.

The model
The model we want to consider is specified by the action with Λ = −3 and equations of motion This admits a homogeneous and isotropic charged black brane configuration with nontrivial scalar fields profiles [24] that was explored in [18] as a simple holographic model for spatial translational symmetry breaking. Such a configuration has scalar fields with a linear dependence on the spatial coordinates x i = x, y common to the dual field theory and translationally invariant geometry and gauge field

4)
The dual field theory state is a thermal state with finite charge density and with translational symmetry breaking, whose properties can be fully specified in terms of T, ρ and k. The chemical potential µ is determined in terms of the charge density ρ by requiring the regularity condition that A should vanish at the horizon of the black brane, leading to µ = ρz 0 , with z 0 being the horizon location where f (z 0 ) = 0. The location of the horizon z 0 is associated with the temperature T of the field theory state through

JHEP07(2018)065
which gives the Hawking temperature of the black brane geometry. Notice that for fixed k, µ and z 0 , the mass parameter m appearing in the gravitational solution is not an independent quantity. It is fixed by the condition f (z 0 ) = 0 and is directly related to the energy density of the dual equilibrium state and the isotropic pressure Finally, the entropy density of this configuration is The reason why such a holographic solution with a completely homogeneous and isotropic geometry can be used to effectively describe momentum dissipation can be grasped from the Ward identities Following the standard AdS/CFT dictionary, the operators O I are dual to the bulk scalars and the couplings ϕ I are directly related to the asymptotic values of the bulk scalar profiles, that is ϕ I = kx i δ i,I in our case. Similarly the U(1) current J µ is dual to the AdS gauge field A µ and the boundary field strength F νµ is determined in terms of the asymptotic value of A µ . Let us first notice that (2.10) implies that the charge density ρ = J t is conserved. From (2.9) instead it follows that spatial momenta T ti will generically not be conserved whenever on the r.h.s. one has non-vanishing vevs. The coupling between the scalar and the gauge field in the bulk is such that the boundary electric field E i = F it induces a non-zero expectation value for O I , and thus Before concluding this section, let us quickly review the holographic result for the equilibrium optical conductivity computed in this model [18], which will be of use in the rest of the paper. The probe optical conductivity σ measures the linear response of the boundary current J x to a boundary probe electric field E x = −∂ t A 0 x . To compute it holographically one can consider the minimal consistent set of "vector" bulk fluctuations around the equilibrium background (2.3)-(2.4), and use the relation between the AdS asymptotic modes of δA x and the boundary quantities to write (2.14) The computation of the optical conductivity therefore amounts to solving the following system of linearized equations for the fluctuations (2.12) subject to appropriate asymptotically AdS boundary condition for the metric, nonvanishing source for the gauge field, and vanishing source for the scalar fluctuation.
Away from the zero-frequency limit these equations can be solved numerically, and one finds that for small enough values of k as compared to the other parameters of the equilibrium solution the resulting optical conductivity has the low frequency Drude form [18,20,27]. In figure 1, we reproduce a sample plot of the optical conductivity for k = 0.2, T = 0.2, and µ = 1.0. The finite value of the zero-frequency DC conductivity can be obtained analytically [18] The relaxation rate τ Q associated to the Drude peak corresponds to the purely imaginary frequency of the lowest lying quasinormal mode of the bulk vector preturbations. In the small k regime we will be interested in, this has been obtained analytically in [27] and reads

JHEP07(2018)065 3 Setup and details on numerical calculation
To study the response of the model of the previous section to the boundary electric field, we go to ingoing Eddington-Finkelstein coordinates and, following [13], consider the ansatz We solve the resulting system imposing appropriate asymptotically AdS boundary conditions under the assumption that for early enough times, when the pulse E x (v) has not been turned on yet, the solution coincides with the equilibrium configuration of section 2. By now, there are standard methods for solving such numerical relativity systems (see e.g. [28][29][30][31]). We will review the main ingredients of this procedure below. Inspecting Einstein's equations, it is convenient to define derivative operators along ingoing and outgoing radial null geodesics, which act on a field X(z, v) as follows: With this notation, the equations of motion become Notice that we are denoting, for example,Φ = ∂ z (Φ), i.e., the dot derivatives are taken before the z derivatives. The above set of equations represents a convenient set of nonredundant equations that can be obtained from all the non identically vanishing equations of motion following from our ansatz. The strategy for numerically solving the equations from the specified initial and boundary conditions proceeds iteratively as follows. The fields (B, Φ, a x ) represents the free initial data. All the other fields can then be solved from the equations of motion on a constant time slice. In more detail, after specifying the initial data at a given lightcone time, we solve (3.4) for the field Σ. This is a non-linear ordinary differential equation as no time derivatives appear. Next, equations (3.5) and (3.6) provide two coupled linear ordinary differential equations for F x and a v . Given Σ, F x and a v on the fixed time slice, one can proceed similarly to solve for the dotted fields. First we can solve (3.7) forΣ. Then, we can solve the two coupled linear differential equations (3.8) and (3.9) forḂ andȧ x and, subsequently, the linear ordinary differential equation (3.10) forΦ. Finally, via the linear ordinary differential equation (3.11) we determine F z . This way we obtain all the fields and dotted fields on the initial time slice. The time evolution is obtained by undoing the definition of the dotted derivative fields, to obtain a set of dynamical equations for (B, Φ, a x ) (3.14)

JHEP07(2018)065
Knowing their time derivative at a given time, we can time-evolve (B, Φ, a x ) to the next time step and we can repeat the above procedure to solve for all the fields on that time step. This way we can iterate the algorithm to the time we want. To solve the equations (3.4)-(3.11), we use the Chebyshev spectral method and introduce a Chebyshev grid z i in the z-direction. This way fields are replaced by vectors, X(z) → X i = X(z i ), derivative operators become matrices acting on the vectors and the differential equations translate into sets of coupled equations for the different field components X i . More precisely, the equation for Σ is non-linear and becomes a set of non-linear coupled algebraic equations for the coefficients Σ i , which we can collectively denote by The index j counts the number of components in the equation, which is the same as the number of variables Σ i . We solve this set of non-linear equations using the Newton-Raphson method. This method finds an approximate solution to f j = 0 as follows. First we start from a guess solution Σ (0) i . Then, use the updating routine Now Σ (1) should provide a vector which is closer to the solution of f j = 0 than our original guess. Repeating the algorithm by taking Σ (1) as a new guess and using (3.16), we again get closer to the correct solution. Iterating this algorithm many times, one should converge to the solution of f j = 0. In practice the number of iterations needed depends on how good the initial guess was. In our numerical algorithm we use the solution from the previous time step as the initial guess. This way we only need a few (typically 3) iterations to solve the equation of motion to the desired accuracy (around 10 −13 accuracy). The rest of the equations (3.5)-(3.11) are all linear in the unknown variables and can be straightforwardly solved by standard matrix inversion methods. We have used the numpy.linalg.solve and numpy.linalg.inv functions, which are included in the Python Numpy package and are based on the LAPACK library. At the practical level, when solving (3.4)-(3.11), in order to simplify the numerics, we also find it convenient to subtract or rescale the near boundary behavior of some of the fields. In particular, we work with the set of regularized fields X r defined as follows

JHEP07(2018)065
The functions F x , a x , a v are left intact. Note that in these redefinitions,Φ r is not the dot derivative acting on Φ r , but a new variable defined through the equation in (3.18). The same applies for all the other dotted fields. In addition to (3.4)-(3.11), there are other components of Einstein's and Maxwell's equations that do not vanish identically for our ansatz. These are in principle redundant with (3.4)-(3.11), but in practice they are useful for testing our numerical solutions.
To evolve (B, Φ, a x ) we use the fourth order explicit Runge-Kutta method. The time domain is divided in discrete time-steps v n and the value of a field at step n + 1, X(v n+1 ), is obtained as the value at step X(v n ), plus the weighted average of four different time increments determined in terms of (3.12)-(3.14).
Using this algorithm, we can solve the full numerical problem modelling the pump probe experiment. In practice we have found it computationally faster to separate the problem of determining the probe conductivity as a separate problem. Thus, we use the above algorithm to solve for the spacetime corresponding to the system subject to the pump electric field. To obtain the probe conductivity, we linearize the equations of motion around the numerically known background spacetime and solve them using a similar numerical procedure as above. The main advantage of this procedure is that when solving the linearized equations of motion, we do not need to use the Newton-Raphson method, but all the equations are solved using linear algebra, which is faster. This becomes particularly useful when we consider several probe "experiments" for the same pump pulse. We have checked the numerical accuracy of solving the linearized system by comparing the results to those obtained using the (slower) full code for both pump and probe parts. Further checks are provided by testing the system on the equilibrium states, in particular by comparing the numerically obtained conductivity with the analytic formula for the DC conductivity. These agree to a very good accuracy (in the cases we have tested they agree up to 10 −7 % accuracy).

Numerical error estimate
There are two sources for the numerical error in our procedure. The first one arises from discretizing the z coordinate, and the second one arises from discretizing the time coordinate. As a measure of the numerical error we use the remaining three redundant equations of motion. For an exact solution, these equations would be automatically solved. Denoting the equations as Eq i = 0 where i = 1, 2, 3, we consider the following quantity We evaluate the equations on the spacetime grid using fourth order finite differences for approximating time derivatives and then find the maximum value of |Eq i | within the grid. Finally we take the root mean square of the maximum error of the three equations. This error measure is displayed in figure 2 as a function of the number of timelike N t and spacelike N z lattice points for fixed timelike and spacelike size of the computational domain. From this error measure, we find that the numerical error approximately first decays as Err ∝ N −4 t as N t is increased and then saturates to a constant value. This is expected  Figure 2. Numerical error as a function of the number of timesteps N t . The different curves correspond to different numbers N z of spatial lattice points. Here we study a Gaussian pulse with the choice of parameters: A = 0.5, t 0 = 3, ∆t P = 1, as there is a remaining error due to finite number of spatial lattice sites N z . Increasing this number then decreases the saturated value approximately exponentially. Thus, as both N t and N z are increased the error is found to decrease rapidly, which gives strong evidence that the numerical calculation is converging towards a solution of the continuum equations of motion. 2 In practice we have found that rather small values of the spatial lattice sites such as N z = 8 or N z = 10 are sufficient to give reliable results. For N t we use the highest values from those shown in figure 2. This is forced by the fact that the probe pulses have to be short in order to reasonably approximate delta functions. On the other hand the timelike computational domain has to be large in order to get a reliable Fourier transform of the differential conductivity, without finite size effects. For example for a computational domain of length of order 10 3 we use N t of the order 10 6 . This results in a computational time of the order of tens of hours on a laptop.

Non-equilibrium background spacetimes
To model the process of applying the pump electric field, we start from an initial state corresponding to an equilibrium black brane dual to a state at a given temperature T I .

JHEP07(2018)065
The time dependent pump electric field then takes the system out of equilibrium to a configuration captured by the ansatz (3.1) to finally reach a new equilibrium configuration at a different temperature T F . Throughout this process we keep k fixed and ρ is conserved, as guaranteed by Ward's identities. In more detail, the starting equilibrium configuration in terms of the regularized fields defined in (3.18) corresponds to setting and all the other regularized fields in (3.18) to zero. The parameters m I and µ I are determined in terms of T I , k and ρ according to the relations of section 2. The specific form we use for the pump field is given by This represents a Gaussian wavepacket with central frequency ω P and width ∆t, centered at t 0 and cut off by a smoothed step function at t end ≡ t 0 + 3∆t (from which time onwards we consider the pumping to have finished). Throughout this section, we choose the parameters t 0 = 50, ∆t = 15 and δ = 0.01. The pulse amplitude A is instead tuned in order to obtain the desired increase in temperature. At the time of the pulse, the metric functions start time-evolving, exciting all the rescaled fields defined in (3.18). Notice that in particular, this will give a nontrivial expectation value for the boundary operators associated to the bulk fields. More specifically, according to our ansatz, the current J µ associated to the bulk gauge field, the operator O associated to the scalar excitation Φ and the non-isotropic stressenergy tensor T µ associated to the bulk metric will acquire a time dependent expectation value. At late times they will all settle to new equilibrium values with and again all other rescaled fields vanishing. The final mass parameter of the black hole will increase throughout the process, m F > m I , consistently with the fact that energy has been pumped into the system. As an example, figure 3 shows plots of some metric function components obtained from the numerical solution.
To obtain boundary theory expectation values from the bulk solution, one has to perform the corresponding holographic renormalization procedure [18]. The resulting one point functions are given in terms of asymptotics of the bulk fields as [13] = T tt = −2F z,r , where the different bulk function appearing on the r.h.s. are all evaluated at the AdS boundary z → 0.

Vanishing pulse frequency
We start by considering the particular case where the pump field is not oscillating, ω P = 0. Figure 4 shows the boundary theory expectation values for the same type of time dependent state represented in figure 3.
In particular one can observe that the pump electric field E x (t) induces an electric current and a momentum current in the field theory. Furthermore, there is a substantial pressure anisotropy induced and, in the case represented here, the energy density increases by more than a factor of two. All one point functions, except for the energy density, seem to have a relaxation time far longer than the time scale of the pump field. Inspecting the logarithmic plots in figure 5 one finds that they decay towards equilibrium exponentially in time. The rates of the exponentials are consistent with where ω * = −iω i is the, purely imaginary, lowest quasinormal mode in the vector channel.
It is important to note that the pressure anisotropy is decaying with double the rate of the other expectation values. We are able to provide an explanation for this, working under the reasonable assumption -supported by the numerical calculations -that the deviations from thermality are sufficiently small at late times, so that the equations of motion can be expanded in powers of the deviations. For this, consider an expansion around the final thermal black brane configuration (4.9) From this we can argue that the decay rate of the source term sets the decay rate of the B field, which is thus twice the decay rate of the vector perturbations. That is, twice the imaginary part of the lowest vector quasinormal mode. This explains the factor of two in the decay rate we see from the numerics in figure 5.

Increasing the pulse frequency
So far we have studied the case of an approximately Gaussian pump electric field (4.2) with a vanishing mean frequency ω P = 0. This has reproduced the by now standard story that the late time relaxation of the black brane solution is dominated by the lowest quasinormal mode that gets excited. A slight subtlety was that some of the metric components relax with a rate given by twice the lowest quasinormal mode from a different sector. Next, we will study how this picture changes as we increase the mean frequency of the pump pulse. In figure 6  Furthermore, the magnitudes of T tx , (T xx − T yy ) , O are decreasing with increasing ω P , while the magnitudes of T tt and J x stay fixed. We note that, strictly speaking, even if the leading QNM has only infinitesimal amplitude, one could still choose to refer to its decay constant as the decay time. If the amplitude of the QNM is below the experimental resolution, however, then it is not measurable, and in that sense irrelevant. In this paper, we therefore refer to thermalization as instantaneous or very fast if the slower decay modes have zero or negligible amplitude.
These observations suggest that the spacetime could be approximated with the Vaidya spacetime, with an appropriately chosen time dependent mass function, at large enough ω P . In the rest of this section we provide evidence in support of this claim. First we will show that working in the limit of a large pulse frequency the leading order solution is exactly of the Vaidya form. We obtain this result working analytically in the large frequency expansion. Next, we analyze the amplitude associated to the quasinormal mode decay described above to show how this is determined by the relation between the power spectrum of the pump pulse and the quasinormal mode frequency.

Large frequency solution
In the regime where the pump frequency is very large compared to the other parameters of the gravitational background the bulk solution can be studied analytically. 3 We assume the electric field is of the simple oscillating form E x (t) = cos(ω P t)Ω(t) , (4.10)

JHEP07(2018)065
where the enveloping function Ω(t) is assumed to have compact support and slow variation compared to the cosine. Using the knowledge obtained from the numerical solution and inspecting the equations of motions, we formulate an ansatz for the 1/ω P expansion of each field and for the type of time dependence (rapidly or slowly varying) for each term in the expansion, and proceed to solve the resulting system of equations order by order. The details of the analysis are reported in appendix A. For the different fields and metric components, the leading correction to the unperturbed solution induced by the rapidly varying source E x takes the form At leading order in the frequency expansion only F z gets corrected by 4 This directly shows that the response to the rapidly oscillating electric field takes at leading order the Vaidya spacetime form with the mass function M (v) given by the background m value plus the contribution coming from F (4.14) The first correction to the Vaidya form of the geometry comes from the F x component of the metric at order 1/ω P In the limiting approximation where one can treat the function Ω(t) as a constant under the integral, we would simply have F x ≈ ρz sin(ω P v)Ω(v)/(3ω P ). Notice however that in our case for those times v where E x (v) has no support, that is times where the pump pulse has been turned off, the suppression of this correction is even stronger. In fact, with a 4 Strictly speaking here and in the expression (4.14) below we are assuming that on the r.h.s. we are consistently taking only the leading contribution from the integral 1 2 v −∞ dv Ex(v ) 2 , which in general will also have subleading terms in 1/ωP (see appendix A). choice of E x (v) of the form (4.10) and for Ω any smooth function, F x is suppressed more strongly than any inverse power of ω P , as follows from basic Fourier analysis. At order 1/ω 2 P also Φ and a x get their leading contribution from the pulse, which we report here for comparison with the result obtained from the numerics while the order 1/ω 3 P gives the leading corrections to the background values of a v and B small value of ω P = 0.5. Thus, the two results do not agree quantitatively very precisely, although the qualitative form of the solutions is already very similar. Furthermore for T tx , the difference between the full solution and the approximate analytic one is already surprisingly small. To test the convergence of the approximate analytic solution to the full numerical solution at large ω P , we define the subtracted one point functions These are plotted in figure 8, where we have multiplied them with appropriate powers of ω P in order to make the corresponding expectation values order one in the large ω P limit. It can be seen that the full numerical solution is converging well to the approximate analytic one.

Estimating the size of non-thermality from linear response theory
In this subsection, we estimate the size of the quasinormal mode contributions to the one point functions (and thus, to the gravitational background solution) from linear response theory. That is, we assume that the electric field E x is small and we evaluate at linearized level its effect on one point functions. This assumption is clearly not a priori valid for our setup, but the final result we obtain this way is surprisingly close to the exact numerical results. Linear response theory tells us that the leading contribution to the expectation value of an operator χ(t) due to the presence of an external electric field is given by 5 where The expectation value is taken in the final equilibrium thermal state, and in writing (4.24) we have taken into consideration the fact that we are considering a spatially homogeneous configuration and consequently G χ,Jx R is independent of the spatial position. We will be considering the operator χ to be one of the vector sector T tx , J x , O.
The retarded correlator can be expanded at late times (i.e. large |t − t |) in terms of quasinormal modes where ω n are the quasinormal modes frequencies shared by the correlators of the vector sector operators and g n are the residues of the quasinormal mode poles in the Fourier transformed correlator.
In the situation we are interested in, A x vanishes after the pulse is turned off. For late enough times t compared to the time where the source has been turned off, we can reliably substitute the quasinormal mode expansion inside the integral (4.23). This way we obtain the late time expression (4.26) Thus, χ(t) decays at late times with a rate set by the quasinormal modes and an amplitude set by the integral involving A x . When t is large enough to be outside the support of A x , the time integral with upper limit is formally equal to the integral over the entire temporal domain. Integrating by parts we can express this in terms of the electric field 5 Here for simplicity we work in the equivalent gauge where the electric field is generated by Ax, that is Ex = −∂tAx. In writing the linear response (4.23) we are also assuming that the expectation value of the operator χ(t) is zero when the electric field is absent, which is the case for the operators we will be considering.  which gives for the linear response of the operator χ at late enough times where the electric pump field has been turned off This shows how the amplitude associated to each quasinormal mode contribution is determined by the Fourier transform of the pulse electric field evaluated at the frequency of the quasinormal mode itself. At the times we focus on in our analysis, the lowest lying QNM with purely imaginary frequency generically dominates over the others. There can be in principle cases where other QNMs, with frequency having a non-vanishing real part, may be in resonance with the pump pulse and compete with the leading QNM. Nonetheless, these contributions would decay fast in time and quickly become negligible. In figure 9 we plot the expectation value T tx at the time t = t end when the pump pulse turns off. In the large ω P approximation the expectation value is immediately zero at this time. In the numerics we see significant deviations of T tx from zero for small ω P . The deviation mainly arises from the lowest quasinormal mode contribution, whose size is given by the Fourier transform (4.27) in the linear response approximation. The blue curve in figure 9 is obtained by fitting (4.27) with a constant coefficient in front as a fitting parameter. As the figure shows, the linear response curve fits the numerical data very well with a root-mean-square error 6 of RM SE ≈ 0.0036. Similarly, also the other one point functions J x and O can be fitted with the form (4.27). An exception is the one point function (T xx − T yy ) which is better fitted to (E x (ω = ω * )) 2 . In the right part of for (T xx − T yy ) the relation is approximately quadratic. The quadratic approximation appears for the same reason as the factor of two in the decay of (T xx − T yy ) , as the corresponding field is sourced by the squares of the vector sector fields.

Driven oscillator toy model
As we will show in this section, two signatures in the evolution of our holographic setup can be captured and explained by a simple toy model given in terms of a driven damped harmonic oscillator. These are the instantaneous relaxation at large driving frequency ω P → ∞ and the smallness of the quasinormal mode contributions when ω P is separated from the real parts of the quasinormal mode frequencies. The equation of motion of the driven oscillator is given byẍ where ω 0 is the undamped oscillation frequency and κ the damping strength. We will choose the driving force to have the form where Ω(t) varies much more slowly than the cosine. Without a driving force, displacements of x decay back to zero exponentially in time as These complex frequencies represent the analogue of the quasinormal mode frequencies in our system. The driven system (5.1) is solved explicitly as where the time dependent amplitudes associated to each mode are In writing the solution we have assumed that x(t) vanishes at early times before the driving force has been turned on. Let us consider the large ω P limit of the solution. The basic intuition is that in this case x(t) oscillates fast in which case the equation of motion can be approximated as x(t) ≈ f (t), with the approximate solution x(t) ≈ − cos(ω P t)Ω(t)/ω 2 P . We can see this from the exact solution using integration by parts twice, 7 leading to

JHEP07(2018)065
Substituting into (5.5), we arrive at Thus, we see that the oscillator follows the driving force instantaneously and, in particular, it relaxes instantaneously as the force vanishes. This behavior is similar to the instantaneous thermalization of Vaidya spacetimes where one point functions relax to their thermal values as soon as the boundary source turns off. Notice also that the series expansion in 1/ω P can be computed to arbitrary orders by repeatedly integrating by parts. This way the instantaneous behavior is seen to hold to all orders in the 1/ω P expansion. The next question we want to address is that what happened to the quasinormal mode contributions and whether a large ω P is a necessary condition to have instantaneous thermalization.
At late times, that is for times after the driving function has been turned off, it is apparent that the solution written in the form (5.5) takes the "quasinormal mode" decay form (5.3). In fact, A ± (t) become time independent as long as one considers values of t where the driving force has been turned off, and therefore at late times give the amplitude associated to the quasinormal modes. For concreteness, let us assume that Ω(t) has support on a compact region so that the driving function is turned off after some time t end . For all those times t > t end where the driving force has been turned off, one can formally replace the integrals in (5.6) with integrals over the entire time range, that is where we denote the Fourier transforms aŝ For any smooth choice of Ω(t), the coefficientsΩ of the quasinormal modes will be suppressed more strongly than any inverse power of the arguments ω ± ± ω P . To build some more explicit intuition, we can specialize to an example close to our holographic calculation by choosing a Gaussian envelope Strictly speaking with this choice the forcing pulse is never turned off. However, for all practical purposes, at large enough times compared to the Gaussian width we can consider the driving force to be vanishing. With this choice, the amplitudes of QNMs take the form Notice that, modulo a factor 2ν, these are nothing else than the Fourier transform of the driving forcef evaluated at the QNM frequencies p = ω ± . This shows explicitly how the amplitudes associated to QNMs depend on the overlap between the spectrum of the driving force and the real part of the QNM frequencies: if the driving frequency ω P coincides with the real part of a QNM frequency there will be no Gaussian suppression of the amplitude of the corresponding mode. Conversely, for Re(ω ± ±ω P ) = 0 the QNM amplitude is exponentially suppressed in the square of this combination. Let us further specialize to the over-damped case, where κ > 2m. This exactly mimics our holographic setup in the fact that the relevant QNM frequencies are purely imaginary.
Since ω ± have no real part, the amplitudes of the QNM excitations are exponentially suppressed in ω P as (5.14) Thus, for large ω P ∆t the quasinormal mode contributions are very strongly suppressed.

Out of equilibrium conductivity
Next, we want to study the conductivity properties of the non-equilibrium states we have prepared with a pump pulse. For this purpose one introduces a smaller "probe" electric field δE x (t). This electric field induces a change in the current δ J x (t) . This way we can define a real-time conductivity called differential conductivity σ(t, t ) through the relation Although there is no standard definition of frequency space conductivity out of equilibrium, we will follow [33] and define σ(ω, t) = tm t dt e iω(t −t) σ(t , t) , where t m is the time at which the experiment ends. The conductivity defined this way can be related to the current two point function as discussed in appendix B. In thermal equilibrium, the conductivity (6.2) approaches the standard definition of optical conductivity at frequency ω, and spatial momentum k = 0 as the observation time t m is sent to infinity. This is also shown in appendix B. In practice we calculate the conductivity in two steps. First, we calculate the differential conductivity appearing in (6.1), and then perform a Fourier transform to obtain (6.2). By choosing a probe field δE x (t) = δ(t − t 0 ), (6.1) becomes This way we obtain the differential conductivity directly from the knowledge of δ J x (t) .
In the numerical implementation, we actually use a smoothed version of a delta function, which we choose to be a Gaussian Re(σ) a δt = 0.00 initial thermal final thermal Figure 10. The optical conductivity at different times and in the initial and the final thermal equilibrium states. Left: the pump frequency is ω P = 0.5, and the optical conductivity right after the pump pulse has ended coincides with the final equilibrium one. Right: the pump frequency is ω P = 0, and the optical conductivity is seen to interpolate in time between the initial and final equilibrium values.
In the limit δt → 0, this approaches a delta function, while in practice we keep δt nonvanishing but small enough that it does not affect our results considerably. Smoothing out the delta function over a small scale δt affects the conductivity at frequencies ω ∝ 1/δt and larger, while the main interesting time-dependence in the conductivity is seen for frequency ω = O(1). In our analysis, we have taken δt = 0.05, which has only small (order 1%) effect on the conductivity in the regime of interest.

Numerical results
First, we consider the pump pulse profile (4.2) with different values of the pump frequency ω P . Figure 10 shows the optical conductivity as a function of ω for different values of times. The time δt is measured from the time t m = 3∆t + t 0 at which the pump pulse practically turns off (due to the smoothed theta function in (4.2)), Figure 10a shows the conductivity for the large frequency pump pulse which appears thermal immediately at time δt = 0. This is consistent with the results of the previous section which showed that the larger ω P , the closer to the Vaidya spacetime we get. On the other hand figure 10b shows the conductivity for ω P = 0, in which case the conductivity deviates significantly from the final thermalized conductivity for all times displayed. Next, we study how the conductivity approaches its final thermalized value. We will focus on the DC conductivity, that is, the optical conductivity σ(ω) at ω = 0. Figure 11 shows σ DC as a function of time. In the previous section we saw that the background spacetime approaches a static black hole with a rate set by the lowest quasinormal mode. Thus, we might expect that the conductivity approaches its final thermalized value with the same rate. This is where we find a slightly surprising result. The conductivity approaches its thermalized value with a rate 2Im(ω * ). We will come back to this factor of two in the next subsection. Thus, the DC conductivity is consistent with the approximate form where σ th DC is the final thermalized value of the DC conductivity. The coefficient C quantifies how far out of equilibrium the conductivity gets. Defining with the approximate form (6.6) we have δσ DC = −C.
In figure 12 we show how δσ DC behaves as a function of ω P . As before, the initial and final temperatures are kept fixed while varying ω P . Clearly the largest deviation appears as ω P = 0. Recalling that the amplitude of the deviation of the background spacetime from the equilibrium one was well approximated by the Fourier transformed electric field E x (ω * ), we are motivated to also attempt a similar fit to the deviation of the conductivity from its equilibrium value. Fitting δσ DC with E x (ω * ) does not give a good fit but instead E x (ω * ) 2 does. A fit to E x (ω * ) 2 is shown in figure 12 as the solid green curve. This is a one parameter fit with the overall amplitude being the only parameter. Thus, our results for the time dependence and ω P dependence of the non-equilibrium conductivity can be summarized in  Figure 12. The maximum deviation of the DC conductivity from its thermalized value as a function of the pump frequency ω P . The solid green curve is a one parameter fit to a form (E x (ω * )) 2 Finally we study how the deviation from thermality behaves as we change the difference between the initial and final temperatures. We have chosen to keep the final temperature T f fixed and to vary the difference ∆T = T f − T i by changing the initial temperature T i . A plot of δσ DC as a function of ∆T is shown in figure 13. This confirms the intuition that the more we increase the magnitude of the pumping electric field, the further from equilibrium the conductivity deviates (in units of the final thermalized conductivity σ th DC ).

A symmetry argument for the thermalization rate
We have given numerical evidence that the conductivity thermalizes with a rate e −2iω * t where ω * is the lowest vector quasinormal mode. The conductivity deviates from thermality simply because the background spacetime deviates from a static black hole. Thus, it is not surprising that the thermalization time scale of the conductivity is related to that of the background spacetime. What is somewhat surprising instead is the factor of 2 relating the two rates. Here we provide a symmetry argument for it.
Let us start considering the symmetries of the bulk action (2.1) of the holographic model we are studying. These include the subgroup SO(2) × SO(2), where the first factor represents the rotations M acting on the spatial coordinates x i = (x, y) common to the boundary, and the second factor the global rotations R that act non-trivially only on the scalar duplet φ I = (φ 1 , φ 2 ), rotating the two fields into one another. The equilibrium solution, and more specifically the scalar field configuration (2.3) explicitly breaks this symmetry as SO(2) × SO(2) → SO(2) res , with the residual SO(2) being the subgroup that leaves the scalar field configuration invariant, The bulk metric and the bulk gauge field configurations (2.4) do not break any of the original SO(2)×SO(2) symmetry, as they are isotropic (and homogeneous) in the boundary spatial coordinates and they do not transform under the global SO(2) symmetry. Hence the residual SO(2) is automatically preserved by the bulk metric and gauge fields, and the entire equilibrium solution is a scalar of SO(2) res . One can conveniently organize deviations from the equilibrium background solution according to representations of the SO(2) res . Starting with the ansatz (3.1) for the nonequilibrium solution, and employing the definitions given in (4.6), we can split the fluctuations as δF z , δΣ, δa v scalar, δF x , δa x , δΦ vector, δB symmetric traceless tensor.

JHEP07(2018)065
As we described, at sufficiently late times the spacetime is close to thermal equilibrium and one can to a good approximation expand the equations of motion in perturbations around the equilibrium spacetime. At the linearized level the three type of perturbations completely decouple from each other. This can be immediately understood by thinking about the equation of motions in terms of the SO(2) res symmetry. In general, the linearized equations for the fluctuations can be schematically written as T T δT = 0.
Here δS, δV, δT collectively indicate fluctuations belonging the scalar, vector and traceless tensor sector respectively. L XY indicates operator constructed from the equilibrium solution and derivative operators, which acting on linearized fluctuations establish a map from the sector Y to the sector X. Since we only consider spatially homogeneous background fields and fluctuations, if follows that the operators L effectively do not contain any derivative operator in the boundary spatial coordinates -or rather these have a trivial effect. This, together with the fact that the equilibrium solution completely belongs to the singlet representation of SO(2) res , implies that the operators L T T δT = 0 , (6.13) as we explicitly discussed in section 4.1. To linearized level the only non-vanishing perturbations we considered were the vector fluctuations, which decayed towards equilibrium with a rate q = e −iω * t set by the lowest vector quasinormal mode ω * . Working at the quadratic level in the fluctuations, the equations of motion can now be written according to the structure given by SO(2) res in the form 14) T T δδT = J T .
Now δδX indicates quadratic fluctuations and, in a similar way as above, L XY indicates an operator constructed from the equilibrium solution and derivative operators. The sources J X collect terms that are quadratic in the linearized fluctuations δS, δV, δT . Using the same symmetry argument as above, together with the fact the we only have vectorial perturbations at the linear order, we can conclude that the equations for the quadratic fluctuations relevant to our case take the form T T δδT = J T . (6.15)

JHEP07(2018)065
Again, this consistently reproduces what we observed in section 4.1. In particular, as the leading excitations of the scalar and tensor sector are sourced by expressions that are quadratic in δV ∝ q, we have that δδS, δδT ∝ δV 2 ∝ q 2 . All in all, this analysis teaches us that the non-equilibrium solution written in a perturbative expansion has the schematic structure with S ∝ q 0 indicating the equilibrium, scalar sector, solution. Given this, we can proceed to analyze the linearized perturbations that compute the optical conductivity in the background (6.16). We indicate these new sets of linearized fluctuations collectively asδs,δv,δt to distinguish them from the background ones. Similarly we useL to indicate the operators that act on these fluctuations in the equations of motion, which take an analogous form as (6.12).
We are interested in how much these fluctuations deviate from the form they would take around the background equilibrium, O(q 0 ), solution. For this we organize our computation in an expansion in q. Making this perturbative structure explicit, we have up to second order in q included L ss ≈L ss,0 +L S ss,2 δδS,L sv ≈L V sv,1 δV ,L st ≈L T st,2 δδT , L vv ≈L vv,0 +L S vv,2 δδS +L T vv,2 δδT ,L vs ≈L V vs,1 δV ,L vt ≈L V vt,1 δV , (6.17) L tt ≈L tt,0 +L S tt,2 δδS ,L ts ≈L T ts,2 δδT ,L tv ≈L V tv,1 δV .
Similarly we write the fluctuations asδs =δs 0 +δs 1 + . . . , withδs i being order q i , and analogously for the other sectors. Solving the resulting equations order by order in q, the zeroth order problem corresponds to the study of linearized fluctuations around the thermal equilibrium value. All three sectors remain decoupled and to compute the optical conductivity one only excites the vector sector, so onlyδv 0 is non-vanishing.
At the next order in q we havẽ L XY,0δ Y 1 +L XY,1δ Y 0 = 0, (6.18) and using the fact thatδs 0 =δt 0 = 0 and the explicit form ofL XY,1 that can be read form above, it is easy to see that the equation forδv 1 reduces tõ Thus, there is no order q contribution to the optical conductivity. At the next order instead 20) and generically all the sectors get sourced by the lower orders solution. Concentrating oñ δv 2 , which is the relevant sector for the optical conductivity,

JHEP07(2018)065
We therefore conclude that for the vector fluctuations the leading deviation from their thermal value is of order q 2 . That is, we showed that the rate of thermalization of the optical conductivity corresponds to e −2iω * t .

Discussion
In this paper we have analyzed the pattern of conductivity thermalization in a minimal holographic setting that includes finite charge density and a mechanism of weak momentum dissipation. We have shown that, when quenched by a laser pulse with a significant DC component, the equilibration time of the conductivity is given by half the lowest-lying imaginary bulk quasinormal mode: τ = −1/(2Im ω * ). The appearance of the QNM governing momentum relaxation can be understood from the fact that the DC component of the electric field sets the charges of a finite-density system in motion, thus inducing finite momentum densities in the system. We have also provided a symmetry argument for the factor of two. If the mean frequency of the wave packet is large, the pulse lacks resonance with the corresponding QNM, and the thermalization is effectively instantaneous. The latter result is surprising from different points of view and gives rise to a number of questions. The instantaneous thermalization of the conductivity is surely not intuitive from the boundary perspective, but is somewhat natural from the bulk point of view once we are given a background dynamics of the Vaidya form. However, when going to finite density, a priori we would not have expected the bulk dynamics to remain (even approximately) of the instantaneously thermalizing Vaidya form as in the zero density case [23,24]. It is then interesting to ask to what extent this result is generic and, especially, whether it still holds in non-relativistic systems. A natural option would be to consider Einstein-Dilaton-Maxwell holographic models that explicitly exhibit Lifshitz scaling and hyperscaling violation, which have also been generalized to include momentum dissipation [34]. However, a closer look reveals that they do not allow for a dynamical electric field, as this will generically result in a modification of the dilaton profile and affect the scaling exponents. On the other hand, models that are Lorentz-invariant in the UV and where the relevant scaling properties only emerge in the IR (see, e.g., [35][36][37]) seem to be less constraining in this sense, and could represent an interesting context where to explore this question further.
From a broader perspective, our observation poses a question about the underlying principles governing such ultrafast equilibration. One might ask to which extent this surprising prediction depends on the precise details of the UV description of the field theory, and more specifically what is the role played by the large N limit, on which our classical gravity computation relies. Recently it has been shown that the causal behavior of certain observables associated with the Vaidya geometry, in the zero density case, can be explicitly recovered from conformal field theory structures in the limit of infinite central charge [38]. It has also been suggested that 1/N -corrections deflect the system from this regime [38]. This might suggest that the Vaidya-like response of the holographic strange metal to the laser pulse is simply an artifact of the regime of classical gravity. The eigenstate thermalization hypothesis provides us with another way to think of this instantaneous equilibration, relating rapid thermalization of local observables to the dense entanglement JHEP07(2018)065 within the full many-body quantum system [39,40]. From this perspective, it is possible that the seemingly universal behavior of holographic systems is not an artifact of specific limits, but rather a generic manifestation of the entangled nature of those states that have holographic duals.
In the same way that the holographic predictions have been shown to enjoy some UV independence through the minimal viscosity [41] and the early onset of the hydrodynamic behavior [28] that are reflected in heavy ion collisions experiments [42], it would be interesting to look for signatures of this instantaneous thermalization in condensed matter systems. This is the perspective we highlighted in the companion paper [25].

A Large ω P solution
In order to work out the solutions in a 1/ω P expansion, based on what we observed from the numerical results for the bulk solution, we formulate the following ansatz for the expansion of the different metric components

JHEP07(2018)065
The response to the time dependent pulse order by order can have contributions with or without large frequency. This is encoded in the two terms in the last sum in (A.1). The functions X n appearing in the expansion are assumed to vary slowly in time, whileX n are quickly varying functions and bring with them a factor ω every time they are hit with a time derivative. To make this manifest, we explicitly included a "time derivative counting factor" ω P in their argument, which will eventually be set to one. Notice that when solving the resulting system of equations order by order, X n will enter at the same order as single time derivatives of the fieldsX n+1 . We will be working under the assumption that in the pulse (4.10) the enveloping function Ω(t) has always negligible variation. At leading order the non-identically vanishing and linearly independent equations are the E vv , E vz and E vx components of Einstein's equations, which read The equation forF x,0 is readily solved as Similarly for F z we have Notice that the integral on the r.h.s. will in general give an order ω 0 P contribution and an order ω −1 P contribution as reflected on the l.h.s. of the equation. That is (A.14) For later convenience we define the order ω 0 P quantity ∂ v Φ x,0 (z, ω P v) − ∂ zΦx,0 (z, ω P v) + kz 2 F x,0 (z, ω P v) = 0 , (A.16) 2 ω P ∂ v ∂ zãx,0 (z, ω P v) − 3ρzF x,0 (z, ω P v) = 0 , (A.17) 4 ω P ∂ v 2F x,1 (z, ω P v) + z∂ zFx,1 (z, ω P v) + (3k 2 z + 2z 3 ρ 2 )F x,0 (z, ω P v) = 0 . In thermal equilibrium σ(t, t ) can only depend on the difference t − t due to time translational symmetry of the thermal density matrix, so that σ(t, t ) = σ(t − t , 0). Changing the integration variable t into t = t − t , the two integrals separate and we obtain The second term on the righthand side can be recognized as the frequency space nonequilibrium conductity (6.2) evaluated at t = 0. Thus, we see that in thermal equilibrium it satisfies which is the standard definition of optical conductivity in thermal equilibrium. Since σ(ω, t) is independent of t in thermal equilibrium, (B.8) holds for arbitrary t. This establishes the claim that the non-equilibrium conductivity reduces to the equilibrium conductivity in thermal equilibrium and as t m → ∞. If t m is finite, the step that fails above is the factorization of the integrals after the change of variables t = t − t , since this changes the integration limits. But as long as one considers values of t t m , the integral defining σ(ω, t) can be safely extended to infinity.
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.