Thermalization of the spectral function in strongly coupled two dimensional conformal field theories

Using Wigner transforms of Green functions, we discuss non-equilibrium generalizations of spectral functions and occupation numbers. We develop methods for computing time-dependent spectral functions in conformal field theories holographically dual to thin-shell AdS-Vaidya spacetimes.


JHEP04(2013)069
1 Introduction The problem of thermalization of a quantum system that is prepared in a highly excited state is ubiquitous in physics. In contrast to situations where the equilibrium state of a quantum system is slightly perturbed and its return to equilibrium can be described by linear response theory, there is no general approach to the problem of thermalization starting far from equilibrium. In some cases, e. g. when the system is dilute and composed of weakly coupled quasiparticles, kinetic theory provides a general solution. However, many quantum systems of interest are either dense or strongly coupled, or both, and kinetic theory does not apply. A widely studied example is the strongly interacting matter produced in high energy nuclear collisions [1]. When thermally equilibrated, this matter is known as a quark-gluon plasma, but the collision is thought to produce an initial state composed of densely packed, highly excited, but locally coherent gauge field domains, often referred to as a "glasma" [2]. The thermalization of this initial state is an open problem of intense phenomenological interest [3].
The late time evolution of the energy-momentum tensor of a thermalizing quantum system is described by the effective low-energy theory of energy and momentum, fluid dynamics. Chesler and Yaffe [4] and Heller et al. [5,6], have shown that fluid dynamics becomes a good approximation rather quickly. But this "hydroization" of a quantum system does not imply full or even approximate thermalization of the system, because the anisotropy of the energy-momentum tensor may remain quite large, and there is no assurance that the expectation values of other observables are near their thermal values.
Time evolution from near-equilibrium initial states in a strongly coupled gauge theory with a holographic AdS/CFT dual has been extensively studied (see the review [7]). More recently, focus has moved to thermalization from far off-equilibrium initial states. One interesting quantity to consider in such studies is the holographic entanglement entropy [8,9], which was generalized to time-dependent backgrounds in [10], in particular the collapse of a homogeneous massless shell in the AdS space-time, leading to black hole formation in the bulk, which corresponds to thermalization in the dual field theory [11]. Thermalization was studied in more detail in [12][13][14][15][16][17] following the time evolution of entanglement entropy, Wilson loops and equal-time correlation functions. The most nonlocal quantity, the entanglement entropy, reached thermal equilibrium last, but still very fast, on a time scale saturating the causality bound [18,19]. A complementary approach to thermalization was investigated in [20][21][22], who considered thermalization after an anisotropic deformation of the boundary spacetime [23], in particular how fluctuations created near the horizon and dissipation come to a balance to satisfy a generalized fluctuation-dissipation theorem. Isotropization and thermalization were also investigated in [24], by considering a large set JHEP04(2013)069 of bulk initial conditions. Rotationally symmetric isotropization including radial flow was recently investigated in [25].
Here we extend this previous work by considering time-like correlation functions. This allows us to generalize concepts usually studied in equilibrium to non-equilibrium states, and to compute them holographically. Two important quantities of this kind are the spectral density and the occupation number distribution of states. Time-dependent versions of these quantities are discussed and illustrated in a simple example in section 2.
(1.1) which models time evolution following the sudden injection of energy in the dual field theory. We will develop two approaches for computing two-point functions with timelike separations. The first approach, described in section 3, uses the geodesic approximation, which has been employed for computing two-point functions of spacelike separated heavy operators in [12,14,15,26]. There are significant new challenges in the present casereal timelike geodesics do not reach the boundary of AdS space (and therefore cannot be used to connect timelike separated boundary points). An alternative approach could have been to use analytic continuation from Euclidean signature, but AdS-Vaidya does not allow a standard continuation of this kind. Indeed, the absence of a conventional Euclidean continuation with which to define Green functions is a general challenge in the study of non-equilibrium systems.
Here we propose novel approaches to overcoming these obstacles within the AdS/CFT correspondence: (i) by defining a new, non-standard Euclidean continuation, (ii) by using complex geodesics, (iii) by developing a prescription for splicing propagators across the infalling shell. In a saddle point limit appropriate to heavy operators, we obtain an analytic expression for equal-space correlators that agrees precisely between all three approaches. Method (iii) does not require the probe to be heavy, and hence in section 4 we use it to obtain numerical results for retarded two-point functions in momentum space.
A complementary approach to the holographic non-equilibrium spectral function and the fluctuation-dissipation relation has been developed in [27,28]. They consider perturbative expansions (derivative and amplitude) away from thermal equilibrium, working first at leading order [27] (where they also study non-equilibrium shift effects) and then extending to higher orders [28]. In particular [28] also derives a generalized (holographic) fluctuationdissipation relation for Wigner transformed spectral function and statistical function (the anticommutator Green function) which holds for non-equilibrium states which are perturbatively connected to thermal equilibrium. It would be interesting and important to compare the approach of [27,28] to the methods and context of this paper.
2 Time-dependent spectral functions and occupation numbers

Spectral function and occupation numbers in equilibrium
Consider a system prepared in a state described by a density matrixρ. In the interaction picture or in the Heisenberg picture, the density matrix is time-independent while the JHEP04(2013)069 operators evolve in time. In order to study how thermalization proceeds, we will look at the evolution of correlation functions. The relevant correlators are the time-ordered (Feynman) correlation function and the retarded correlation function If the system is translation invariant in space and time, we can consider the Fourier transforms of the correlation functions G F (k, ω) and G R (k, ω). In Fourier space, the retarded correlation function is complex valued, and its imaginary part defines the spectral function 1 The imaginary part of the time-ordered correlator contains additional information, which depends on the mean occupation number density n(k, ω) of the Fourier modes. When the system is at thermal equilibrium, the correlation functions are related by the Fluctuation-Dissipation Theorem (FDT) with the occupation number equal to the Bose (or Fermi) distribution n(k, ω) = (e ω/T ∓ 1) −1 . The authors of [20] have emphasized the usefulness of (2.4) as a measure of the degree of thermalization in holographic theories, in the context of a model where fluctuations are created near the stretched horizon and transported to the boundary. In these works the fluctuation-dissipation theorem has also been expressed as a relation between the symmetrized and antisymmetrized correlation functions

