Time-dependence of the holographic spectral function: Diverse routes to thermalisation

We develop a new method for computing the holographic retarded propagator in generic (non-)equilibrium states using the state/geometry map. We check that our method reproduces the thermal spectral function given by the Son-Starinets prescription. The time-dependence of the spectral function of a relevant scalar operator is studied in a class of non-equilibrium states. The latter are represented by AdS-Vaidya geometries with an arbitrary parameter characterising the timescale for the dual state to transit from an initial thermal equilibrium to another due to a homogeneous quench. For long quench duration, the spectral function indeed follows the thermal form at the instantaneous effective temperature adiabatically, although with a slight initial time delay and a bit premature thermalisation. At shorter quench durations, several new non-adiabatic features appear: (i) time-dependence of the spectral function is seen much before than that in the effective temperature (advanced time-dependence), (ii) a big transfer of spectral weight to frequencies greater than the initial temperature occurs at an intermediate time (kink formation) and (iii) new peaks with decreasing amplitudes but in greater numbers appear even after the effective temperature has stabilised (persistent oscillations). We find four broad routes to thermalisation for lower values of spatial momenta. At higher values of spatial momenta, kink formations and persistent oscillations are suppressed, and thermalisation time decreases. The general thermalisation pattern is globally top-down, but a closer look reveals complexities.