Time-dependent generalizations
In order to determine the degree of thermalization we would like to define time-dependent generalizations of the occupation number density n(k, ω) and the spectral function ρ(k, ω). By definition, a time-dependent system in not time translation invariant, so the Feynman and retarded correlators explicitly depend on two instants of time: Space is still homogenous, so the correlators depend only on x 1 − x 2 which can be Fourier transformed to momentum k. But a simple transform from time to frequency space ω is no longer meaningful as we seek a time-dependent notion of occupation number and spectral density. We therefore resort to a well-known alternative -the Wigner transform.

JHEP04(2013)069
First introduce an average time T and a relative time t by and define T -dependent Green functions by keeping the average time T fixed and Fourier transforming in the relative time t: where we have suppressed the spatial coordinates. The T -dependent spectral function is then defined as In this paper, our goal is to study how the spectral function evolves as a function of the average time T . In order to focus on information near T , we may consider small intervals of relative time t around t = 0. To this end, we add a Gaussian window function into the integrand: 2 where σ controls the width of the time window, and G is a correlator. In the presence of a finite time window the T -dependent spectral function becomes We can now define a generalized, time-dependent notion of occupation number density by the relation In an interacting field theory we expect the occupation number density of the nonequilibrium system to evolve towards thermality (e βω −1) −1 (for bosons), so that the system eventually satisfies the fluctuation-dissipation theorem and becomes fully thermalized. We note that the notion of occupation number density is only meaningful for momenta and frequencies where the spectral function has support: if a mode does not exist 2 This is similar in spirit to the Husimi distribution, which is a Gaussian smearing of the Wigner distribution on phase space [31,32]. It is also related to a well known method in signal analysis. For a signal s(t) with time dependent frequency, it is convenient to add a Gaussian weight factor into the integrand with a small window about a mean time instant t , Then by studying sσ(t , ω) as a function of t one can track changes in frequency. This procedure is known as the Gabor transform, and it is a special case of a more general class of wavelet transformations (with more customized window functions called wavelets). In our case, we are studying functions of two time variables, and we are primarily interested in tracking changes as the untransformed variable T changes, hence we are setting t = 0. It would be an additional possibility to perform a full Gabor transform, in which case varying t would allow us to study information at different time separations, but that is beyond the scope of this paper.

JHEP04(2013)069
in the spectrum, it is meaningless to ask whether it is occupied. In physical quantities, the occupation number density will always be accompanied by the spectral function, so that no contributions originate from regions in momentum space where the spectral function vanishes. This is reflected in (2.4). The issue does not arise in the conformal field theories considered later in the paper but the simple quantum mechanical example in the following section has a discrete spectrum and in that case expressions like (2.13) need to be interpreted with care. As explained in appendix A, we can now follow the analysis in [42] to define a timedependent entropy density where s(k, ω, T ) = (n(k, T, ω) + 1) ln (n(k, T, ω) + 1) − (n(k, T, ω)) ln (n(k, T, ω)) .

(2.15)
This definition is a natural extension of the conventional microcanonical definition of von Neumann entropy of a system, and reduces to the latter in equilibrium. It should be contrasted with more specialized notions of entropy such as entanglement and Renyi entropies, which have been the subject of many studies of holographic thermalization.

A simple example: the quenched harmonic oscillator
To illustrate the out-of-equilibrium notion of a time-dependent spectral function in (2.9), and in (2.12) when a finite time window is introduced, we first consider the harmonic oscillator with characteristic frequency ω 0 . As we will eventually be interested in applying our methods to quantum field theories, it is useful to think about the harmonic oscillator as a 0 + 1 dimensional quantum field theory. This quantum field theory has a single one particle state a † |0 with energy ω 0 . Multiparticle states then appear as the excited states (a † ) n |0 . We will study the single particle spectral function obtained from the two point function of x(t) and the occupation number of the single energy eigenstate. It should be also noted that all of our results in this section generalize trivially to free quantum field theories with the substitution This is because free field theories are just a set of decoupled harmonic oscillators labeled by a spatial momentum k, and because the observables we consider are related to the two point correlation functions, which are diagonal in k-space. The retarded Green function for the oscillator in its ground state is Writing (2.19) the spectral density is given by Following the general discussion of section 2.2, we now introduce a time window into the Green function, Using (2.17), (2.18) and the known Fourier transform of a Gaussian, we find (2.23) The spectral function is then and using (2.20), we find The Gaussian time window simply results in a Gaussian smearing of the delta function peaks in the spectral density. Next we study a quenched harmonic oscillator, with frequency ω i for t < t 0 and frequency ω f for t > t 0 . The time of the quench can be taken to be t 0 = 0 without loss of generality. The retarded Green function G R (t 2 , t 1 ) now depends on both t 1 and t 2 , not only on their difference. It is given by (2.26) The matching conditions at t 0 = 0 are such that for t 1 < 0 both G R and ∂G R /∂t 2 are continuous across t 2 = 0, and for t 2 > 0 both G R and ∂G R /∂t 1 are continuous across t 1 = 0.

JHEP04(2013)069
In terms of the average time and relative time coordinates, t 1 ≡ T − t/2, t 2 ≡ T + t/2, the retarded Green function reads We can then study its Wigner transform (2.8) and obtain the time-dependent spectral function (2.9). The main features of ρ(T, ω) can be uncovered by inspecting our expression (2.27) for G(T, t). Consider, for instance, the very early time regime, T → −∞. In this limit, G(T, t) is given by the first line in (2.27), and ρ(T, ω) reduces to the spectral function for a time-independent harmonic oscillator with frequency ω i . Similarly, in the late time limit T → ∞, we recover the spectral function of a harmonic oscillator with frequency ω f . The time-dependent spectral function for a quenched harmonic oscillator thus interpolates between those of harmonic oscillators with the initial and final frequencies.
To get a first idea of the intermediate-time regime, we can take the average time to be the transition time, T = t 0 = 0. In this case (2.27) reduces to (2.28) and we see that at T = 0 the retarded Green function executes beats between the initial and final frequencies in relative time.
The Fourier transform of the retarded Green function with respect to relative time can be carried out explicitly for general values of the average time. Taking the imaginary part then gives the following time-dependent spectral function, We observe that the time-dependent definition (2.9) allows for negative values of the spectral function for ω > 0. This is to be compared with the equilibrium spectral function that is positive definite for positive ω.  Although it is not necessary, it may be useful to introduce time windows as in (2.11) in order to define generalized time-dependent spectral functions that only depend on the behavior near T . This will have a smearing effect, for instance on the various delta functions appearing in our quenched harmonic oscillator example as can be seen in figure 1. Observe in particular that the spectral function can be made positive definite for ω ≥ 0 for an appropriate value of the width σ (figure 1 right). This is similar to the situation with conventional phase space Wigner density which attempts to define a notion of particle density in the phase space of a quantum theory. The Wigner density is sometimes negative, but the associate Husimi density, obtained by smearing with an appropriate Gaussian is everywhere positive [31,32].
Examining the smeared spectral function for different times (figure 2) shows that ρ σ evolves between the spectral functions of harmonic oscillators with the correct initial and final frequencies. quickly evolves to its asymptotic value, which is given by the expectation value of the harmonic oscillator number operator (for oscillator frequency ω f ) evaluated in the ground state of an oscillator of frequency ω i , which is the state of the system at the time of the quench. The oscillations around the final value arise from the smearing due to the finite time window -their size depending on the choice of the parameter σ.

The geodesic approximation in the AdS/CFT correspondence
To use the methods developed above to discuss time dependent spectral functions in the AdS/CFT correspondence we need to first compute two point functions in the latter setting. In this section, we use a geodesic approximation to compute two-point functions with timelike separation in a two dimensional CFT dual to the AdS 3 -Vaidya spacetime (1.1). This approximation is accurate for operators of large conformal dimension in the CFT. We will find simple analytic results for the special case of equal-space two-point functions, but our methods generalize to generic separations. These results allow us to compute an integral over spatial momenta of the time-dependent spectral function of high-dimension operators, as well as a weighted integral of time-dependent occupation numbers.
Consider a scalar operator O(x, t) with conformal dimension ∆ in the dual CFT. The time-ordered two-point function O(x 1 , t 1 )O(x 2 , t 2 ) is given by a path integral over all paths P that connect the two insertion points (x 1 , t 1 ) and (x 2 , t 2 ) on the boundary: For large ∆, we can use a saddle point approximation, in which the sum over all paths can be approximated by a sum over all geodesics connecting the boundary endpoints [33] T

JHEP04(2013)069
where L denotes the geodesic length. However, due to contributions near the asymptotically AdS boundary, the geodesic length between two boundary points contains a universal divergence and needs to be renormalized. Throughout this section, we will define a renormalized length δL ≡ L − 2 ln(r 0 ), in terms of the bulk cut-off r 0 , by removing the divergent part of the geodesic length in pure AdS (see eq. (B.16)-(B.17) in appendix B). The renormalized two point function can then be approximated by where δL is the renormalised length of the geodesic that connects the points (x 1 , t 1 ) and (x 2 , t 2 ) on the boundary. The geodesic approximation works most straightforwardly in Euclidean spacetimes, or in static spacetimes that can be Wick rotated to Euclidean signature. In [14,15], this approach was used in a time-dependent setting, to determine equal-time two point functions in the field theory dual of the Vaidya spacetime. This work was further extended in [26] to non-equal-time two-point functions with spacelike separations. In all these cases, the relevant bulk geodesics P were spacelike, and the geodesic length L(P) was real.
New complications occur when the points on the boundary are timelike separated, in which case we might expect to connect them using timelike geodesics. The main problem is that (real) timelike geodesics in an asymptotically AdS spacetime never extend to the boundary. For example, in the case of empty AdS 3 with metric the timelike geodesics (e 2 > j 2 ): satisfy r(λ) ≤ e 2 − j 2 , ∀λ. Thus they do not reach the boundary, where r → ∞.
In fact, this is a specific realization of the general difficulty of computing correlation functions in time-dependent spacetimes in a semiclassical approximation. Given that a general prescription is not available, we will here propose and investigate three potential approaches: 1. For time-independent spacetimes the standard prescription for vacuum correlators is to calculate in Euclidean signature and then analytically continue to Lorentzian signature. Since there is no obvious Wick rotation of the AdS-Vaidya spacetime, we will propose and use a non-standard Euclidean continuation in section 3.1.
2. We will consider complexified geodesics in section 3.2. In this approach, geodesics between timelike separated boundary points make excursions in a complexified AdS-Vaidya geometry before returning to real points on the boundary.

JHEP04(2013)069
3. In section 3.3 we will directly connect the AdS and black hole retarded bulk-toboundary propagators across the infalling shell, and will apply a saddle point approximation to the resulting integral expression. This approximation does not directly rely on the geodesic approximation but shares the same regime of validity of large scalar masses.
For equal-space correlators, we will find the same simple analytic expression in all three approaches.

Continuing to Euclidean signature
In appendix B, we review how AdS and BTZ can be Wick rotated to Euclidean signature, and how various two-point correlation functions of high-dimension operators in the dual field theory can be computed using bulk geodesics. Here we develop a similar approach for AdS-Vaidya.

Wick rotation of AdS-Vaidya
Recall that the thin-shell AdS-Vaidya metric represents an infalling shock wave of 'null' dust in AdS 3 . Outside the shock wave (v > 0), the metric is given by the BTZ solution while inside (v < 0), we recover the pure AdS 3 metric While for static spacetimes, such as AdS 3 or BTZ, a Wick rotation of the metric is straightforward, this is no longer the case for an explicitly time dependent geometry. Thus we propose the following procedure: • We formally take the Vaidya metric to be a limit of a family of "spacelike Vaidya" metrics. The latter need not be solutions to the equations of motion with conventional matter -they are auxiliary objects used to regularize the computation.
• This spacelike Vaidya metric can be Wick rotated to Euclidean signature by a simultaneous rotation of the time coordinate, as well as of an additional parameter.
• In this auxiliary Euclidean spacetime, correlation functions can be computed in a geodesic approximation. The result can then be Wick rotated back. Taking the limit we find our proposed result for the conventional AdS-Vaidya space-time.
While the results we find will look reasonable, it is not a priori obvious that this procedure should be valid. This is why, after discussing the present method in detail, we will also develop two alternative prescriptions. All three prescriptions match, suggesting that the method of considering a time-dependent space-time as a limit of others with standard Euclidean continuations may be more generally useful.
Euclidean Vaidya. On the Lorentzian 'spacelike' Vaidya metric (3.9), we can now perform an analytic continuation, on the time coordinate z = iv, as well as on the parameter Z = iE. Without loss of generality, we can take Z > 0 and find the metric where as before θ(z/Z) = θ(z). Note that the radial coordinate r now runs from the metric becomes which is manifestly positive.