Introduction
Holographic strongly interacting large N field theories present us a great opportunity to learn about non-equilibrium dynamics of quantum many-body systems. This is because the time-evolution of multi-point Schwinger-Keldysh correlation functions, which form the quantum Martin-Schwinger hierarchy, can be expected to be solvable exactly via classical gravitational dynamics in one higher dimension involving a few fields coupled to gravity. It is indeed hard to find examples of such potentially exactly solvable interacting nonequilibrium quantum systems, unless we impose many additional symmetries or complete integrability [1,2].
At weak coupling and large N , the study of time dependence and thermalisation of two-point correlation functions has been done in a few special examples, where thermalisation has been successfully demonstrated [3] using the two-particle-irreducible (2PI) effective action technique (see [4] for a review). Nevertheless, such methods are restricted to special initial conditions and generic couplings to the environment also cannot be easily introduced. Understanding of such generic non-equilibrium evolutions is indeed desirable for many reasons. One such reason for instance, is to obtain the quantum generalisation of the theory of fluctuations in non-equilibrium classical statistical systems (see [5,6] for reviews). Another reason is to understand if there could be a sufficiently large measure of initial conditions in special quantum many-body systems which could resist thermalisation and decoherence.
The power of holography suggests that we may be able to get some hints for answers to such questions, because both generic initial conditions and couplings to environment can be realised using present numerical gravity techniques (see [7] for a review).
By and large, the studies on non-equilibrium aspects of the holographic correspondence, have so far been focused on the time-evolution of expectation values of gauge-invariant operators, i.e. on one-point functions. This has already led to many interesting results, for instance an explicit demonstration of emergence of hydrodynamic behaviour (see [8] for a review) -which in turn has been instrumental in understanding qualitatively some tantalizing aspects of quantum critical systems and the quark-gluon plasma formed by heavy-ion collisions. Relatively less attention has been devoted to explicit calculations of time-dependence of multi-point correlation functions for generic initial conditions and couplings to environment.
We will focus here on obtaining a method for finding the time-dependence of the holographic spectral function explicitly in generic non-equilibrium states. In our concrete examples, we will find radically different qualitative behaviours of the time dependence of the spectral function which reveal characteristic features of the dynamics of the non-equilibrium state, even when the one-point function of the corresponding operator does not show any distinctive transition in its behaviour as we vary external parameters. For a large class of states with homogeneous time-dependent dynamics, we will show that the holographic spectral function can follow four different patterns of thermalisation at low (spatial) momentum. We will also examine the patterns for higher momentum.
Previously, a proposal has been given in [9] for obtaining the Schwinger-Keldysh holographic partition function, by extending the Schwinger-Keldysh contour in the bulk spacetime. This proposal can in principle lead us to obtain the Schwinger-Keldysh correlation functions in non-equilibrium states of a d−dimensional holographic quantum field theory that can be explicitly constructed by performing a d−dimensional Euclidean path integral with operator insertions. This path integral leads to an explicit (d+1)−Euclidean dual bulk spacetime via the holographic dictionary, which should then be glued to a (d+1)−Lorentzian spacetime that allows us to obtain the dynamics in the real part of the Schwinger-Keldysh contour in the dual field theory. 1 At the level of two-point Schwinger-Keldysh correlation functions, it has been argued in [13] that this prescription is equivalent to first obtaining the bulk Schwinger-Keldysh correlation functions in the dual (non-equilibrium) spacetime, and then extrapolating the latter to the boundary as suggested earlier in [14,15]. Similar arguments have also been made in [16]. The method of extrapolating non-equilibrum bulk correlation functions to the boundary has been applied in several other earlier and later works [17][18][19][20][21][22] for special non-equilibrium states. The major focus has been to understand the time-dependence of holographic correlations functions in a simple non-equilibrium state that undergoes instantaneous thermalisation, and which corresponds to a bulk dual represented by an AdS-Vaidya geometry that is formed by an ultra-thin homogeneous null shell of infalling matter (for applications to the quark-gluon plasma formed by heavy ion 1 For the thermal state, a very different realization for the holographic Schwinger-Keldysh contour has been given earlier in [10][11][12] using the causal geometry of the eternal anti-de Sitter black hole. collisions see [23][24][25][26]). It can be argued that the non-equilibrium retarded correlation function for a probe operator with a large scaling dimension can be found using the geodesic approximation in the bulk (see [13] for some not so obvious subtleties). Several works using this geodesic approximation [20,[27][28][29] have provided some interesting hints for mechanisms of thermalisation of the spectral function. The general features of thermalisation of the retarded propagator as revealed by the geodesic approximation have been pointed out particularly in [27].
However, the method of extrapolating bulk correlations to the boundary is inadequate from two points of view. Firstly, it is desirable to have a method for obtaining non-equilibrium Schwinger-Keldysh correlation functions in holography without using any additional input other than those which directly follow from the holographic correspondence itself. Secondly, it is also desirable that such a method can be used in practice to make actual computations in arbitrary non-equilibrium states. In general bulk spacetimes, one needs additional inputs for defining the bulk correlation functions themselves (as proposed in [17,18] for instance), whose general applicability far away from equilibrium is unclear. On the other hand, the bulk extension of the Schwinger-Keldysh contour as proposed in [9] is not so easy to be implemented in practice, because determining the initial data for which a Lorentzian spacetime can be smoothly continued to an Euclidean spacetime is not always straightforward.
In an earlier work [30], some of the authors have co-developed a method to obtain non-equilibrium spectral function in generic near-equilibrium states, which can be described holographically as black branes with decaying quasinormal mode fluctuations. It was shown that using the perturbative amplitude and derivative expansions developed in [31,32] one can systematically extract the time-dependence using an universal generalisation of the Son-Starinets prescription, which guarantees regular backreaction of the probe field dual to the operator under consideration. Indeed this method is exactly equivalent to the method to be proposed in the present work. 2 We will find that one can indeed extend the method of [30] even far away from equilibrium to obtain the holographic retarded propagators of arbitrary operators in generic non-equilibrium states using only the fundamental state/geometry map as the guiding principle. Numerical implementation of our method will be based on [7].
The organisation of the paper is as follows. In Section 2, we will describe our method for obtaining the holographic retarded propagator in detail and how it can be implemented in practice numerically. We will show explicitly that our method agrees with the Son-Starinets prescription at thermal equilibrium. Furthermore, we discuss the numerical implementation for a specific class of geometries that capture holographically broad features of non-equilibrium states driven by a homogeneous quench of an arbitrary duration.
In Section 3, we will present our results. Following the numerical computation, we find four broad patterns of thermalisation of the spectral function at spatial momenta lower than the initial temperature. These patterns have very novel features which cannot be found in the geodesic approximation. As the modulus of the spatial momentum (|k|) is increased we observe smooth changes in the thermalisation patterns. For instance, the thermalisation time corresponding to any quench duration decreases with increase in |k|.
In the concluding section, we discuss general lessons that can be learnt from our results. We will also formulate some broad questions, and discuss how we can address them in future works. The appendix provides some supporting results.
2 Methodology and a simple implementation

The prescription
The retarded propagator is simply related to linear response even when the system is far away from equilibrium. Let us consider a d-dimensional system perturbed from an initial state |Ψ by an external perturbation. We denote by H T the total Hamiltonian given by where H is the unperturbed Hamiltonian and ∆H is a perturbation that involves coupling of an external source f (x, t) to an operator O(x, t): with γ being a small dimensionless parameter denoting the strength of the perturbation. It is assumed that ∆H(t) can be switched on adiabatically so that it does not affect the initial state |Ψ . Then it follows that the time evolution of the expectation value of the operator O (in the Heisenberg representation with the original Hamiltonian H) is given by: where the retarded propagator G R (x, t; x , t ) is defined as: Our method for computing the holographic retarded propagator implements the computation of causal response described above using the state/geometry map of holography, which is the tenet that corresponding to every solution without naked singularity in the theory of gravity there exists a state in the dual field theory, and vice versa. The chief novelty of our method involves • implementation of adiabatic switching of the source, and • understanding of the right initial conditions to be imposed in the bulk when this additional source is switched on.
Regarding the first point, we show that we can take a suitable numerical limit where the source is a delta function in time, implying that not only the source but also all its time derivatives vanish within desired numerical accuracy at the initial time. This automatically implies right initial conditions in the bulk. The response then directly leads us to the retarded propagator. Our method thus differs from previous ones which employ calculation of the bulk-to-bulk retarded propagator and extrapolating it to the boundary (as followed in [15], [13], etc.), or establishing an analogue of Schwinger-Keldysh contour in the bulk (as proposed in [9]). Nevertheless, our method which is generally applicable to any nonequilibrium state should be equivalent to any other alternative method whenever the latter is applicable. In particular, we will be able to reproduce the retarded propagator obtained via Son-Starinets prescription [34] in thermal equilibrium.
In what follows, we will review the state/geometry map as it is implemented in numerical holography, and state our prescription precisely with full generality. The expert may like to skip the remaining part of this discussion and directly move to the next subsection, where we implement our prescription in a specific simple example, and then come back to the present subsection as the generalisation of our method could be obvious. However, a reader who is not fully familiar with the methods of characteristics as it is applied to numerical holography can benefit from the following brief review and the description of how our prescription can be implemented in arbitrary non-equilibrium states.
In order to make a precise statement of our prescription for the holographic retarded propagator, it is useful to recall the state/geometry map of holography first and how it is implemented via method of characteristics [35]. We will follow the notations of [7]. For the sake of simplicity, let us assume that we are discussing a sector of states at strong coupling and large N limit in which the energy-momentum tensor t µν and a scalar operator O are the only single-trace operators that have non-vanishing expectation values. In this case, we may restrict ourselves to studying the dynamics of classical Einstein's gravity and a minimally coupled scalar field in the (d + 1)−bulk dual.
Let X denote the (d + 1)-bulk coordinates collectively, and then we split them as a radial coordinate r, a time coordinate t and spatial field theory coordinates x i , with i = 1, · · · , d − 1. For the sake of convenience, we choose an Eddington-Finkelstein type gauge in which we can write the bulk spacetime metric in the form: The action of the Einstein-scalar theory is The scalar field Φ(X) is generically a function of all the bulk coordinates. In an asymptotically AdS spacetime, the leading asymptotic (r → 0) behaviour of the bulk metric specifies the background metric g µν on which the field theory is living: In this paper, for the sake of simplicity we will assume this to be flat Minkowski space, and therefore: g Furthermore, the asymptotic behaviour of Φ(X) gives the source for the operator, J(t, x i ), via: Above ∆ is the scaling dimension of the single-trace operator O, and this determines the mass of the dual field Φ. If J(t, x i ) is non-vanishing, it implies a Hamiltonian perturbation of the dual field theory by JO -this is often how a non-equilibrium state is realised experimentally. The asymptotic fall-off of A(X), F i (X) and H ij (X) gives us the various components of t µν (t, x i ) , i.e. the expectation value of the energy-momentum tensor in the dual state, and similarly that of Φ(X) gives us O(t, x i ) , i.e. the expectation value of the dual scalar operator. The latter can be extracted from the gravitational solution using the standard holographic dictionary. In order to specify the dual state, we furthermore require to provide initial conditions in addition to the sources (g µν and J) so that we can determine the gravitational solution uniquely. These are given by the following functions at the initial time hypersurface t = t in (see [7] for example): 3 For appropriate initial conditions, the solutions of classical gravity will not have naked singularities in the bulk. Such solutions of classical gravity will correspond to states in the dual field theory. In case of most such initial conditions, the exact analytic solutions cannot be readily obtained. Nevertheless, one can numerically obtain them [7]. In order to check if a given initial condition is admissible or not, meaning that whether it does or does not lead to naked singularities, we will require to obtain the corresponding numerical solution over a suitable computational domain (i.e. for 0 < r < r c , with r c being a suitable cut-off value of the radial coordinate) and then check if a regular apparent horizon lies within this computational domain (i.e. within 0 < r < r c ). 4 Let us now consider how we should obtain the retarded propagator of an operator in a given non-equilibrium state that corresponds to a solution in gravity with specific initial data (2.10). To begin with, let us consider the case when the concerned operators are t µν and the scalar operator O, i.e. those operators which have non-vanishing expectation values 3 Note that the (d + 1)−variables A(X), Fi(X) and det(H) can be regarded as auxiliary variables -their complete radial profile on the initial time hypersurface are also generated from the initial data. Also the initial radial profile for d+Φ ≡ (∂t + A(X)∂r)Φ is generated by the Klein Gordon equation from that of Φ(X). 4 The apparent horizon always lies behind (i.e. farther from the boundary than) the event horizon coinciding with the latter only at late time (see Figs. 1(a) and 1(b) for the case of AdS-Vaidya geometries). Therefore, finding the apparent horizon within the computational domain ensures absence of naked singularities.
in the dual state. Obviously, we need to then perturb the field theory with additional weak sources: 5 Holographically, this implies that the leading asymptotic behaviour on the gravity side should be modified as follows: and furthermore lim r→0 Our crucial observation, which we will confirm numerically in the next subsection is that we can indeed implement adiabatic switching of the sources above. 6 This means that at initial time t = t in , the additional sources h µν (t, x i ) and f (t, x i ) and their time derivatives to all orders can be adjusted to be arbitrarily small with desired numerical accuracy. As a result, the initial data (2.10) can be kept unchanged so that they will lead to gravitational solutions without naked singularities, even in the presence of the additional sources.
Note that, it is not sufficient to make h µν (t, x i ) and f (t, x i ) vanish at t = t in , we also need to ensure their time-derivatives to all orders vanish too with desired accuracy. Recall that the Fefferman-Graham expansion makes the higher order time derivatives of sources appear in the radial profile of the bulk fields on-shell (on the equal-time hypersurface). Therefore, it can be expected that only if the additional sources and their time derivatives vanish at initial time, we can be sure about the original initial conditions to still remain admissible even though sources do get turned on in the future. The latter simply follows from causality, because the additional sources, which are adiabatically switched on in the future, cannot affect the initial conditions. We will show in the following subsection that the above adiabatic turning on of sources for switching on the Hamiltonian perturbation ∆H(t) can be achieved in holography using a sufficiently narrow Gaussian time-profile (for h µν and f in the specific case under current discussion). Furthermore, we will also show that when the sufficiently narrow Gaussian time-profile is appropriately normalised, it mimics a delta function source to any desired degree of accuracy for computing the response at any arbitrary time. This will simplify the numerical computation of the retarded propagator.
In order to extract the retarded propagator, we simply compute (numerically) the new gravitational asymptotically AdS solution with same initial conditions (2.10), and sources 5 We assume that there is no matter anomaly in the ground state (i.e. no contribution to conformal anomaly coming from non-trivial source J(t, x i ) coupling to the scalar operator) so that tr(t) = t µ µ = 0otherwise there will be an additional contribution to ∆H(t) proportional to tr(H). 6 This adiabatic switching of the additional sources f and hµν appearing in ∆H(t) should not be confused with the (non-)adiabaticity of the non-equilibrium state in which the retarded propagator is to be evaluated. The latter is decided by the external source J in (2.13) which is part of the Hamiltonian H itself -J can be designed to drive instantaneous/slow quench dynamics as we will discuss in the next subsection with an example.
(2.12) and (2.13) where the additional pieces h µν and f are adiabatically switched on as mentioned above. From the new solution, we can extract the new expectation values t µν (t, x i ) γ 1 ,γ 2 and O(t, x i ) γ 1 ,γ 2 using the standard holographic dictionary. Then it follows that: where G R are the following retarded propagators in the dual field theory: Since we know t µν γ 1 ,γ 2 and O γ 1 ,γ 2 , and have already specified h µν and f , we can readily extract the dual retarded propagators in the small γ 1 and γ 2 limit. The explicit numerical procedure will be discussed in the following subsection. It is quite clear that going to higher orders in γ 1 and γ 2 , we can extract fully causal 3-point and other multi-point correlation functions. Equivalently, we can also solve for the linearised fluctuations of the bulk fields δG M N and δΦ with initial conditions: at t = t in , and boundary sources: about the background gravitational solution given by initial conditions (2.10) in the presence of the sources (2.8) and (2.9). In that case, we can readily extract the dual retarded propagators using (2.14) by setting γ 1 = γ 2 = 1 (which will be admissible because we are restricted now to the linear limit). The detailed numerical method will be presented in the following subsection. The above linearisation is certainly good enough for extracting the causal 2-point correlation function, i.e. the retarded propagator.
Clearly we can apply the same procedure for obtaining the retarded propagator of an operatorÕ that has vanishing expectation value in the non-equilibrium state dual to the gravitational solution obtained from the initial conditions (2.10) in the presence of the sources (2.8) and (2.9). We simply study the linearised equation of motion for the dual bulk fieldΦ with initial conditionΦ (r, that is switched on adiabatically. From the solution we extract Õ γ and then the retarded propagator GÕÕ R (x, t; x , t ) as discussed before. Although we will restrict ourselves to the retarded propagator, it is to be noted that the fully causal multi-point correlation functions involving operatorsÕ which have vanishing expectation values and those which do not can also be extracted from the full non-linear solution of the gravitational equations. In this case we switch on all sources adiabatically and keep all initial conditions unaltered. The latter implies that the initial conditions for the bulk fields that vanish in the original solution should be kept trivial, and those for the other bulk fields are kept the same as in the original solution dual to the state in which we are evaluating the correlation functions. We will call operatorsÕ which have vanishing expectation values probe operators from now on.
The other advantage of our method is that we can implement it even when the gravitational solutions are not known analytically, but numerically. In the following subsection, we will discuss a simple example involving a class of background gravitational solutions which are time-dependent and can be known analytically. The retarded correlation function of a probe scalar operator will be found to exhibit diverse qualitative behaviour that will reveal the characteristic features of dynamics of the bulk spacetime.

A simple example
Non-equilibrium states which can be easily realised in experiments are those which result from an external perturbation that is switched on and off adiabatically, resulting in the system being driven out of equilibrium and then eventually settling down to a new equilibrium. Gravitational solutions dual to such non-equilibrium states can be also obtained numerically. In this work, we will study a simple class of gravitational backgrounds which give a rough approximation to such dual non-equilibrium states [36]. In the concluding section, we will point out that indeed our results here are also necessary to develop intuition about the time-dependence of the retarded propagator in the exact non-equilibrium background, to which we will turn our attention to in the future. We will focus on 2 + 1-dimensional holographic quantum many-body systems, and therefore 3 + 1-dimensional bulks.
The simple gravitational backgrounds on which we will focus our attention here are given by the metric 7 where A(r, t) has the form of the standard black brane blackening function with a time-dependent mass function This corresponds to the case when an ultra-thin homogeneous null shell infalls into a black hole changing its mass from M in to M f . From the dual field theory point of view, thus the limit α → 0 corresponds to a very rapid quench. The time-dependence of the retarded propagator in such backgrounds has been studied in the literature mostly for the case M in = 0 (see the discussion in the introduction). When α is finite, the gravitational background (2.20) is then generated by an infalling homogeneous null shell of effective width 2α. The limit α → ∞ corresponds to an adiabatic process of transition of the black hole mass. We call α the quenching parameter. We recall that in the static background where M (t) = M , the horizon of the black brane is located at r H = l 2/3 M −1/3 . The event and apparent horizons typically do not coincide in non-stationary gravitational backgrounds. In the AdS-Vaidya spacetimes (2.20), the location of the apparent horizon r AH (t) is given by the solution of A(r AH (t), t) = 0. The location of the event horizon r EH (t) is obtained by solving the null geodesic equatioṅ which ensures that in the far future it reproduces the event horizon of the final black hole. In Fig.1 we show the evolution of the apparent and the event horizon for various values of α. In these plots, and for the remaining paper we set units l = 1 and choose M in = 1, M f = 8 (i.e. ∆M = 7). The final temperature T f is then twice higher than the initial temperature T in . Note that both in the far past and far future, the event and apparent horizons of (2.20) coincide.  The Hawking temperature of the static black brane is given by We may thus think about (3M (t) 1/3 )/(4πl 2/3 ) as an effective instantaneous temperature T (t) of the AdS-Vaidya spacetime. However, such a description is strictly valid only in the adiabatic α → ∞ limit.
We are going to study the time-dependence of the retarded propagator of a probe scalar operator O in the non-equilibrium state dual to the AdS-Vaidya gravitational background (2.20) using the prescription described in the previous subsection. This operator O does not have a spontaneous expectation value in the dual state, and is dual to a bulk scalar field Φ that trivially vanishes in the background (2.20) before turning on the boundary source perturbation. We will also consider the case when this operator is relevant with scaling dimension ∆ = 2. The latter translates to the mass of the dual bulk scalar field being given by m 2 l 2 = −2. 8 We also assume that the bulk scalar field Φ is minimally coupled to gravity and does not couple to the matter composing the null shell that generates the AdS-Vaidya background (2.20). The bulk action for the scalar field is then Since the background (2.20) preserves spatial homogeneity, the solution for the Klein-Gordon equation for Φ in this background can be written in the form The Klein-Gordon equation then reduces to where ≡ ∂ r , and is the derivative along the outgoing null geodesic. Any solution of (2.27) can be written asymptotically in the form Applying d + to the equation above, we readily obtain the asymptotic form of d + Φ 0 as below: where˙≡ ∂ t . According to the holographic dictionary, the external source f (t, x) = f 0 (t, k)e ik·x leads to the perturbation of the Hamiltonian of the system by In practice, the Klein-Gordon equation (2.27) can be solved using the method of characteristics. At initial time t = t in , we specify Φ 0 (r, t, k). We then radially integrate (2.27) to generate d + Φ 0 (r, t, k) at t = t in using the boundary condition d + Φ 0 (r = 0, t in , k) = −f 0 (k, t in )/2 that follows from the asymptotic expansion (2.30). Then we obtain ∂ t Φ 0 from d + Φ 0 using ∂ t Φ 0 = d + Φ 0 + (A/2)Φ 0 (recall the definition (2.28) of d + ), and update Φ 0 (r, t, k) at the next time step t = t in + ∆t. We further continue at the next time step by performing the radial integration in (2.27) to obtain d + Φ 0 (r, t in +∆t, k) using the boundary condition provided by the given source f 0 (k, t in + ∆t), and so on. Thus, the unique solution for Φ 0 can be obtained numerically for a given radial profile at initial time and a boundary source f 0 specified at all times. Here, we will use the spectral method with 30 grid points to perform the radial integration for obtaining d + Φ 0 (r, t, k) at each time step. In order to update in time, we will use a fourth order Adams-Bashforth time stepper with ∆t = .0005.
As discussed in the previous subsection, in order to obtain the retarded propagator G OO R (x, t, x , t ) in the non-equilibrium state of interest, we need to switch on the source 9 We obtain the result below by changing to Fefferman-Graham coordinates (z,t, x, y) first using After doing this change of coordinates, the asymptotic expansion in the Fefferman-Graham coordinates Φ0(z,t, k) is written in terms of the coefficients in the Eddington-Finkelstein coordinates (2.29) as The one-point function O f is minus one times the coefficient of z ∆ of the Fefferman-Graham expansion times (2∆ − d) when log r terms are absent (see [38] for instance). Since in our case ∆ = 2 and d = 3, we obtain (2.34).
f (t, x) adiabatically by designing an appropriate profile for it, then obtain O(t) f by solving (2.27) in the AdS-Vaidya background (2.20). We need to impose the initial condition Φ(r, x i ) = 0 in the far past t = −∞, and compute the retarded propagator from the relation: Note that due to spatial homogeneity of the background (2.20), the retarded propagator can depend on x − x , and not x and x separately. It is therefore convenient to study the Fourier transform of G OO R (x, t; x , t ) with respect to x−x , which we denote as G OO R (t, t ; k) (with a slight abuse of notation). Clearly, Remarkably, we find that it is possible to implement the delta function limit of the source, so that f 0 (t , k) = δ(t − t 0 ). Practically, we implement this delta function as a narrow (σ → 0) limit of the Gaussian function 10 which is centred at t = t 0 and is normalised such that We will soon show that numerically it is indeed possible to take the limit σ → 0 smoothly in (2.37) for evaluating O 0 (t, k) f (and we will also make error estimates when σ is finite). Indeed if f 0 (t ) is given by δ(t − t 0 ), then we can readily read off G OO R (t, t ) from (2.36) using Although the integral over t runs from −∞ to t, we can extend this integral over the entire real line, because O 0 (t, k) f and hence G OO R (t, t 0 ; k) vanishes for t > t by causality. Thus the above result indeed follows if f 0 is given by δ(t − t 0 ). Therefore, we can directly read off the retarded propagator G OO R (t, t 0 ; k) from O 0 (t, k) f (given as discussed earlier by (2.34)) for a fixed value of t 0 . We need to vary t 0 though to obtain complete information about G OO R (t, t 0 ; k). Using (2.34), we can also write (2.39) as where f 1 is evaluated in presence of the source f 0 (t, k) = δ(t − t ) and with initial condition Φ 0 (r, k, t in ) vanishing in the far past. 10 Boundary source perturbation by a narrow but finite width Gaussian profile has been often used to drive the bulk scalar in numerical time evolution and compute the one point function [39][40][41][42]. Scalar quenches have also been studied with other boundary source profiles [43][44][45][46][47][48].
There is another advantage of using the narrow Gaussian approximation for the Dirac delta function source (2.37), namely that it readily implements adiabatic switching of the source that is necessary for our prescription for obtaining the retarded propagator discussed in the previous subsection. The time t in when the initial condition Φ(r, x i ) = 0 is set can be chosen well before t = t 0 (i.e. t in t 0 − σ) when the narrow Gaussian profile (2.37) of the source has a significant non-trivial value -the source and its time derivatives to all orders vanish to desired level of accuracy at t = t in . Furthermore, as we will show quantitatively below, indeed for the causal response O 0 (t, k) f for t > t 0 (when G OO R (t, t 0 ; k) should be non-vanishing), the convergence to this delta-function limit σ → 0 is very rapid. Thus we can readily obtain G OO R (t, t 0 ; k) using (2.40) for any given t 0 by appropriately solving the bulk Klein-Gordon equation for the bulk scalar field Φ with initial condition Φ(r, t in , x i ) = 0 at an appropriate value of the initial time t in .
To obtain better physical intuition about G OO R (t, t ; k) and also to make connections with experiments, it is instructive to study the Wigner transform of G OO R (t, t ; k) which is defined by below: The Wigner transform is thus a Fourier transform with respect to the relative coordinate t rel = t − t and then after this transformation G OO R depends on the frequency ω and the average time t av = (t + t )/2. Again for the sake of notational simplicity, we denote the Wigner transformed G OO R as G OO R (ω, t av , k). When t av → ∓∞, G OO R (ω, t av , k) does not depend on t av , and coincides with the thermal retarded propagator G OO R (ω, k) with temperatures (3M 1/3 in,f )/(4πl 2/3 ) corresponding to the black brane backgrounds with masses M in and M f respectively, as discussed before.
We will focus mostly on the spectral function ρ OO (ω, t av , k) which is experimentally measurable for example in solids using angle resolved photo electron spectroscopy (ARPES) and is defined as Causality implies that G OO R (ω, t av , k) is analytic in ω in the upper half complex plane, therefore the Kramers-Kronig relation is valid, with P denoting the Cauchy principal value. Thus causality implies that the spectral function encodes complete information of the retarded propagator. Because G OO R (t, t ; k) is real, we also note that, where G A is the advanced propagator. This implies: It is therefore sufficient to restrict ourselves to the region ω > 0. 11 Since we are studying the case of a Hermitian operator dual to a real scalar field, G OO R (t, t ; k) is real as well. Therefore, it follows from (2.40) that where f 1 is obtained from solving the Klein-Gordon equation in the background (2.20) in the presence of the source f 0 = δ(t − t ) and with initial condition Φ 0 (r, t in , k) = 0 in the far past.
In order to get a first look at how well we can approximate the delta function source by a narrow Gaussian, we study the behaviour of f 1 (t, k). We set k = 0 and locate a narrow Gaussian perturbation at t 0 = −9. 12 In Fig. 2(a) and 2(b), we show the result of f 1 for σ = 0.01 when α = 0.05, k = 0. When the shell is thin, the background immediately approaches the thermal equilibrium, and hence a sharp transition from the initial to final temperatures is observed in Fig. 2(b). The oscillating exponential decay of the one point function away from the shell is controlled by the lowest lying quasi normal mode frequency. For ∆ = 2, k = 0 in AdS 4 , the value is 3ω/(4πT f ) = 1.228 − 1.533i as can be obtained using the method of [49]. Our numerical results agree with this value.
By varying σ, we find that the behaviour of f 1 (t, k) converges as we decrease σ. See Fig. 2(c) comparing σ = 0.1, 0.01, 0.001. At early time, different curves of f 1 do not coincide due to different widths of the corresponding sources, but at late times the ring down parts lie almost on top of each other. However, for sufficiently small σ, the curves of f 1 coincide well at all times. We can take the limit of delta function source by computing for several small σ and extrapolating to σ → 0. Using a small σ is sufficient for practical numerical purpose.
One may easily estimate the influence of the finite width σ of the Gaussian profile for the source on the numerical error in finding the Wigner transform. Note that the Fourier transform of the Gaussian profile (2.37) for f 0 (t, k) is: (2.48) Following the above it is not hard to see that when we extract spectral function treating the source f 0 as a delta function in real time using the Fourier transform (2.47), we are making significant errors only for very high values of ω, i.e. for ω σ −1 . However, even before reaching such high values of ω, we will find that ρ(ω, t av ; k) asymptotes the vacuum form (with ∆ = 2): (2.49) 11 The reader can readily recall that for a free bosonic field χ with mass m, the vacuum spectral function is given by: with k = √ k 2 + m 2 . 12 Throughout this paper, we choose Min = 1 and M f = 8 for the AdS-Vaidya background geometry as mentioned before. for all values of t av when k is not too large. This also requires that we should take the initial and final black hole masses M in,f σ −1 , because the mass parameters control the frequency above which time-dependent spectral function asymptotes to the vacuum form (2.49), at least for low values of |k|.
Note in a conformal field theory there is no intrinsic scale, so only ratios of various mass parameters, time scales, etc. affect physical results. Therefore, without loss of generality we can choose M in = 1. This however prevents us from studying the limit when we start from zero temperature, i.e. pure AdS background in the far past. Furthermore, the approximation of the delta function by a narrow Gaussian implies that we cannot study the case when M f is too large as we will be interested in the region ω/M f ≈ 1 too. Despite these limitations, we are equipped to explore a diverse range of qualitatively different behaviour of the time-dependence of ρ(ω, t av , k).

The thermal spectral function
Perhaps the simplest test of our prescription is to reproduce the well-known holographic thermal spectral function. Due to time translation invariance, ρ(ω, t av , k) does not depend on t av in thermal equilibrium. The dual gravitational background is simply the background (2.20) with M (t) = M 0 , a constant. Thus this is just a special case of the discussion in the previous subsection.
It is useful to study where ρ AdS = ρ T =0 is given by (2.49), in order to extract the purely thermal contribution. Furthermore, it is also useful to studyρ since 2ω is the ubiquitous state-independent contact term as evident from (2.47). Note that this term cannot be removed by counterterms. When k = 0, ρ AdS = 2ω. Therefore, ρ sbt =ρ when k = 0. We use σ < 0.01 in our computations.  In Fig. 3, the thermal spectral function for k = 0 obtained using our method is shown. In the plots, we define a dimensionless quantitiesX from its dimensionful counterpart X asX ≡ (3/4πT H ) n X with n being the mass dimension of X. In Fig. 3(a), we plot the thermal spectral function, and shown in Fig. 3(b) is ρ sbt in which the linear growth in the large ω region gets subtracted away. Similarly, in Fig. 4 results fork = 4.5 are shown. From Fig. 4(a), we find the spectrum is exponentially suppressed in ω < k. In Fig. 4(b), the function in ω < k is the same as that in Fig. 4(a), but there is a sharp peak at ω = |k| because of the subtraction. The spectral function in ω > k has the qualitatively same distribution as in Fig. 3(a). Besides ρ and ρ sbt , we also plotρ in Fig. 4(c) which is determined by f 1 contribution as in (2.47). For small ω, the spectral function depends linearly on ω because of the 2ω subtraction.
Using the well known method by Son and Starinets [34] for obtaining the holographic thermal spectral function (see [50,51] for examples), we have checked that the plots in Figs. 3 and 4 and similarly those for other values of |k| are reproduced. However, actual procedures in calculations are quite different, 13 and our method thus surely passes the simplest test of reproducing the holographic thermal spectral function.

Results and a classification of routes to thermalisation
In the previous section we have discussed a general prescription for finding the timedependent retarded propagator using holography in generic non-equilibrium states of a strongly coupled large N theory, and have also shown that it reproduces the well-known thermal case. Furthermore, we have discussed how we can adapt our method to the specific AdS-Vaidya gravitational backgrounds ((2.20) and (2.22)) that holographically capture features of spatially homogeneous non-equilibrium dynamics driven by a quench to a certain degree of approximation. Here we will present our numerical results and classify various routes to thermalisation of the time-dependent spectral function by varying the time dura- 13 The main difference between our method and that of Son and Starinets is that we solve an initial value problem where we impose regularity at the horizon via the choice of initial conditions, whereas in the latter case a boundary value problem is solved with the infalling boundary condition at the horizon. tion of the quench. We will find various qualitative patterns of thermalisation of the spectral function as we interpolate between adiabatic and instantaneous quench limits (these limits are described below). Our results will be specifically for a scalar operator with ∆ = 2 in a (2 + 1)−dimensional holographic CFT -in this case we will obtain exact results going beyond the geodesic approximation examined extensively in the literature.
The quenching parameter α in (2.22) labeling the (3 + 1)−AdS-Vaidya gravitational backgrounds gives us the effective transition time from the initial to the final thermal equilibrium represented by AdS-Schwarzschild black holes of masses M in = 1 and M f = 8 respectively (note M in = 1 also sets our units for field-theory observables. From now on we remove hats from the dimensionless quantities in the plots). The Bondi mass function M (t) of these geometries are plotted in Fig. 5 for various values of α. Clearly α → 0 is the instantaneous quench limit and α → ∞ is the adiabatic limit. The mid-point of the transition is always chosen to be t = 0. We may readily identify a time t α i for a given α such that for t < t α i , the mass function M (t) (or the effective instantaneous temperature) deviates from M in (the initial temperature) by less than 1 percent, and similarly a time t α f such that for t > t α f the mass function deviates from M f (the final temperature) by less than 1 percent.
In Section 2.2, we have illustrated our method for computing G OO R (t, t 0 ; k) (from now on we drop the superscripts OO). We present the three dimensional plots for |G R (t, t 0 ; k = 0)| in case of rapid quench (α = 0.05) and in case of slow quench (α = 2.0) in Figs. 6(a) and 6(b) respectively. The contact termδ(t − t 0 ) in (2.40) has been subtracted in these plots, which therefore simply show |f 1 (t, k = 0)| that results from causal response to the source f 0 (t , k = 0) = δ(t − t 0 ). Clearly, we see that G R vanishes for t < t 0 as expected from causality. Across t 0 = 0, when the background transits from a black brane with temperature T in to that with temperature T f , the decay and oscillation rates of |G R | change by a factor of T f /T in = 2 as we should expect from quasi-normal mode analysis. Also the wake up of |G R | at t = t 0 gets enhanced by a factor of (T f /T in ) 2 = 4 as we move from t 0 = −∞ to t 0 = ∞. These transitions of decay and oscillation rates, and also wake up amplitudes are rapid in the case of α = 0.05, but slow when α = 2.0. Our strategy for finding ρ(ω, t av , k) has also been discussed in detail in Section 2.2. We will actually study the subtracted spectral functioñ ρ(ω, t av , k) = ρ(ω, t av , k) − 2ω, (3.1) with the state-independent contact term 2ω removed. In order to evaluateρ(ω, t av , k) therefore, we just need to perform the integral in (2.47) for various values of t av . This is done in two steps. First, for a fixed t av we obtain the response f 1 (t, k) for the source effectively f 0 (t ) = δ(t − t 0 ) with t 0 = 2t av − t. Both t and t 0 are varied keeping t av fixed. Note causality implies that f 1 (t, k) vanishes for t < t 0 . Furthermore, if t and t 0 are far apart, i.e. if t rel = t − t 0 is large, f 1 (t, k) decays exponentially at a rate controlled by the final black hole temperature -so the integral does not receive significant contribution for large values of t rel . In the second step, the variable t rel = t − t 0 is replaced by the variable t so that the integration in (2.47) can be performed numerically by discrete summation over f 1 (t, k) for various values of t (we use step size of 0.01). We thus obtainρ(ω, t av , k) as a function of ω for various values of t av at a fixed value of |k|. Clearly for any value of α, the spectral function ρ(ω, t av , k) will coincide with that of the thermal form corresponding to the initial black hole as t av → −∞ and with that of the thermal form corresponding to the final black hole as t av → ∞. Nevertheless, the effective transition time will depend on α. We denote the time t α av(i) as that for which ρ(ω, t av , k) becomes independent of t av when t av < t α av(i) and coincides with the thermal form of the   initial black hole up to 1 percent numerical difference 14 (to be precise we mean that the maximal deviation is less than 1 percent). Similarly we define t α av(f) such that for t > t α av(f) , ρ(ω, t av , k) becomes independent of t av coinciding with the thermal form of the final black hole up to 1 percent numerical difference (to be precise we mean that the maximal deviation is less than 1 percent again). Interestingly, t α av(i) and t α av(f) need not coincide with t α i and t α f respectively, where the mass function M (t) coincides with the initial and final values up to 1 percent numerical difference. Indeed, this turns out to be true specially when α is small, α ≤ 0.1 and very large α ≥ 2 as shown below. We divide our study in three broad categories, namely, |k| < M in , |k| ∼ M in and |k| > M in . The first category is represented by the choice |k| = 0, the second one by |k| = 1 and the third by |k| = 4.5. The regime |k| < M in for M in > 0 particularly shows qualitatively new patterns of thermalisation distinct from what has been found in the literature before.
Tables 1, 2 and 3 summarise the various routes of thermalisation at a glance for appropriately chosen values of α and |k| by specifying the presence/absence or degrees of  appearance of various possible features that are going to be described below. We find distinct patterns of thermalisation developing as we increase α and |k|. Each pattern will be discussed at length in the following subsections. It is to be noted that the transitions between these patterns as we increase α and |k| are not very sharp. We have also not attempted to specify the boundaries of α at fixed |k| where each pattern is observed even approximately, however we do find smooth transitions from one pattern to another as we vary α and |k|.

Numerical Results
In what follows, we will present our numerical plots for the time evolution of the spectral function starting with the lowest momentum |k| = 0 and going all the way to |k| = 4.5. In order to make the dynamics of thermalisation visible from the figures, we have chosen the following colour codes: green when the spectral function coincides with the thermal form of the initial black hole at t av = t av(i) , red when it coincides with the thermal form of the final black hole at t av = t av(f) and the other colours are used in the window t av(i) ≤ t av ≤ t av(f) . We will consider three different values of |k| representing three different cases, namely |k| < M in , |k| ∼ M in and |k| > M in respectively.
3.1.1 |k| = 0 representing |k| < M in As shown in Table 1, this regime is the richest with four distinct routes of thermalisation each corresponding to a different value of α, namely 0.05, 0.2, 0.5 and 2. We study these cases so that we we can interpolate from instantaneous quench to adiabatic transition in the dual non-equilibrium state. Although we examine them for |k| = 0, these also hold when |k| M in , as for instance when |k| = 0.5. At early and late times, the spectral function has peaks corresponding to the thermal quasinormal mode of the initial/final black holes. The location and width of these peaks are proportional to the temperature. In all the four patterns of thermalisation, the width of the spectral peak broadens uniformly reflecting rise in the effective temperature (note the mass function M (t) increases monotonically). Nevertheless, there can be various forms of distortions.
(a) α = 0.05: The most dramatic route to thermalisation is obtained for α = 0.05. The plots are in Fig 7. There are three distinct features.

1.
Advanced time-dependence: The spectral function starts getting time-dependent much before t av = t α i , which is when the mass function (effective temperature) picks up a significant change -i.e. t α av(i) < t α i . This is not in contradiction with causality because the spectral function is non-local in time by definition. However, this advanced time-dependence is pronounced only for low values of α. Here t α av(i) = −1.0 whereas t α i = −0.17. A similar advanced time-dependence has been observed in [22] in the fast quench limit. It is worth noting at this point that for large values of α corresponding to the adiabatic limit, we will find the opposite behaviour, namely delayed acceleration, i.e. t α av(i) > t α i .

2.
Kinkiness: Around the time t av = 0, the spectral peak starts getting distorted for ω > M in with significant spectral weight transferred to the higher frequencies. At t av = 0.18, 0.22 and 0.26, we can see a sharp kink in the spectral function in the high frequency region. This feature of large spectral weight transfer to higher frequencies ω > M in will also be present for somewhat higher values of α. More precisely, the spectral weight transfer to higher frequencies is defined here in comparison with the form of the thermalised spectral function at the instantaneous effective temperature.

Persistent oscillations:
For higher values of t av , the kinks disappear, but smaller peaks start appearing and disappearing thus forming transient ripples before the spectral function settles down to the final thermal form at t av(f) . As visible in Fig. 7(c), the ripples propagate towards lower frequencies. The amplitude of these ripples decrease with time, but the number of visible peaks increase (indeed zooming inside each peak we find secondary peaks as well -this is however not shown in the plots). 15 As is clear, these oscillations persist long after t av = t α f so that t α Note that the height of the spectral function peak first decreases as we increase the time from t av(i) . After some time, t av ≈ 0.09, the peak height starts increasing with time until the spectral function attains its final form.
We also did the numerics for α = 0.1 and the qualitative features remain exactly the same as that for α = 0.05.
It is to be noted though that the advanced time dependence and persistent oscillations are not mere artefacts of the Fourier transform involved in the definition of the Wigner transform of the commutator. In fact, these will be absent as we increase the duration of the quench, for larger values of which advanced time-dependence will be replaced by delayed time-dependence.
(b) α = 0.2: The pattern of thermalisation found for α = 0.2 is depicted in Fig 8. Here, we find that the advanced time-dependence is suppressed -in this case t α i ≈ −0.66 and t α av(i) ≈ −0.8. Around t av = 0, the spectral function develops kinks involving the 15 We speculate that a fractal structure may form at lower values of α. This may be related to the nonanalyticity at the onset of thermalisation found in [27] with the geodesic approximation. In order to explore lower values of α, we need to decrease our time-step. This translates to increase in computational time. In future work we will investigate if it is indeed feasible to take the instantaneous quench limit continuously within our approach.
transfer of spectral weight towards higher frequency modes ω M in . Also note that for t av > t av(i) , the height of the spectral peak initially decreases, but then it starts increasing monotonically again as the spectral function assumes the final thermal form. The spectral function still takes slightly longer time to thermalise than the background mass function M (t), i.e., t α av(f) > t α f -but the relative difference between the two is much smaller than that in the previous case. (c) α = 0.5: Going to a slightly longer duration of the quench by choosing α = 0.5, we find that the thermalisation remains non-adiabatic but the qualitative features are very different from the previous two cases as shown in Fig. 9. There is no significant advanced time dependence. However, we find a very tiny transfer of spectral weight around t av = 0 and we call such a thermalisation as mildly non-adiabatic because there are only mild distortions in comparison with the form of the thermalised spectral functions at the instantaneous effective temperatures. As in the previous case, the peak height decreases first a bit before it starts increasing monotonically again. The time scale of thermalisation of the spectral function to the final form is almost same as that for the mass function, i.e. t α av(f) ≈ t α f . (d) α = 2.0: From the trends discussed above, it can be expected that tuning α to even higher values will lead to an adiabatic thermalisation. This indeed shows up in Fig.  10 for α = 2.0. We find that both the height and width of the spectral peak monotonically increases. Due to the slow rate of change of background mass function, the spectral function approximately assumes the thermal form at the instantaneous effective temperature. However, the notion of instantaneous effective temperature is strictly not valid here as the mass is a function of time. Therefore, adiabatic thermalisation more precisely stands for uniform deformation of the spectral function without distortion.
There are two distinct features worth noting nevertheless: 1. Delayed time-dependence: Significant time dependence (by 1 percent) in the spectral function commences with a slight delay from that for the mass function M (t).
Comparing Fig. 10 with Fig. 5, one can clearly see that indeed t α av(i) (≈ −5) > t α i (≈ −6.6). This feature has been found in the geodesic approximation before [27]. 2. Slightly premature thermalisation: Furthermore, the final thermal form of the spectral function is attained (with 1 percent maximal deviation) slightly before the time-dependence of the mass function becomes insignificant, i.e. t α av(f) (≈ 5) < t α f (≈ 6.6), as also evident from Figs. 10 and 5.
We find that the qualitative features of spectral function thermalisation for |k| = 0.5 are the same as that for |k| = 0 for the same values of α. We present the plots in Appendix A for completeness.  Table 2 summarises the patterns of thermalisation for |k| = 1 = M in for the same set of values for α as in the previous case. We find lesser spectral weight transfer to UV modes as compared to the transfer in the cases for |k| = 0 and |k| = 0.5 at lower values of α. Even though some distortions are noted in the curves, kink formations and persistent oscillations are relatively suppressed. The plots in Fig. 11 for α = 0.05 are almost similar to those in Fig. 7 nevertheless. However, we find that the effective time for thermalisation of the spectral function, i.e. At this point, we do see that the thermalisation of the spectral function has a very distinct behaviour compared to that of the one point function. It is known that the imaginary part of the quasinormal-mode frequency decreases as we increase |k| (see [49]), indicating that the one-point function takes longer to thermalise as we increase |k|. We comment on this behaviour in the next subsection.
As we increase the quenching parameter to α = 0.2 and 0.5, we find that mildly non-adiabatic features show up with the absence of prominent kink formation as shown in Figs. 12, 13. Similar to the zero momentum thermalisation patterns, we find the adiabatic (distortion-free) route at α = 2, this is shown in Fig. 14. However, note that in the adiabatic limit the height of the spectral peak initially decreases as we go above t av(i) unlike the zero momentum case. Afterwards, the height increases monotonically. This is actually consistent with the behaviour of the thermal spectral functions at the instantaneous effective temperatures. Both delayed time-dependence and slightly premature thermalisation appear in the adiabatic limit.  Table 3. As observed in the previous case, adding momentum suppresses kink formation in the spectral function. We now also find that it suppresses late time persistent oscillations at α = 0.05. The pattern is mildly non-adiabatic for α = 0.05 (see Fig. 15) and adiabatic for higher values (see Figs. 16, 17 and 18). Advanced time-dependence still appears at α = 0.05. Once again, by adiabatic we actually mean distortion-free, i.e. the spectral functions at intermediate times take almost thermal forms at instantaneous effective temperatures. Since, we cannot explore lower values of α with sufficient numerical accuracy, we cannot assert surely that the persistent oscillations never occur even in the instantaneous quench limit. Nevertheless, we can expect that at higher values of |k| we only have mild nonadiabatic behaviour for very small values of α and the adiabatic pattern for higher values of α.
It is to be noted that in the case of |k| = 4.5 the height of the spectral peak always decreases. This is consistent with the behaviour of the thermal spectral functions at the instantaneous effective temperatures.

Top-down or bottom-up?
It is interesting to ask the question whether the patterns of thermalisation can be classified as top-down or bottom-up. At a fixed value of |k| and α, this amounts to asking whether higher or lower frequency modes thermalise faster. In the previous subsection, our operational definition of the time of thermalisation t α av(f) has been such that the maximal deviation of ρ(ω, t av , k) from its final equilibrium value ρ f (ω, k) should be less than 1 percent for t av > t α av(f) . This definition is blind to the ω−dependence of ρ(ω, t av , k). In order to understand ω dependence of thermalisation time, we take a closer look at ρ(ω, t av , k) at late times, i.e. when t av ≈ t α av(f) . We investigate plots of |ρ(ω, t av , k) − ρ f (ω, k)| as a function of ω in log scale for various values of t av in the instantaneous quench limit α = 0.05 for |k| = 0 (see Fig. 19(a)), |k| = 1 (see Fig. 19(b)) and |k| = 4.5 (see Fig. 19(c)). In case of α = 2 (see Fig. 20), we have plotted the same quantity for |k| = 0.
For α = 0.05, the pattern is globally top-down, but a closer look reveals complexities. There is uneven distribution of spectral weight for frequencies 0 < ω < 10. This is readily visible from the plots for t av = 0.7, 1.0 in Figs. 19(a), 19(b) and 19(c). However, the overall spectral weight seems to be decreasing slowly (with a very mild slope in the log plot which is −0.17 for |k| = 0 and 1, and −0.21 for |k| = 4.5) at higher frequencies ω > 10.  Such log-plots of the Wigner transformed spectral function are not readily available in the literature. 16 For extremely high frequencies, the overall spectral weight decays surely because in this region the spectral function always remains time-independent and is as in the vacuum. When α = 2, i.e. when we are in the adiabatic limit, the oscillations are suppressed and there is more rapid decay of the spectral weight at higher frequencies.
In this case, we may proclaim the behaviour to be more of the top-down type, although strictly speaking it depends on which window of frequencies we are looking at. In Figs. 21(a), 21(b) and 21(c), we have plotted |ρ(ω, t av , k) − ρ f (ω, k)| as a function of t av for various values of ω at |k| = 0, 1 and 4.5 respectively for α = 0.05. From these plots, it is directly visible that the patterns of thermalisation vary with ω in a complex way, and do not coincide even at late time as we would have naively expected. It is also clear from these figures that even in the high frequency region we do not get a quasinormal mode type ring-down as present in the one-point function. The reader may readily contrast these plots with Fig. 2(c) where the quasinormal mode ring-down of the one-point function has been presented. Note in the latter case, the envelope of decay in time has a linear slope in the log plot and there are also oscillations with uniform periods -both the decay slope and oscillation period are controlled by the lowest quasinormal mode of the final black brane. In case of ω = 10, 20 and 30, both the period and amplitude of the oscillations in the spectral function get modulated in complex ways even at late time t av ≈ t α av(f) as visible from the plots.
Although we find complicated patterns of thermalisation of the spectral function at various values of ω, it may still be possible that large scale features of time-dependence are controlled by the lowest quasinormal mode of the final black brane for all ω. We need more numerical accuracy, data and computational time in order to see such features, if they exist.
Interestingly, if we plot the modulus of the full subtracted retarded propagator |δG R (ω, t av , k)| on a log scale as before, where δG R = G R − G R f with G R f being the final thermal form of the retarded propagator, then even at α = 0.05 and |k| = 0 the oscillations are suppressed (see Fig. 22). Furthermore, at late time the contribution remains more or less uniform up to frequencies ω ≈ M f and then seems to decrease exponentially. Thus how close the pattern is towards a top-down description may depend on our actual observable as well as on the specific window of frequencies at which we are making our observations.
In order to understand better how different frequency bands thermalise, the Gabor transform as utilised in [19] could be a valuable tool. 17 Although the conventional Wigner transform sheds more light on the global aspects of thermalisation, and spectral weight distribution across all frequencies as a function of time, etc. and gives direct connections with experimental results (as discussed below), other representations like Gabor transform can shed light on the time-dependent dynamics of individual frequency bands. We would like to explore such questions in the future. Therefore, this should follow the quasinormal mode behaviour as clearly demonstrated in Fig. 6. However, after the Wigner transformation, the behaviour is not expected to be the same. 17 We thank the referee for pointing this out.

Conclusions
Let us sum up lessons learnt from the results presented so far and point out some interesting questions that merit further investigation. The most important lesson is that we see clearly distinct patterns for the thermalisation of the non-equilibrium spectral function as we vary parameters determining the dynamics of the non-equilibrium state even in cases where the one-point functions thermalise in a universal way via a thermal relaxation mode. Indeed, it has been observed before that even if we go beyond linear response regime in holographic theories, the non-equilibrium evolution of the one-point function of a scalar operator is very well approximated by thermal linear response to a very good degree of approximation even close to initial stage when the system is far from equilibrium [39,52]. Furthermore, in [53] it has been shown that even at early times, i.e. before the endpoint of the quench, the one-point functions and space-like separated correlations behave in certain universal ways irrespective of the quench protocol in general conformal field theories. In our example constituting of states approximated by AdS-Vaidya spacetimes, the one point function do not show visibly distinct pattern as it makes transition from the thermal response of the initial black hole at early time to that of the final black hole at late time. However, the non-equilibrium spectral function knows about what is going on in the bulk more clearly.
Depending on the compactness of the shell driving the expansion of the non-equilibrium horizon (reflecting the duration of the quench), it does show very different patterns of thermalisation.
The second lesson learnt is that in order to find different patterns of thermalisation of two-point functions, we should go beyond geodesic approximation in the bulk, and of course also vary the parameters determining the non-equilibrium dynamics in the state. Going beyond the geodesic approximation, we find kink formation and advanced timedependence in particular (as discussed before) and also finer structure in the late-time oscillations when they exist. Our method is certainly well suited to investigate beyond the simple example considered here. It can be readily applied even to situations where the holographic geometry dual to the nonequilibrium state can be found only numerically. It is worth noting that some of the most generic features of non-equilibrium dynamics including quasinormal modes in the background gravitational geometry at late times are not present in the simple examples we have investigated here.
The most experimentally relevant set-up that can be modeled via holography is related to solid state pump-probe spectroscopy (see for instance [54] for results of such experiments on cuprates). In this case, the material is driven away from equilibrium by a strong laser pulse, and a weak laser pulse is sent after an adjustable delay (which can be as short as an attosecond). The latter gets modulated by the medium revealing information about the time-dependence of the current-current retarded propagator. Holographically, the nonequilibrium state can be readily modelled by a dual geometry that is excited by a specific source at the boundary (which needs to be obtained numerically). Our prescription can then be applied to obtain the non-equilibrium retarded propagator. We would like to investigate such holographic set-ups in the future.
Recently in [55], it has been shown that even in the geodesic-like approximation the time-dependent behaviour of non-local observables like entanglement entropy depend strongly on the initial conditions given by the geometry of colliding shockwaves in holographic models of heavy-ion collisions. Such qualitative features are not prominently visible in the one-point functions, as for instance in the energy-momentum tensor. It should be interesting to understand if the time-dependence of the holographic spectral function can indeed reveal the detailed geometry of the colliding shock waves. In that case, it would also be interesting to understand how these time-dependent spectral functions can be measured in heavy-ion collision experiments.
A bigger question could be as follows: can we make a broad classification of patterns of thermalisation of the non-equilibrium Schwinger-Keldysh propagators in strongly interacting many-body systems? This is analogous to classifying hydrodynamic flows qualitatively based on the Reynolds number. The analogue of the Reynolds number in this case should be appropriate combinations of external parameters driving the nonequilibrium dynamics and intrinsic parameters (like masses and couplings of the bulk theory, or equivalently the spectrum of primary operators and their three-point functions in the dual CFT) specifying the system. Indeed holography can lead us to such a broad classification which may be compared with experimental results.
Finally, we have proposed a method only for obtaining causal holographic Schwinger-Keldysh propagators at and away from equilibrium. We need to extend our method further for obtaining the other Schwinger-Keldysh propagators (specifically the Feynman propagator). This will give us insight into generalisation of fluctuation-dissipation relations away from equilibrium, and complete information about thermalisation and decoherence. In order to achieve this, we need to go beyond purely causal response by making the boundary sources in holography dynamical in a self-consistent way as in semi-holography (see [56,57] in particular). We will like to generalise our method using the latter framework in future.