Geodesic length in Euclidean Vaidya
We now compute the length of geodesics that start at a point (r 0 , x 1 , y 1 ) and end at (r 0 , x 2 , y 2 ), where r 0 denotes the location of the regularized AdS boundary. We assume ∆y ≡ y 2 − y 1 > 0 and, for computational simplicity, mostly focus on geodesics with ∆x ≡ x 2 − x 1 = 0. Geodesics in the Euclidean Vaidya geometry can be constructed by gluing together at the shell, with the appropriate refraction condition, geodesics in Euclidean AdS

JHEP04(2013)069
and Euclidean BTZ. The latter have been worked out explicitly in appendix B and read for z > 0, andx for z < 0. Here ± denote two separate branches (which appear above and below the turning points of AdS or BTZ geodesics) and x 0 , Γ ± , y 0 , λ 0 ,x 0 , Σ, σ,ȳ 0 ,λ 0 are constants. We can distinguish three cases.
• If both endpoints occur before shell injection (y 1 < y 2 < 0), there is a geodesic connecting them which is entirely within Euclidean AdS. Any additional geodesics would require the existence of geodesics in BTZ that connect two points on the shell, and these do not exist.
• If both endpoints occur after shell injection (0 < y 1 < y 2 )), there is a geodesic connecting them which is entirely within Euclidean BTZ. 3 3 We can check that no new geodesics are admitted in the presence of the shell as follows. Consider a geodesic with both endpoints in Euclidean BTZ, i.e. after shell injection (0 < y1 < y2). If the geodesic has to extend into the z < 0 region, it will have to cross the shell twice. It will therefore consist of two z+(r) branches, which are connected by a geodesic in AdS3.z−(r) thus needs to have a minimum, located at For the special case of geodesics that connect points on the boundary with ∆x = 0, one has cos(σ) = 0, such that r = 0. However, since r > √ Z 2 + R 2 ,z−(r) reaches no minimum and the geodesic is thus completely in Euclidean BTZ. The renormalized geodesic length is worked out in eq. (B.18) and reads

JHEP04(2013)069
• Finally consider geodesics with y 1 < 0 and y 2 > 0, and such that ∆x = 0. They cross the shock wave once and hence we need to determine how they refract at the shell location in order to compute their length.
Boundary conditions at the shell. Consider a point P in (P out ) just inside (outside) the shell, and let the coordinate difference between P in and P out be ∆X = (∆z, ∆r, ∆t). Take another point M on the shell z = 0, and let the coordinate difference between P in (P out ) and M be ∆X in = (∆z in , ∆r in , ∆t in ) (∆X out = (∆z out , ∆r out , ∆t out )), so that ∆X in + ∆X out = ∆X. Then the distance from P in to P out via M is where f in and f out where defined in (3.12). We want to find the point M that minimizes ∆s. Extremizing this with respect to ∆r in and ∆x in we find and therefore In the ∆X → 0 limit, we obtain the refraction law that is supplemented by the continuity conditions x in (r * ) = x out (r * ) and z in (r * ) = z out (r * ) = 0, where r * is the value of the radial coordinate at which the geodesic reaches the shell.

JHEP04(2013)069
Geodesics crossing the shell. If y 1 < 0 < y 2 , the geodesic starts in Euclidean AdS, crosses the shell once and ends in Euclidean BTZ. The geodesics are given respectively by the expressions (3.15) and (3.14). Imposing the boundary conditions results in the following relations Together with the refraction condition (3.21) that leads to and the continuity conditions at the shell these fix all unknowns (x 0 , y 0 ,x 0 ,ȳ 0 , Σ, σ, Γ + , Γ − and r * ) and determine the renormalized geodesic length where in the last line we have used the continuity condition λ + (r * ) =λ + (r * ).

Two-point functions of thermalising CFT
As reviewed in appendix B for the cases of AdS and BTZ, by Wick rotating the Euclidean two point function G E (x 2 , y 2 ; x 1 , y 1 ), one finds the time-ordered (Feynman) twopoint function 4 The retarded two-point function can be obtained from the time-ordered one by

JHEP04(2013)069
From (3.33) we then find (the geodesic approximation of) the Feynman and retarded two-point functions in the thermalizing CFT dual to 3-dimensional AdS-Vaidya: . (3.37) In the limit R → 0 or t 2 → 0 + , G R (t 2 , x; t 1 , x) reduces to the retarded vacuum two-point function (B.27) with ∆t = t 2 −t 1 or ∆t = −t 1 respectively, while for t 1 → 0 − to the thermal correlator (B.28) with ∆t = t 2 . The above formulas are valid for non-integer values of ∆. For a discussion of integer values of ∆, see appendix C.
To compute time-dependent spectral functions and occupation numbers, we would need to Fourier transform the relative times and positions in the retarded and time-ordered two-point functions. This requires knowledge of the position space two-point functions for arbitrary separations, which are more cumbersome to obtain than the equal-space two-point functions we have explicitly computed so far. We will therefore focus on the information we can extract from the equal-space two-point functions, i.e. the time-dependent spectral function integrated over all spatial momenta.
Following the discussion in section 2, we introduce relative and average time coordinates t ≡ t 2 − t 1 and T ≡ (t 1 + t 2 )/2 and Fourier transform G R (t 2 , x; t 1 , x) with respect to t. The imaginary part of the result is where the additional factor of 2π with respect to the definition (2.9) appears because the retarded two point function is here given in position space, rather than momentum space. Numerical results for the spectral function integrated over momenta are plotted in figure 4 for different times. As time evolves, the curves interpolate between the vacuum result (in dotted blue in panel (A)) and the thermal one (in continuous red). As the spectral function is odd in ω, we only plot the range ω ≥ 0. Panel (B) plots, as a function of time and for different values of scaling dimension ∆, the slope β of the spectral function in the range of linear growth close to ω = 0. For a fixed frequency ω, ρ(ω, T ) does not increase monotonically from the vacuum to the thermal value, but it oscillates.
As in the example of the quenched harmonic oscillator we discussed in section 2.   ducing an appropriate time window as in eq. (2.11)-(2.12), we can smear the curves and make the spectral function positive definite for ω > 0 (see figure 5). From the time-ordered two-point function, using definition (2.13), we can also obtain time-dependent occupation numbers weighted by the spectral function and integrated over all spatial momentañ(T, ω). These are plotted in figure 6 as a function of time, without (in (A)) and with (in (B)) the Gaussian smearing.

Complexified geodesics
To investigate the validity of the non standard Euclidean continuation used in section 3.1, we here repeat the computation of the geodesic length in the 3-dimensional AdS-Vaidya spacetime using a different approach. While real timelike geodesics can not connect timelike separated points on the boundary, complex geodesics can. The latter have been studied in [34][35][36] and we review their construction in both AdS 3 and BTZ in appendix D. As we show below, the geodesic length obtained with this second method coincides with (3.33).
As before, we can construct Vaidya geodesics that cross the shockwave by joining together AdS and BTZ geodesics, with an appropriate refraction condition at the shell.
while for v < 0, v is defined as v = t − 1 r and the geodesics are given by The refraction condition at the shell has been computed in [14,15] and reads where the second condition is trivially satisfied for equal-space geodesics. Continuity at the shell fixes where we have defined Λ 1 ≡ λ * − λ 1 + iβ 1 and Λ 2 ≡ λ * − λ 2 + iβ 2 . The boundary conditions are and using the continuity condition (3.47), the refraction condition can be written as which corresponds exactly to the result (3.33) that we found using the analytic continuation to Euclidean signature.

A saddle point approximation for joining propagators
In this section we consider a third way of calculating the equal space retarded propagator in the AdS-Vaidya background. Specifically we will consider the saddle point approximation, which is valid in the limit of large scaling dimension ∆, of the following integral formula for the CFT correlation function 56) where t 1 < 0 < t 2 , and G B,R andG B,R are propagators that we will describe in more detail below. First, we have switched to coordinates (x, v, z) in which the AdS-Vaidya metric reads The time coordinate t appearing in (3.56) is the boundary value of the Eddington-Finkelstein time v.

JHEP04(2013)069
In Lorentzian signature one has to distinguish between bulk-to-boundary propagators and boundary-to-bulk propagators. The retarded boundary-to-bulk propagator G B,R having its initial point at the boundary and the final point in the bulk, evolves initial data from the boundary into the bulk. On the other hand, the retarded bulk-to-boundary propagator G B,R takes bulk data and evolves it towards the spacetime boundary. Below we give a precise definition of these two-point functions in terms of the bulk-to-bulk propagator.
In (3.56), the two sided derivative is understood to act only inside the square brackets, G AdS B,R is the retarded boundary-to-bulk propagator in AdS 3 andG BTZ B,R is the retarded bulkto-boundary propagator in BTZ. The integral in (3.56) is performed at the position of the shell v = 0 and g is either the BTZ or the AdS 3 metric. 5 The basic logic is as follows. A delta-function source inserted on the boundary before the shell injection produces a disturbance that propagates to the shell according to G AdS B,R . After passing through the shell, the disturbance propagates back to the boundary according toG BTZ B,R . In the following we provide a derivation of (3.56).
Derivation of eq. (3.56): first consider the retarded bulk-to-bulk propagator G BB,R associated to a scalar field of mass m in an asymptotically AdS spacetime. It solves the equations vanishes when v < v and satisfies the standard normalizable boundary condition at the AdS boundary. 6 Given the initial data Φ(x , v = 0, z ) = Φ 0 (x , z ), for v > 0 we will show that the unique solution to the massive wave equation satisfying the boundary condition of vanishing non-normalizable mode is given by For v > 0, it is clear from the definition of the bulk-to-bulk propagator (3.58) that (3.59) solves the massive wave equation. By integrating (3.58) with respect to v from (3.60) 5 For the above formula to hold it is necessary for √ −gg vµ to be continuous at the junction v = 0. As this is the case in the Vaidya spacetime, using either the BTZ or the AdS3 metric gives the same result. 6 Here we focus on scalar fields satisfying the standard normalizable boundary conditions Φ ∝ z ∆ as z → 0, but generalizations to more general boundary conditions are in principle straightforward.

JHEP04(2013)069
The equations of motion (3.58) imply that G BB,R is a regular function of v near v = v . Thus, the mass term as well as the ∂ x ( √ −gg xx ∂ x G BB,R ) term can be assumed to be regular near v = v leading to order contributions to the integral in (3.60). The terms singular near v = v in the integrand come from the derivative terms which can be written as where we have dropped the explicit coordinate dependence. Using the fact that √ −gg vz is a continuous function of v, and assuming that ∂ z G R is a regular function of v, 7 the last term above gives a contribution of order to the integral in (3.60). Finally the integral (3.60) can be performed to give where we used the fact that G R vanishes when v > v . Thus, we see that the bulk-to-bulk propagator satisfies as v → 0; so that indeed (3.59) satisfies the correct initial condition Φ(x, v = 0, z) = Φ 0 (x, z).
In particular we can use (3.59) to join two retarded bulk-to-bulk propagators

(3.65)
To obtain the boundary theory retarded correlation function G R from a bulk-to-bulk propagator we can use the well known relation (see e.g. [11,36]) where t is the boundary value of v. From (3.65) we then obtain 67) where we have identified the retarded boundary-to-bulk propagator as [11] G B,R (x, v = 0, z; x 1 , t 1 ) = 2ν lim and the retarded bulk-to-boundary propagator as

JHEP04(2013)069
Application to AdS-Vaidya: in the specific case of AdS 3 Vaidya, (3.67) reduces to (3.56). As reviewed in appendix E, the retarded boundary-to-bulk propagator in AdS 3 spacetime is given by and the retarded bulk-to-boundary propagator in BTZ is given bỹ in terms of the coordinates (3.72) The above coordinate transformation maps BTZ locally into AdS 3 , as reviewed in appendix E. In the new coordinates, the shell is located on a surfacev = 0 where the coordinate relations reduce to x = 1 R log Rx and z =z/(Rx). Substituting (3.70) and (3.71) and performing the change of coordinates (3.72), the integral (3.56) becomes to leading order in ∆. For simplicity here we consider non-integer ∆ -the extension to integer ∆ is described in appendix C. The integration contour is determined by the i factors in (3.70) and (3.71) and in general it can be rotated to pass through the saddle point For the equal space correlation function (x 1 = x 2 = 0) the saddle point is at .

JHEP04(2013)069
which agrees (up to an undetermined normalization) with the geodesic approximation result (3.37) for t 1 < 0 < t 2 . To make this argument precise one would need to check that the saddle point is on a steepest descent contour, which is nontrivial in the presence of branch cuts.

Beyond the geodesic approximation
In this section we will develop an alternative method for computing propagators in AdS-Vaidya that goes beyond the geodesic approximation and enables us to compute the result for all spatial momenta k and for any conformal dimension ∆. This method allows us to compute the full spectral function but only gives numerical results.

Holographic dictionary without time translation invariance
Usually the holographic dictionary is discussed in translation invariant backgrounds. The Green functions then depend only on the difference x − x of two points. When translation invariance is broken in at least one direction, the Green functions depend explicitly on two points. A standard holographic method to compute the retarded Green function in the field theory at the boundary is to study the dual bulk field: choose an infalling boundary condition at a (coordinate) horizon, solve the field profile from the equation of motion, and then extract the Green function from the asymptotic behavior. So we will begin with a brief discussion of the holographic dictionary for one-point and two-point functions in the case of broken translation invariance. Consider a free boson in 3d asymptotically AdS (AAdS) space of action where the counterterm action is given by and γ is the induced metric on the z = surface. 8 In the field theory on the boundary, using standard results of linear response theory applied to AdS/CFT [38][39][40] and large-N factorization, the response of an operator at (x 2 , t 2 ) in the presence of a delta function source at (x 1 , t 1 ) is given by the retarded correlation function 9 The details of the derivation of (4.3) are included in appendix F. On the other hand, according to the standard holographic dictionary, the one-point function in the boundary theory is given by The counterterm (4.2) is applicable for the range of scaling dimensions ∆ < 2, to which we specialize in the following numerical calculations. For larger ∆, and for more complicated bulk field content, one needs to add more counterterms (see e.g. [37]). 9 Here we are assuming that O does not have a vacuum expectation value.

JHEP04(2013)069
where J and φ + denote respectively the leading and subleading asymptotic behavior of the bulk scalar field φ as z → 0 where ∆ = 1 + ν, ∆ − = ∆ − 2ν and ν = √ 1 + m 2 . Comparing (4.3) and (4.4), we obtain the two point function where φ +,δ is the subleading mode of a bulk solution that satisfies the boundary condition J(x 2 , t 2 ) = δ(x 2 − x 1 )δ(t 2 − t 1 ) for the non-normalizable mode. A similar formula for the retarded correlator in momentum space in a time translationally invariant state was obtained in [40]. For more discussion on the relationship between the one point function procedure and the more conventional prescription for obtaining retarded correlators [39], see [40]. To construct such a solution, it is more convenient to specify the boundary data and reconstruct the full bulk field profile using the retarded boundary-to-bulk propagator rather than specifying the boundary condition at the Poincaré or event horizon and integrating towards the boundary.
In the case of AdS 3 spacetime, the boundary-to-bulk propagator is given by (4.8) where C ∆ = Γ(∆)/(πΓ(∆ − 1). Using (4.7) we obtain the solution φ, whose normalizable mode can be read off, which after using (4.6), leads to the CFT retarded two point function This also gives us the two point function in the CFT state dual to the Vaidya spacetime before the shell for t 1 < 0 and t 2 < 0.
In the next section we apply this method to the Vaidya background, the only difference with respect to the pure AdS setup being that to construct the bulk field profile one needs to continue the solution across the null shell.

Numerical construction of the CFT retarded two-point function
As discussed above, to obtain the retarded boundary correlator, we need to solve the scalar equation of motion (2 − m 2 )φ = 0 in the Vaidya metric background. We parametrize the geometry using ingoing Eddington-Finkelstein (EF) coordinates and the metric (3.57) and define h(z, v) = 1 − θ(v)R 2 z 2 . (4.10)

JHEP04(2013)069
Since the metric is translation invariant in the spatial coordinate direction x, it is convenient to perform a Fourier transform. The equation of motion for the scalar field then becomes This equation is first order in the EF time v, and the horizon is a characteristic surface. This means that a solution to the equations of motion is uniquely specified by a single initial condition at a constant v slice, together with a boundary condition at the AdS boundary. Also, this means that the only condition one has to impose at the position of the infalling shell v = 0 is the continuity of the field φ. This makes the collapsing null shell simpler to analyze than the case of timelike collapse, where one needs to specify two independent junction conditions at the (time-dependent) position of the shell. Due to its causal nature, for t 1 > 0 the retarded Green function is automatically thermal. Here we focus on correlation functions in the boundary theory with the initial point (x 1 , t 1 ) in the zero temperature (AdS vacuum) region (t 1 < 0). Note that it is simple to start with the initial condition in the vacuum region: the bulk scalar field solution before the collapse is given by the standard retarded AdS vacuum boundary-to-bulk propagator. It solves the equation of motion (4.11), it approaches a delta function at the boundary and vanishes outside the future lightcone of (x 1 , t 1 ).
The Fourier transform to mixed (k, v, z) space of the retarded AdS 3 boundary-to-bulk propagator is (see appendix G) . (4.13) In particular, at the shell v = 0 the field profile is given by and (4.12) indeed satisfies the asymptotic behavior (4.5) with J(t) = δ(t − t 1 ): where we recall v reduces to t on the AdS boundary. In order to construct the full field profile everywhere, we must extend it into the region after the null shell collapse (v > 0). As we know the behavior of the field at v ≤ 0

JHEP04(2013)069
analytically, the calculation reduces to solving the first order equation of motion (4.11) in the BTZ background v > 0, with initial data specified by (4.14). In the following we solve (4.11) numerically. For convenience we choose the radius of the event horizon to be R = 1 by a conformal rescaling.
For v > 0, the scalar field equation (4.11) can be solved noting that it can be written as which can be integrated with respect to z to give For computational convenience we define a rescaled field P (k, v, z) as in terms of which (4.17) becomes To solve this equation numerically, we discretize both space and time where (h v , h z ) = (T /N v , 1/N z ), N v , N z are the numbers of lattice sites and T is the range of the time coordinate. The field P (v, z) is replaced by discrete variables derivatives are replaced by discrete derivatives 22) and the integral by a sum By knowing the value of the field P at some initial constant v slice, we can use the discretized version of (4.19) to obtain the value of the field at the next time step. Iterating this procedure we obtain the time evolution of the scalar field. Alternatively, a similar discretization and integration procedure is efficiently implemented by the Mathematica command NDSolve. The results of the two computations agree and the plots we present in the next section have been produced using Mathematica.

Retarded two-point function at fixed momentum
The left panel of figure 7 shows the time evolution of the retarded two-point function for fixed momenta (in blue continuous), together with the time evolution of the vacuum (dashed red) and thermal (dot-dashed black) correlators. The two-point function in the thermalizing CFT decays faster in time than the vacuum correlator.
A logarithmic plot reveals that the absolute value of the correlator oscillates, with the enveloping curve decaying exponentially in time (right panel of figure 7). Ingoing field fluctuations in a black hole background are well known to have a quasinormal behavior with a set of discrete oscillation frequencies with negative imaginary parts. These modes have been identified in the AdS/CFT correspondence as determining the approach to thermal equilibrium in the level of linear response [41]. Even though in our case we are considering dynamics far from equilibrium, it seems reasonable that at least at large times the dynamics of the field is dictated by the lowest quasinormal modes. The right hand side of figure 7 compares the exponential decay of the field modes with those of the lowest quasinormal mode, which in the black brane background is known to be given by [39] ω quasi = k − i∆. (4.24) Also the oscillation frequency of the correlator is given by the real part of the lowest quasinormal mode, here k.

The momentum dependence of the retarded autocorrelator
Numerical results for the thermalizing retarded correlation function as a function of momentum, at fixed time, are shown in figure 8. The plot reveals a modulation by the power law k ν at large momentum which is also present in the conformal field theory vacuum For small momenta the thermalizing retarded correlator in figure 8 differs significantly from the thermal correlator. In particular, the thermal correlator is a non-vanishing function of k at small k, while the thermalizing correlator vanishes at small k. This non thermal form at small momenta is in qualitative agreement with the "top-down" thermalization pattern obtained from many of the holographic thermalization studies.

Spectral function
As in sections 2.2 and 3.1.3, we introduce a mean time T = t 1 + t 2 and a time separation t = t 2 − t 1 . Fourier transforming G R (k, T, ω) with respect to t we obtain the time-dependent spectral function ρ(k, T, ω) = −2 ImG R (k, T, ω) . (4.27) In figure 9, we display the frequency dependence of the spectral function for fixed momentum and various average times T , obtained by Wigner transforming our numerical retarded correlator. Unlike in the harmonic oscillator example, the spectral function here is a smooth function with no quasiparticle peaks. In the vacuum the spectral function has the simple form (B.37): ρ ∝ θ(ω − k)(ω 2 − k 2 ) ν . The momentum k determines the frequency when the spectral function obtains support. In the thermal state, the spectral function has support all the way down to ω = 0 and at large ω approaches the vacuum form. The time dependent spectral function in figure 9 smoothly interpolates between the vacuum and thermal forms.

Conclusions and outlook
We have studied some aspects of thermalization following a holographic quench in the strongly coupled quantum field theory dual to the AdS 3 geometry. In this scenario, energy is suddenly injected in the ultraviolet in the form of a thin shell of null energy that falls in the bulk AdS 3 space eventually forming a BTZ black hole. In particular, we have investigated unequal-time two-point functions of the boundary field theory via their holographic duals, either in the semiclassical approximation by calculating the length of the bulk geodesic connecting the endpoints of the boundary correlator, or by solving the appropriate wave equation in the time-dependent bulk geometry. One quantity of particular interest is the spectral function of the correlator, which changes from the vacuum spectral function to the thermal spectral function during thermalization. While causality dictates that the vacuum spectral function ρ(ω, k) exhibits a gap ω min = k caused by the absence of space-like excitations, the thermal spectral function exhibits no gap due to the presence of real excitations in the thermal state. In the case of the quench followed by thermalization, the time-dependent spectral function smoothly changes from the vacuum spectral function to the thermal one. Because it involves a Fourier transform covering an infinite time range, the time-dependent spectral function deviates from the vacuum spectral function even before the quench, especially for low frequencies, where it does not exhibit a gap.
The techniques developed in this article enable the calculation of many quantities of interest in the thermalization process. For example, one could study the time dependence of the occupation number of field modes in the boundary gauge theory.

A A notion of time-dependent entropy
Here we summarize the arguments in [42] that lead to the entropy formula given in section 2. Let us consider particles with Bose statistics indexed by a macroscopic quantum number I. Also suppose that the states index by I have a degeneracy g I , and that there are n I particles with quantum number I. In both quantum physics and in statistical physics we would really have a probability distribution over n I , but we will assume that n I 1 and the fluctuations are much smaller than the mean so that we can ignore them. Thus we will take n I as given and will not integrate over the distribution of this variable.
Then, the number of ways of distributing the n I particles over the g I microstates is The Now let α be a particular micro state with macro occupation numbers {n I }. Assuming that all the microstates are equally probable, is the probability of the microstate α Thus we can write a formula for the von Neumann entropy where we used the fact that α P {n I } (α) = 1. Assuming that n I , g I 1 we can use Stirling's formula to write S = I (n I + g I ) ln(n I + g I ) − g I ln g I − n I ln n I , (A.5)

JHEP04(2013)069
where we also dropped the −1 in n I + g I − 1 and g I − 1. We can simplify further by factoring out g I : There are some cancellations, and then we can write Here n I /g I is the average occupation number of microstates.
To make contact with the notions introduced in section 2, we interpret g I = ρ(k, ω, T ), i.e., we take g I to be the time dependent density of states with macroscopic quantum numbers k, ω. We also take n I = n(k, ω, T )ρ(k, ω, T ), i.e. the total occupation number of modes with k, ω. We then get where s(k, ω, T ) = (n(k, ω, T ) + 1) ln (n(k, ω, T ) + 1) − (n(k, ω, T )) ln (n(k, ω, T )) . (A.9) To illustrate this formula, consider the simple case when the spectral density is a sum of (on-shell) delta functions: where the sum may account for any degeneracy by having several ω i (k) have the same value. Then in equilibrium: where n B (ω) = (e βω − 1) −1 is the Bose distribution. Thus n I = n B (ω) in equilibrium only depends on ω, not on k, and it is independent of the level density and the degeneracy of levels. The formula for the entropy then becomes: The quantity s accounts for the entropy associated with the different possible occupation numbers of one bosonic mode of energy ω.

B.1 Geodesic length in Euclidean AdS 3 and BTZ
We here explicitly derive the renormalised length of a geodesic between two boundary points in Euclidean BTZ and AdS 3 . This serves both to introduce the Wick rotation in this setting, as well as to derive a number of expressions that will be used in the main text.

JHEP04(2013)069
We parametrize Euclidean BTZ using Poincaré coordinates in which the metric reads where R < r, the surface r = R is the black hole horizon and r → ∞ represents the spacetime boundary. Since Euclidean AdS 3 is recovered from (B.1) by sending R → 0, we will work out explicitly the BTZ case alone and take the AdS limit only at the end. The Euclidean BTZ geodesics are given by for constants 0 < Γ 2 − < R 2 < Γ 2 + and x 0 , y 0 , λ 0 . As a function of r, we find the following branches (r ∈ [Γ + , ∞[), The separation between two boundary insertion points (x 1 , y 1 , r = ∞) and (x 2 , y 2 , r = ∞) is given by 9) and the length of the geodesic can be found to be in terms of a regularized AdS boundary at r = r 0 . In the limit r 0 → ∞ this reduces to ∆λ = 2 ln   2r 0

B.2 Two point functions of vacuum and thermal CFT
The Euclidean two-point functions can be computed in the geodesic approximation from the renormalised geodesic lengths that were obtained in the previous section. A suitable Wick rotation leads to the time-ordered and retarded two-point functions in Lorentzian signature. Below we also briefly comment on how to extract the spectral function and occupation number from these results. where the first can be obtained from the second in the R → 0 limit of vanishing black hole mass. Despite the fact that we have used a geodesic approximation, the result for the vacuum and thermal two-point function is exact, being fully constrained by conformal invariance (apart from an overall scaling).

B.2.2 Lorentzian two point functions
To go from Euclidean coordinates (x, y, r) to Lorentzian coordinates (x, t, r), we perform the Wick rotation y = it. It is well known that by Wick rotating the Euclidean two point function, one finds the time-ordered (Feynman) two point function For non integer scaling dimension ∆, this leads to 10 The identity, where G R/A denotes the retarded (advanced) two-point function, allows to deduce In our case, this results in where we take x = ∆x = x 2 − x 1 and t = ∆t = t 2 − t 1 , we can calculate the spectral function ρ(k, ω) and occupation number n(k, ω) from the identities ( and agree with those computed in [39] up to an overall normalization. From these expressions, we can deduce that 11 For half integer values of ∆, we can use the following expressions: and Γ(x + 1) = xΓ(x). (B.36)

JHEP04(2013)069
Integrating over all momenta gives the results while for the occupation number, using (B.31), we recover Bose-Einstein distribution C Notes on the case of integer Delta

C.1 Euclidean signature
Consider a conformal field theory in d dimensions. The Euclidean two-point function of a scalar operator with conformal weight ∆ is given by The Fourier transform of this function is where we have used the Gaussian integral the relation (valid for α > 0) and made the substitution s = 1/u This derivation is only strictly valid for d/2 > ∆ > 0, which is not of physical interest to us since we always have ∆ ≥ d/2. However, one can analytically continue the result in ∆ to all positive real values for which ∆ − d/2 is not zero or a positive integer. We could also use a regulator ξ, such that

C.2 Minkowski signature
We can perform a Wick rotation of the Euclidean two-point function to Minkoswki signature by taking y = e iθ t and continuously take θ from 0 to π/2− (where > 0 and → 0 is implied). We then find the time-ordered two-point function For non-integer ∆, we can rewrite this as For integer ∆, one needs to pay attention to contact terms that may arise as for example in (C. 23) In the mean text of this paper, we have assumed that ∆ is non-integer. One can find the results of integer ∆ by using the regularisation ∆ − δ as in the previous section, where the contact terms result in analytic terms in k 2 . One can show that this will change the expressions for the real part of the time-ordered and retarded propagators, but not the imaginairy part. Therefore, the expressions for the spectral function and occupation number are valid for all values of ∆.

D Complexified geodesics in AdS 3 and BTZ
In this appendix, we review the construction of complexified geodesics in AdS 3 and BTZ.

JHEP04(2013)069
Note that for complexified r-coordinate, going from the boundary (at infinity) to the boundary (at infinity) corresponds to |r(λ)| 2 going from the point at infinity in the complex plane back to the point at infinity. The direction in which this happens is determined by the parameter β. We now rewrite the affine parameter λ as a function of |r|: λ ± (|r|) = λ 0 ± arccosh |r| 2 e 2 − j 2 + cos 2 (β) , (D.9) so that the geodesic length can be expressed as ∆λ(|r|) = λ + (|r|) − λ − (|r|) = 2 ln |r| 2 e 2 − j 2 + cos 2 (β) + This result is independent of the regularisation parameter β, and it agrees with our expectation for the scaling of two-point functions in a vacuum CFT.
To be able to take the AdS limit R → 0, we have to consider the first class of solutions with (J 2 − E 2 + R 2 ) 2 > (2JR) 2 . In this first class, we can adopt the convenient notation 15) such that Λ + Λ − = J 2 R 2 and (Λ + − R 2 )(Λ − − R 2 ) = E 2 R 2 . We then find the complex solutions .
(D. 16) The separation on the boundary is given by where we have suppressed the spatial coordinates as they play no crucial role in the following analysis. For concreteness we work in the Schrodinger representation where the only operator depending on time is the density matrix. In the following we will take H 0 to be a generic Hamiltonian which is allowed to have explicit time dependence. Let us assume we have solved for the density matrix ρ when J = 0, meaning that we know the general solution to the equation Next we want to compute the retarded propagator, which can be obtained from the Feynman one through the identity G B,R (k, t, z; t ) = −θ(t − t ) G B,F (k, t, z; t ) + G * B,F (k, t, z; t ) . (G.8) To determine the Feynman propagator in the future lightcone, we should continue (G.7) to real time as a = z 2 +(y −y ) 2 → −(t−t ) 2 +z 2 +i . Since we need fractional powers of a we need to take into account the i factor to determine the phase a α → (−(t−t ) 2 +z 2 +i ) α = e iπα b, where b is the positive quantity b = (t − t ) 2 − z 2 . In continuing to real time the Feynman propagator is multiplied by a factor of i as a part of its definition. In this way we obtain . (G.10) The retarded bulk-to-boundary propagator can be obtained from the Feynman propagator through the identity where the θ(t − t − z) appears because G AdS B,F is purely imaginary for b < 0, which is outside the future lightcone.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.