Thermalization with chemical potentials, and higher spin black holes

We study the long time behaviour of local observables following a quantum quench in 1+1 dimensional conformal field theories possessing additional conserved charges besides the energy. We show that the expectation value of an arbitrary string of {\it local} observables supported on a finite interval exponentially approaches an equilibrium value. The equilibrium is characterized by a temperature and chemical potentials defined in terms of the quenched state. For an infinite number of commuting conserved charges, the equilibrium ensemble is a generalized Gibbs ensemble (GGE). We compute the thermalization rate in a systematic perturbation in the chemical potentials, using a new technique to sum over an infinite number of Feynman diagrams. The above technique also allows us to compute relaxation times for thermal Green's functions in the presence of an arbitrary number of chemical potentials. In the context of a higher spin (hs[\lambda]) holography, the partition function of the final equilibrium GGE is known to agree with that of a higher spin black hole. The thermalization rate from the CFT computed in our paper agrees with the quasinormal frequency of a scalar field in this black hole.


Introduction and Summary
The study of thermalization in closed interacting quantum systems has a long history (see, e.g. [1] for a review). It has been known ever since the celebrated work of Fermi, Pasta and Ulam (FPU) that interacting classical systems need not necessarily equilibrate. The question of finding sufficient conditions for thermalization in quantum systems is also an open one. Recently, the advent of holography has linked the issue of thermalization in strongly coupled quantum field theories to another important, classical, problem of black hole formation (see, e.g. [2,3] and references therein). In the latter setting too, the issue of gravitational collapse of a given matter distribution is rather nontrivial; indeed there is an interesting debate in the current literature (see, e.g., [4]) regarding the fate of perturbations in anti-de-Sitter spacetimes.
In this paper, we will focus on two-dimensional conformal field theories (CFTs) on an infinite line σ ∈ R. We will consider the system at t = 0 to be in a "quenched state" 4 Here |Bd is a conformal boundary state; the exponential factors cut off the UV modes to make the state normalizable. W n denote the additional conserved charges in the theory. 5 This choice of the quenched state is a generalization of that in [5] for which ǫ n = 0, for n > 3. The wavefunction for t > 0 is given by The questions we will explore, and answer, are: what is the long time behaviour of various observables in |ψ(t) ? In particular, does the expectation value of an operator (or a string of operators) approach a constant? If so, (i) is the constant value characterized by a thermodynamic equilibrium, and (ii) what is the rate of approach to the constant value? More generally, we would also address, and partially answer, the questions: how does the existence and rate of thermalization depend on the initial state and the choice of observables?
Thermalization We find in this paper that the expectation values of local observables (supported on a finite interval A : σ ∈ [−l/2, l/2]) asymptotically approach (see (12) for the precise statement) their values in an equilibrium ensemble, whose temperature and chemical potentials are related to the cutoff scales in (5) as follows β = 4ǫ 2 , µ n = 4ǫ n , n = 3, 4, ...
The relations (4) are uniquely dictated by the requirement that the expectation values of the conserved charges H, W 3 , W 4 , ... in the initial state match those in the mixed state (3) (see (52)). In the absence of the extra parameters ǫ n , n = 3, 4, ... this result is derived by the elegant method of conformal transformations [5]. In the presence of these parameters, this method is not available; in this paper, we deal with the extra exponential factors in terms of an infinite series and do a resummation. 4 In the original sense of the term, a quantum quench is defined as a sudden change from a hamiltonian H 0 to a hamiltonian H which governs further evolution for t ≥ 0. The system is assumed to be in the ground state of H 0 at t = 0, which serves as an initial state for subsequent dynamics; the dynamics is nontrivial since the initial state prepared this way is not an eigenstate of H. In this paper, as in [5], we will mean by a "quenched state" simply a pure state which is not an eigenstate of the Hamiltonian H. The kind of quenched state defined in (1) is sometimes said to describe a global quench or a homogeneous quench, as the state is translationally invariant. We will briefly mention inhomogeneous and local quenches in Section 6. 5 For the purposes of this paper, we will identify them with W n -charges of 2D CFT, n = 3, 4, ... (with W 2 ≡ H), although much of what we say will go through independent of this specific choice as long as these charges mutually commute and are defined from currents which are quasiprimary fields of the conformal algebra.
We emphasize that the thermalization we found above persists even when we have an integrable model with an infinite number of conserved charges. Relaxation in integrable systems has been found in recent years in the context of, e.g., (a) one-dimensional hardcore bosons [6], (b) transverse field Ising model [7], and (c) matrix quantum mechanics models [8]. The equilibrium ensembles in this context have been called a generalized Gibbs ensemble (GGE). Our present result on integrable conformal field theories adds to the list of these examples. Interestingly, the thermalization we find works even for free conformal field theories, e.g. a free scalar field theory. 6 With the above identification of parameters, we will rewrite the initial quenched state (1) henceforth as We find the following specific results: 1. Thermalization time scale for single local observables: We find that at large times ψ(t)|φ k (σ)|ψ(t) = Tr (φ k (0)ρ eqm (β, µ i )) + a k e −γ k t + ...
where φ k (σ) is an arbitrary quasiprimary field (labelled by an index k). Below we compute the thermalization exponent γ k in a perturbation in the chemical potentials and to linear order it is given by Here ∆ k = h k +h k is the scaling dimension and Q n,k are the (shifted) W n -charges (see (37) for the full definition) of the field φ k (in case of primary fields) or of the minimum-dimension field which appears in the conformal transformation of φ k . To obtain this result, we perform the infinite resummation mentioned below (4). At large times, the perturbation series for the one-point function in the chemical potentials exponentiates (see (37)), to give the corrected exponent in the above equation. In various related contexts, finite orders of perturbation terms in chemical potentials have been computed before [9,10,11]. Our finding in this paper is that at large times, there is a regularity among the various orders leading to an exponential function as in (6)

(see Section 2.2.2 for details).
Universality: In the case of zero chemical potentials, it has been noted in [12], that although the relaxation time τ k = πǫ 2 /(2∆ k ) = 2πβ/(∆ k ) is non-universal (in the sense that it depends on the specific initial state (1)), the ratio of relaxation times for two different fields, namely, τ k 1 /τ k 2 = ∆ k 2 /∆ k 1 is universal (it depends only on the CFT data and not on the initial state and is hence expected to be valid for a general class of initial states). In the presence of the additional cut-off parameters ǫ i , i = 3, ... in the initial state (1), the ratio τ k 1 /τ k 2 = γ k 2 /γ k 1 = (∆ k 2 + nμ n Q n,k 2 )/(∆ k 1 + nμ n Q n,k 1 ) is, however, not independent of the initial state. However, as we will briefly discuss in Section 6, for a large class of quench states (e.g. where the energy density is uniform outside of a domain of compact support) the β-dependence of τ k , in the absence of chemical potentials, can be understood as the dependence on the uniform energy density (see a related discussion in [13]). The time scales τ k , therefore, do have a limited form of universality in the sense that it depends on a rather robust feature of the initial state. Our calculations in this paper leads us to believe that this feature will continue in the presence of chemical potentials, in the sense that the additional dependence of the time scales 1/γ k on the µ n is fixed by the charge densities corresponding to the additional conserved charges. We hope to address this in [14].
2. Multiple local observables, reduced density matrix: Besides the one-point functions discussed above, it turns our that we can demonstrate thermalization of all operators in an interval A of length l. It is convenient to define a 'thermalization function' I A (t) [15] as Hereρ = ρ/ Trρ 2 denotes a 'square-normalized' density matrix. 78 We show below that at large times the thermalization function has the form where γ m refers to the exponent (7) for the operator φ m with minimum scaling dimension. 9 α(l) is computed as a power series inl which we find using the short interval expansion, valid forl ≪ 1, i.e. l ≪ β. Two immediate consequences of (9) are (i) Thermalization of an arbitrary string of operators: Note, from (9), that Since the square-normalized density matrices can be regarded as unit vectors (in the sense of footnote 7), and I A (t) can be regarded as the scalar productρ dyn,A (t)·ρ eqm,A , (10) clearly impliesρ 7 Note that operators in a Hilbert space H can themselves be regarded as vectors in H × H * ; under this interpretation Tr(A B) defines a positive definite scalar product. With this understanding, we will regard the hatted density matrices as unit vectors. 8 Throughout this paper, we will consider field theories with an infinite spatial extent. The entire Hilbert space is assumed to be of the form H A ⊗ HĀ. TrĀ implies tracing over HĀ. 9 We will assume here that the spectrum of such ∆'s is bounded below by a finite positive number. In case of a free scalar field theory, we can achieve this by considering a compactified target space.
This implies the following statement of thermalization for an arbitrary string of local operators (with σ 1 , σ 2 , ... ∈ A) (12) (ii) Long time behaviour of reduced density matrix: Carrying on with the interpretation of I A (t) as a scalar product, we can infer following asymptotic behaviour ofρ dyn (t) from (9): where Tr(Q 2 ) = 1, Tr(ρ eqm,A (β, µ i )Q) = 0. We will specify further properties ofQ later on.
Importance of local observables: In case of a free massless scalar field, it is easy to show that quantities like ψ(t)|α 2 1 α † 1 |ψ(t) perpetually oscillate and never reach a constant (see a related calculation in [8]). The modes α n represent Fourier modes and are non-local. Indeed, as [15,16,12] showed, in the absence of chemical potentials, the exponential term in (9) is e −2γm(t−l/2) and the thermalization sets in only after t exceeds l/2. Thus, for l = ∞, there is no thermalization, which is consistent with the above observation about perpetual oscillations. We expect the form e −2γm(t−l/2) to continue to hold in the presence of chemical potentials 10 , since the effect of the chemical potentials on the exponent γ k can be viewed as a shift of the anomalous dimension δ∆ k = nμ n Q n,k + O(μ 2 ) (see, e.g. (64)). This shows that, as in the case of zero chemical potentials, equilibration sets in only after t exceeds l/2. We will see a similar phenomena next in the context of a decay of perturbations to a thermal state.
3. Decay of perturbations to a thermal state: We compute (see Section 4 for details) the time-dependent two-point Green's function G + (t, l; β, µ) for two points spatially separated by a distance l. We find that for t, l, t − l ≫ β, the time-dependence is exponential, with the same exponent as in (6): Note that the above thermalization sets in for t > l. For t < l, the two-point function has an exponential decay in the spatial separation (see Section 4 and Figure 3). The computation of the above relaxation times in the presence of an arbitrary number of chemical potentials uses the technique, described above, of summing over an infinite number of Feynman diagrams, and is one of the main results of our paper.

Collapse to higher spin black holes:
In [17] the bulk dual to the time-dependent state (2) corresponding to initial condition (5), for large central charges, has been constructed in the case of zero chemical potentials. The dual geometry corresponds to one half of the eternal BTZ (black string) geometry, whose boundary represents an end-of-the-world brane. In [18] the result has been extended to the case of non-zero angular momentum and a Chern-Simons charge. In case of an infinite number of chemical potentials, a bulk dual to the equilibrium ensemble (3) has been identified, in the context of the Gaberdiel-Gopakumar hs(λ) theory [19], as a higher spin black hole with those chemical potentials [20]. It is natural to conjecture [18,8] that the time-development (2) should be dual to a collapse to this higher spin black hole. At late times, therefore, the thermalization exponent found above should correspond to the quasinormal frequency of the higher spin black hole. We find that (see Section 5 and [21]) this is indeed borne out in a specific example.
The plan of the paper is as follows. The results 1, 2, 3 and 4 above are described in Sections 2, 3, 4 and 5, respectively. The resummation of an infinite number of Feynman diagrams (corresponding to insertions of arbitrary number of chemical potential terms) is discussed in Section 2.2.2, which uses results in Appendix A. The calculation of the overlap of reduced density matrices in Section 3 needs the use of the short-interval expansion, which is described in Appendix B. In Section 6 we present our conclusions and make some remarks on inhomogeneous quench [14].

One-point functions
In this section we will consider the behaviour of the following one-point functions of a quasiprimary field φ k (σ) We will briefly recall how these are computed in the absence of the chemical potentials [5,22]. The first expectation value corresponds to the one-point function on a strip geometry, with complex coordinate w = σ + iτ , σ ∈ (−∞, ∞), τ ∈ (−β/4, β/4) where τ is eventually to be analytically continued to τ = it. This can be conformally transformed to an upper half plane by using the map For a primary field with h k =h k (of the form φ k (w,w) = ϕ k (w)ϕ k (w)), this procedure gives 11 (for other primary fields, the one-point function vanishes) The subscripts str, cyl will denote a 'strip' and a 'cylinder', respectively.
We have used the following result for the one-point function on the UHP: which follows by using the method of images where the antiholomorphic factor of φ k (z,z) on the upper half plane at the point (z,z) is mapped (up to a constant) to the holomorphic ϕ * k 12 on the lower half plane at the image point (z ′ ,z ′ ) with z ′ =z,z ′ = z [23,22]. In the above a k , A k are known numerical constants. Note that (19) so that in the large time limit we have The second, thermal, expectation value in (15), for µ n = 0, corresponds to a cylindrical geometry in the w-plane, with τ = 0 identified with τ = β. By using the same conformal map (16) this can be transformed to a one-point function on the plane. For a primary field the latter vanishes. Hence (6) is trivially satisfied.
For a quasiprimary field φ k , its conformal transformation generates additional terms, including possibly a c-number term c k (e.g. the Schwarzian derivative term for φ k = T ww ) and generically lower order operators. The c-number term does not distinguish between a plane and an UHP. This leads to the following overall result (for µ n = 0): where ∆ k now is the scaling dimension of the minimum-dimension operator in a T (z 1 )φ k (z) OPE. This is clearly of the general form (6) for µ n = 0.
We now turn to a discussion of these expectation values (15) in the presence of chemical potentials µ n , n = 3, 4, ..., as in (5) and (3). We will denote the new conserved currents as W n (w) andW n (w), n = 3, 4, .... The conserved charge, W n , is defined as Here the contour Γ is taken to be a τ = constant line along which dw 1 = dw 1 = dσ. Under the conformal transformation (16) to the plane/UHP, the holomorphic part of the contour integral becomes 12 We distinguish ϕ * k from ϕ k to allow for charge conjugation.
where the a n,n−2m denote the mixing of W n (z 1 ) with lower order W -currents under conformal transformations [24]. The contour Γ 1 is an image of the contour Γ onto the plane. The expression for the antiholomorphic part W n | antihol is similar. As mentioned before, in this paper we will regard the W n as conserved charges of a Walgebra, although the results we derive will be equally valid as long as these charges, together with H, form a mutually commuting set, and the currents (W n (w),W n (w)) are quasiprimary fields.

One-point function on the cylinder with chemical potentials
For simplicity we first consider the equilibrium expectation value in (15). Unfortunately, unlike the thermal factor above, the factor e − n µnWn in (3) cannot be dealt with in terms of a conformal map. We will, therefore, treat this factor as an operator insertion, and write We will now illustrate how to compute this for a single chemical potential, say µ 3 , using perturbation theory Feynman diagrams: 13 The first term in the above expression is the constant c k that we already encountered in (21). For a holomorphic primary field φ k , the second, O(µ 3 ), term, transformed on to the plane, gives Here we have used the contour representations (22) and (23). The correlator inside the second integral obviously vanishes (it factorizes into a holomorphic and an antiholomorphic one-point functions, leading to a vanishing connected part). The first integral vanishes unless φ k = W 3 (this uses the orthogonality of the basis of quasiprimary fields). In the latter case, using the integral evaluates to c/(90z 3 ); combining with the factor of z 3 outside (h k = 3 in this case) we get a z-independent constant, as we must, because of translational invariance on the plane. With an antiholomorphic primary field φ k , the calculation is isomorphic. For a primary field with nonvanishing h k ,h k the result vanishes. For quasiprimary φ k , as well as for other W n charges, the conformal transformation to the plane additionally generates lower order operators (see, e.g. (23))), each of which can be dealt with as in (26). The result is a finite constant which we will denote as W n φ k (w,w) = c n,k (this will be non-vanishing only for special choices of φ k , e.g. φ k = W n ). As explained above, for n = 3 and φ k (w,w) = W 3 (w), c n,k = −2πc/(90β 2 ).
In a similar fashion, the O(µ 2 3 ) term in (25) can be transformed to the plane. Again, we present the explicit expression for the simple case of a holomorphic primary field φ k .
For holomorphic quasiprimary φ k , additional, similar, terms appear due to the generation of lower order operators under conformal transformation to the plane. Only the holomorphic correlator survives (as in the O(µ 3 ) calculation). Thus, e.g. if φ k = T (z), the stress tensor, we have Again, after performing the integration over z 1 and z 2 , we obtain a z-independent constant, as we must. The analysis of more general fields φ k and two arbitrary W -charges is straightforwardly generalizable. The result is a finite constant (can be zero for a particular φ k ) which we denote as W m W n φ k (w,w) = c mn,k Note that in (27) the result does not depend on the location of the contours Γ 1 , Γ 2 on the plane, since the W -currents are conserved. Summarizing, we get

One-point function on the strip with chemical potentials
Similarly to the previous subsection, we will treat the µ-deformations in (5) as operator insertions: As before, we begin by illustrating the calculation of this quantity with the simplest case of a single chemical potential µ 3 , using perturbation theory Feynman diagrams: The { , } denotes an anticommutator. The operator ordering implies the following: when W 3 appears on the left of φ k (w,w), e.g., in W 3 φ k (w,w) , the integration contour (22) for W 3 on the strip lies above the point (w,w); similarly when W 3 appears on the right of φ k (w,w), e.g. in φ k (w,w)W 3 , the contour for W 3 is below the point (w,w). The first, µ-independent, term in the above expansion is already calculated in (21).

O(µ n ) Calculation
As before, we find it convenient to use the conformal transformation (16). The correlator on the strip then reduces to that on the UHP, as in the µ = 0 case before. For a holomorphic primary field φ k , this gives where the operator ordering explained above implies that the contour Γ 1 lies to the left of the point (z,z) on the UHP. Now, in the analogous calculation (26), the second connected correlator on the complex plane vanished because of factorization into one-point functions.
Correlators on the UHP are, however, related to those on the plane by the method of images (an example of which we saw in (18)). In particular,W 3 at the point (z 1 ,z 1 ) on the UHP becomes the holomorphic operator W * 3 = −W 3 on the LHP at the point (z ′ 1 ,z ′ 1 ) with z ′ 1 =z 1 [23,22]. The contour Γ 1 gets mapped to its mirror image Γ ′ 1 on the lower half plane. With this, we get On the complex plane, the contour Γ 1 on the UHP can be deformed to Γ ′ 1 on the LHP, hence the two contours simply yield a factor of 2. In fact, combining with the other ordering, and applying a similar reasoning, we get an overall factor of 4. Thus, combining with results from Section 2.1, we get, for holomorphic primary fields A similar statement is true for an antiholomorphic primary field. Let us turn now to primary fields φ k (w,w) with h k ,h k = 0 (of the form φ k (w,w) = ϕ k (w)ϕ k (w), as discussed before in the context of (18)). In the cylinder calculation in Section 2.1 the µ-corrections for these vanished. In the present case, they are non-zero for operators of the form φ k (w,w)= ϕ k (w)φ k (w), with h k =h k (as in (17)). After conformally transforming to the UHP, we regardφ k on the UHP as ϕ * k at the image point on the LHP (up to a constant). Combining with the arguments used for the holomorphic operators, we eventually get Figure 1: Various contours needed to compute the W n insertions in (30). At late times, the insertion of each contour, irrespective of the position of the contour, amounts to insertion of a given factor linear in t. This allows to resum arbitrary orders of arbitrary W n -charge insertions, leading to the exponential time-dependence as in (6). See Figure 2 for more.
The ratio of correlators inside the integral is given by where q 3 is the W 3 -charge of the field φ k . Integrals of the kind (34) are discussed in detail in Appendix A.2. The final result (see (71)) is that the O(µ 3 ) correction, in the long time limit (20), is given by (using that all four contours Γ 1 ,Γ 1 , Γ ′ 1 ,Γ ′ 1 contribute equally, cancelling the Up to O(µ 3 ), it agrees with (7).
In case of a quasiprimary field φ k (w,w), it mixes, under conformal transformation to the plane, with lower dimension operators. The most relevant operator among these, which is of the form ϕ k (z)ϕ k (z), is then to be used in (34) for obtaining the dominant time-dependence; in that case ∆ k , Q 3,k refer to this operator (rather than to the original φ k ).
For higher W n charges, the currents W n (w) are typically quasiprimary, and hence they mix with lower order W m (z) under conformal transformation to the UHP. Thus the O(µ n ) correction to the dynamical one-point function φ k dyn is a linear combination of terms of the form (68) (weighted by a set of coefficients a n,m , as in (37) below). Collecting all this, the O(µ) correction with all chemical potentials is given by m=0 a n,m i n−2m (2π) n−2m−2 q n−2m,k = i n (2π) n−2 2q n,k + i n−2 (2π) n−4 a n,2 2q n−2,k + ..., Note that for W 3 deformations, the expression for Q 3 as in (36) corresponds only to the first term in the above series expression for Q n . This is because the W 3 current is a primary field and does not mix with any lower W current under a conformal transformation. From n = 4 onwards, the additional terms in Q n,k 's represent the mixing of W n currents with W n−2m under conformal transformations.

Higher order µ-corrections
Let us first consider that O(µ 2 n ) correction: Again, for holomorphic (or antiholomorphic) primary fields φ k (w), it is straightforward to generalize (33) to this order.
For a primary field of the form φ k (w,w) = ϕ k (w)ϕ k (w), proceeding as in the previous subsection, we get The essential ingredient in this calculation is By repeating the strategy of (75), we get where we have first used the W m (z 2 )ϕ k (z) OPE, and then the W n (z 1 )ϕ k (z) OPE. In a manner similar to that in Appendix A.2, we conclude the following structure of I nm (z, z ′ ): Note that at late times t ≫ β, ([log(−z ′ ) − log(−z)] → 2(2πt)/β and dominates over the constant term (the precise sense is that of (48)). Similar to Appendix A.2, the 4 × 4 = 16 locations of the contour-pairs (Γ 1 , , ...., all contribute equally, therefore converting (µ n /4)(µ m /4) → µ n µ m . Combining all these, we get (40). The charges q n that are defined by the W n ϕ OPE and appear in (42), get multiplied by some constants 14 and shifted by lower W n−2k charges to give the Q n in (42), as in (37).

Arbitrary orders and Exponentiation:
It is straightforward to generalize the above O(μ 2 ) calculation to higher orders in the perturbation in chemical potentials. Thus, at the order r i=1 µ n i , there are r insertions of W-currents, leading to integrals of the form I n 1 n 2 ...nr (z, z ′ |Γ 1 , Γ 2 , ..., Γ r ) ≡ Again, repeating the strategy of (75), we get the following leading (viz. (log) r ) contribution (see (48) for the definition of the leading-log contribution) where we have first used the W nr (z r )ϕ k (z) OPE, then W n r−1 (z r−1 )ϕ k (z) OPE, etc. As in the O(µ 2 ) calculation above, we obtain the following behaviour at late times basic ingredient for the exponentiation we are going to find. Furthermore, it is easy to see that the leading log contribution is the same irrespective of where each contour Γ i is placed (out of 4 possible choices, e.g. Γ 1 , Γ ′ 1 ,Γ 1 ,Γ ′ 1 in Figure 1). As before we must combine the contribution of all positions of the contours, which, therefore, amounts to multiplying the result for (44) by 4 r which converts the original coefficients coming from exp[− n µ n W n /4] as follows This is the second basic factor leading to the exponentiation. Combining all these, and incorporating some additional constants (see footnote 14) we get the following, leading, order (µ n 1 ...µ nr ) contribution Once again, the constants Q n are related to the q n as in (37)) in a manner similar to the O(μ) and the O(μ 2 ) calculation above. We note that the leading log contribution used in this paper can be isolated by considering a scaling The second term in (47), or for that matter, in (40), is subleading at large times in the sense of this scaling. Using the above results, we now have, for primary fields of the form φ k (w,w) = ϕ k (w)ϕ k (w) where we have absorbed some constant factors in a k . γ k is given by (7); Q n,k are the shifted W n charges of φ k as defined in (37). The proof of the above equation for general quasiprimary operators φ k (w,w) works out much the same way as in case of the O(µ) terms, as discussed in Section 2.2.1. We emphasize that it is only the leading contributions at large times which we have proved here to exponentiate. Thus, we do not claim that the constant terms marked "const" in the above equation are all the same. As we have remarked before, the leading contributions can be isolated using the scaling mentioned in (48). The schematics of the above calculation is explained in the Figure 2.  (17) without chemical potentials (the shading indicates the boundary of the upper half plane). The second term represents the O(µ n ) correction, which involves one insertion of a W n -charge (which is an integral over the z 1 -contour. As explained in the text, at long times, this insertion amounts to multiplying the zero order term by a term of the form f n log(z), where f n is described in (37). The third term represents insertion of two such W -charges; as we explained in the text (see (40) and below), each insertion again amounts to multiplying by the factor mentioned above, along with a factor of 1 2! . The pattern continues, to ensure an exponentiation to G 0 (z) z n fn , as in (49). Since at long times G 0 (z) ∼ e −γ (0) k t (see (17)), and z ∼ e −2πt/β , adding the chemical potentials amount to a shift of the exponent γ (0) k → γ k as in (6).

Calculation of I(t)
Let us rewrite the expression for the thermalization function I(t) (8) in the form Z sc ≡ Tr(ρ dyn,A (t)ρ eqm,A (β, µ)),Ẑ sc = Z sc /(Z s Z c ) Z ss ≡ Tr(ρ dyn,A (t)ρ dyn,A (t)),Ẑ ss = Z ss /Z 2 s , Z cc ≡ Tr(ρ eqm,A (β, µ)ρ eqm,A (β, µ)),Ẑ cc = Z cc /Z 2 c , Z s = Tr(ρ dyn (t)) = ψ 0 |ψ 0 , Z c = Tr(ρ β,µ ) In Appendix B we explain how to compute I(t) using the short interval expansion, valid when the length of the interval l is small compared with the other time scales β and t in the problem. We reproduce the main formula (79) for our purpose, where we explicitly denote the dependencies on the length l of the interval, the inverse temperature β and the chemical potentials µ (the dependence on β on the RHS is implicit; the one-point functions depend on both β and µ-see Section 2).
It is understood, for the logic of the short interval expansion to go through, that all contours which represent insertion of the W -charges (see Fig 1) are drawn outside of the small disc-like region of both sheets of

Proof of thermalization
Using the short-interval expansion above, and the long time behaviour of one-point functions from Section 2), it is easy to prove that the system thermalizes in the sense of (10) or (11).
To prove this, note that it is only the holomorphic (or antiholomorphic) fields φ k which possibly have non-zero expectation values in the long time limit (20). For these fields, the one-point functions on the cylinder and on the strip agree (see (21), (33), (39) ). By virtue of (51), we therefore have in the long time limit Z sc = Z ss = Z cc . Hence using the expression (50) for the thermalization function we get I(t → ∞) = 1 which proves (10) and consequently (11).
The above-mentioned equality of one-point functions between the strip and cylinder geometries for holomorphic (or antiholomorphic) fields imply the same for the conserved W n -(orW n )-currents. This, therefore, proves that Note that in proving this, we have used the correspondence (4) between the parameters of the initial state and the putative equilibrium state. The above equation, therefore, proves the correspondence (4).

Thermalization rate
To evaluate the rate of approach of I(t) to its asymptotic value 1, we organize the terms in Z sc ,Ẑ ss ,Ẑ cc as followŝ  (53) where a, b, ... denote descendents of the identity operator, k labels other primaries (than the identity) and their descendents.Ĉ ≡ C/C 0,0 .

µ = 0
Let us first consider the case of zero chemical potentials. Using the results in Sections 2, and Appendices A and B.1, we get To this order, it is easy to see that the contribution to I(t) from descendents of identity, demarcated by a T , a TT , vanishes. The leading contribution to I(t), demarcated by a k , occurs only inẐ ss and comes from ( φ m (z,z) str ) 2 for which h k is the minimum (= h m ) (this could be a field which appears after a conformal transformation of the original quasiprimary field). The time-dependence shown of S ss 2 comes from (17). Using this, we get This is of the form (9) for µ = 0, with The discarded terms in (55) are faster transients. This proves (9) for zero chemical potential. This result has already appeared in [15]. 15

µ = 0
The generalization of the above result to the case of non-zero chemical potentials is straightforward. Once again, the dominant time-dependence arises from ( φ m (z,z) µ str ) 2 in the S ss 2 orẐ ss . The time-dependence (9) follows by using (49) in S ss 2 .

Properties ofQ
From the asymptotic behaviour (9) of the thermalization function we indicated the asymptotic behaviour (13) of the dynamical reduced density matrixρ dyn (t). By using the long time behaviour of the one-point functions (6), we can easily deduce the following dominant behaviour of overlaps ofQ with various quasiprimary fields at late times Tr(Qφ k (t)) ∝ e −(γ k −γm)t , Tr(Qφ m (t)) → constant. 15 Our exponent differs from Cardy's value by a factor of 2.

Decay of perturbations of a thermal state
We found in the previous sections that the long time behaviour of the reduced density matrix ρ dyn,A (t) resembles that of a thermal ensemble plus a small deformation which decays exponentially. We will find in the next section that the thermal ensemble (or more accurately the generalized Gibbs ensemble) corresponds to a (higher spin) black hole geometry in the bulk. The small perturbation of the equilibrium ensemble is thus expected to correspond to a small deformation of the black hole geometry. Consequently, the exponential decay of the deformation in the CFT should correspond to a 'ringing-down' or a quasinormal mode in the bulk. We will address the above issue in the next section which deals with bulk geometry. However, in order to make the correspondence of the above paragraph more precise, in this section we will directly present a CFT computation of the decay of a perturbation to a thermal state. Note that this computation is, in principle, different from the exponential decay of the one-point function in the quenched state, (6). However, what we will find is that the long time behaviour (6) of an operator φ k (0, t) in the quenched state is the same as that of its two-point function (57) in the thermal state (3) (with chemical potentials). The latter measures the thermal decay of a perturbation and is more directly related to a black hole quasinormal mode. Throughout this section, we will assume that the conformal dimensions of φ k satisfy h k =h k .
We define the thermal two-point function as 16 By the techniques developed in the earlier sections, a computation of this quantity amounts to calculating the following correlator on the plane where the µ n -deformations are understood as an infinite series of contours as in the previous section. For µ = 0, the above two-point function is standard. Including the Jacobian of transformation, we get which clearly matches the long time behaviour of the one-point function (6) in the quenched state for µ = 0. Here ∆ k = 2h k .
In the above, we considered the thermal Green's function for two points which are both at the same spatial point σ = 0. It is easy to compute the Green's function when the two points are spatially separated by a distance l, say with σ 1 = l and σ 2 = 0. We get The coordinates of the two points, in the notation of (58) are modified here to z = ie 2π(l−t)/β ,z = −ie 2π(l+t)/β , y = i,ȳ = −i. The prefactor with the square bracket comes from the Jacobian of the transformation from the cylinder to the plane. The behaviour of the Green's function is shown in Figure 3. It is important to note that the exponential decay, found in (6) shows up only for time scales t ≫ l.
|G + (t,l)| l=6 l=8 t Figure 3: Plots of the thermal Green's function G + (t, l; β, 0) for β = 2π, ∆ k = 1.5. The curve on the left (blue) is for l = 6, and the curve on the right (orange) is for l = 8. Note that the exponential decay in time occurs for times larger than l.
The effect of turning on the chemical potentials can be dealt with as in the previous sections. At O(µ n ), we will have, as before, a holomorphic contribution and an antiholomorphic contribution. The former is proportional to As we see, the structure of the integral is the same as in the previous section. As before, logarithmic terms appear in the above integrals which give the leading, linear, t-dependence. Similar remarks also apply to the antiholomorphic contour. Since the calculations are very similar to those in the previous two sections, we do not provide all details. By resumming the series over the infinite number of contours, we find in a straightforward fashion that where b(µ) is time-independent, and is of the form b(µ) = 1 + O(µ). This long time decay is the same as that of the one-point function (6) in the quenched state, as claimed above. For points separated by a distance l, the above exponential decay shows up for t ≫ l, as in (60).
In the above, we have discussed the two-point function in real space. It is straightforward to convert the result (60) without chemical potentials to Fourier space, which develops poles at Our results in (6) can be interpreted as a shift, caused by the presence of the chemical potentials µ n , of the dominant pole ω k,0 | µ=0 to where the notation is the same as that of (6). In this paper we will not address the question of the shift of the subdominant poles ω k,m (for m = 1, 2, ...) due to chemical potentials (the current status of these can be found in [26,9,10]). Two-point functions of the kind (57), for a single chemical potential µ 3 , and up to order µ 2 3 , have appeared earlier in [9] (calculations up to O(µ 5 3 ) have appeared in [10]). What we find in our paper is that at large times, the perturbation series in µ n , up to all orders in all chemical potentials, can be resummed, to yield the leading correction to the thermalization rate in the presence of chemical potentials.
At a technical level, the one-point function in the quenched state corresponds to a onepoint function in a geometry with a boundary, and for operators considered here, these turn into a two-point function on the plane, by virtue of the method of images. The thermal decay naturally involves a two-point function on the plane 17 and agrees with the above two-point function at late times.

Holography and higher spin black holes
Zero chemical potential: As remarked in the Introduction, a global quantum quench described by an initial state of the form (5), for large central charges and zero chemical potentials, has been shown in [17,18] to be dual to one half of the eternal BTZ (black string) geometry, whose boundary represents an end-of-the-world brane.
In an independent development, it was found in [27] that the quasinormal mode of a scalar field Φ k (σ, t, z) of mass m in a BTZ background (dual to a CFT operator φ k of dimension ∆ k ≡ 1 + √ 1 + m 2 ) is of the form exp[−2π∆t/β] at large times. This time-dependence agrees with the CFT exponent in (60) exactly. This shows that the exponential decay of a CFT perturbation to a thermal state corresponds to the decay of the corresponding scalar field in the bulk geometry. This result has been extended to higher spin fields in the BTZ background in [28].
Non-zero chemical potentials: In case the CFT has additional conserved charges, in particular if it has a representation of a W ∞ algebra (and consequently the hs(λ) algebra [19]), then the bulk dual corresponding to those conserved charges have been conjectured to be the conserved higher spin charges of higher spin gravity. In particular, [20] have shown that if one interprets the grand canonical ensemble (4) (more generally, the GGE) in the framework of an hs(λ) representation, then the bulk dual corresponds to a higher spin black hole.
Thus, we would like to conjecture that the bulk dual of the quantum quench with chemical potentials, would correspond to a gravitational collapse to a higher spin black hole.
As an important consistency check, by analogy with the case with zero potential, in the present case too, the leading quasinormal mode (QNM) of a scalar field Φ k (σ, t, z) should have a time-dependence given by (62). Following the results in [26] (see also [9,10,21]) 18 we find that at late times t ≫ β the QNM for the hs(λ) scalar field Φ + behaves, up to O(µ 3 ), as e −iω k,0 t , where where the index k here refers to the operator φ k dual to the scalar field Φ + . Noting that for this operator we have ∆ k = 1 + λ, and Q 3,k = 1 3 (1 + λ)(2 + λ) [9,10], we see that the QNM frequency ω k,0 agrees, to the relevant order, with the pole (64) of the thermal 2-point function which, in turn, is related to the thermalization exponent by the relation ω k,0 = −iγ k , with γ k given in (6).

Discussion
In this paper, 2D conformal field theories were considered with additional conserved charges besides the energy. We probed non-equilibrium physics starting from global quenches described by conformal boundary states modified by multiple UV cut-off parameters (1). It was found that local observables in such a state thermalize to an equilibrium described by a grand canonical ensemble (4) with temperature and chemical potentials related to the cut-off parameters. We computed the thermalization rate for various observables, including the reduced density matrix for an interval. It was found that the same rate appears also in the long time decay of two-point functions in equilibrium (see (6) and (14)). In the context where the number of conserved charges is infinite, and they are identified with commuting W ∞ charges, the equilibrium ensemble (a generalized Gibbs ensemble, GGE) corresponds to a higher spin black hole [20]. We found that the thermalization rate found above agrees with the leading quasinormal frequency of the higher spin black hole; this constitutes an additional, dynamical, evidence for the holographic correspondence between the global quenches in this paper and the evolution into the higher spin black hole.
One of the main technical advances made in this paper is the resummation of leadinglog terms at large times, presented in Section 2.2.2, which leads to exponentiation of the perturbation series, leading to the thermalization rate, presented in (6), (49), as a function of chemical potentials. This allows us to also compute the effect of chemical potentials on the relaxation times of thermal Green's functions. Another technical advance consists of the computation of the long-time reduced density matrix (9), using a short-interval expansion, which allows us to prove thermalization of an arbitrary string of local observables.
One might wonder whether the results presented in this paper are tied to the use of translationally invariant quenched states such as (1), whose energy density and various charge densities are uniform. We will address the question of inhomogeneous quench in a forthcoming paper [14], both in the CFT and in the holographic dual, using the methods of [29] where we create an inhomogeneous energy density by applying conformal transformations. It turns out [14] that if the initial state has inhomogeneities in a compact domain and has uniform energy densities outside, local observables again thermalize asymptotically with exponents governed by the uniform densities. Other important issues involve local quenches (see, e.g. [30]), and compact spatial dimensions. The issue of thermalization when space is compact is quite subtle. It has been shown in [15] that at large times one can have the phenomenon of revival (observables effectively returning to their initial values). The dynamical entanglement entropy for a quantum quench in a space with boundaries is an interesting, related, issue; we hope to come back to this in a forthcoming publication [31].

strip: In this case
where A TT , a TT are constants as in (17) and (18).
Case k = descendent of other primaries: In this case, 1. cylinder: The one-point function vanishes as in the case of primaries. 2. strip: The one-point function can be related to one-point function of primaries which is dealt with above.
In this equation we have displayed the principal value of the relevant integrals (the discontinuity (70) tells us the coefficient of the log term or the linear t term). However, we would like to understand the above results more simply, by using the W n (z 1 )ϕ k (z) OPE which is of the form: where ϕ k,i (z) is of dimension h k + i. 19 Using this, we get an expansion for the connected 3-point function of the form: Performing the integral in (68), The ellipsis in each round bracket represents terms with higher powers of 1/(z 1 − z) (up to a maximum of (z 1 − z) −n ); successive round brackets themselves are arranged in higher inverse powers of z − z ′ . Using the W n (z 1 )ϕ * k (z ′ ) OPE in a similar fashion and using the symmetry property g n (z 1 , z, z ′ ) = (−1) n g n (z 1 , z ′ , z) we can arrive at a general structure where R n (z, z ′ ) = (−1) n−1 R n (z ′ , z) is of the form P n−1 (z, z ′ )/(z − z ′ ) n−1 (P n−1 (z, z ′ ) is a homogeneous symmetric polynomial of degree zero). See the explicit form of R n for n = 3 in (69). The omitted terms are all ratios of homogeneous polynomials in (z, z ′ ) of the same degree in the numerator and in the denominator. This implies that we have, in the long time limit (20) I n (z, z ′ |Γ 1 ) = I 3 (z, z ′ |Γ 1 ) = 2q n,k (2π/β)t + q n,k × const + O(e −2πt/β ) which, of course, agrees with (71). Note that the dominant time-dependence 2q n,k t(2π/β) comes from the long-time limit of the coefficient R n (z, z ′ ) of the log terms, which can be read off from the discontinuity I n (z, z ′ |Γ 1 ) − I n (z, z ′ |Γ 1 ) (see (70)). Now, the contour Γ 1 −Γ 1 dz 1 can be deformed to a very small circle Γ z dz 1 around the point z; therefore the leading long-time behaviour R (0) n (z, z ′ ) can be derived by using the leading OPE singularity in (72) and computing the residue at z 1 = z:

B Short interval expansion
In this section we will explain a formalism suitable for computing partition functions of the kind that appear in (50). For convenience we will first compute these quantities in Euclidean time τ = it and later analytically continue back to Lorentzian time. With this, each of the expressions Z sc , Z ss , Z cc is of the form Tr(ρ A,1 ρ A,2 ) = geometry 1 Dϕ 1 geometry 2 where S[ϕ] represents the action for the CFT (with fields ϕ) and the delta-functional in the measure represents a gluing condition between a geometry '1' and a geometry '2' along a 'cut' which is the location, at a particular time τ , of the spatial interval A : σ ∈ (−l/2, l/2) 20 For Z ss , both geometries are that of a strip of the Euclidean plane described by complex coordinates (w,w) = σ ± iτ defined by boundaries at τ = ±β/4 with boundary conditions determined by the boundary state |Bd introduced in (5). For Z cc , both geometries are that of a cylinder cut of the Euclidean plane with identified boundaries at τ = −β/4, 3β/4. The geometries for both Z ss and Z cc are familiar from calculations of Entanglement Renyi entropy (of order 2) and can be calculated from appropriate correlation functions of twist fields [32] which exchange two identical geometries. For Z sc , the two glued geometries are different (that of a strip and a cylinder), hence the method of twist operators do not apply in a straightforward fashion. (See Figure 4). In this paper, we will therefore, employ the method of the short interval expansion. Figure 4: Two different geometries, the strip and the cylinder, glued along the cut as described in the text.
The method of the short interval expansion allows us to compute the functional integral over this geometry by replacing a small tube enclosing the two glued cuts by a complete basis of operators φ k1 ⊗ φ k2 where the operators live in the two Hilbert spaces.
The idea of the short interval expansion [33] is as follows. To begin, we express the functional integral (76) as an overlap of two wavefunctions in H 1 ⊗ H 2 , as follows Z 12 = Tr(ρ A,1 ρ A,2 ) = ψ out |ψ in = w 1 ∈D 1 Here D 1 (respectively, D 2 ) is a small disc drawn around the cut in geometry 1 (respectively, geometry 2). Note that only |ψ in depends on the gluing condition since the delta functional in the measure does not affect |ψ out . The basic point of the short interval is that in the limit when the length l of the cut is small compared with the characterizing length scale of the geometries (in our case, when l ≪ β), the wavefunction ψ in [ϕ 1 , ϕ 2 ] becomes jointly localized at the centre (w 1 ,w 1 ) of the disc D 1 and at the centre (w 2 ,w 2 ) of the disc D 2 21 , and hence can be expanded in terms of local operators, as follows Here k 1 , k 2 label a complete basis of quasiprimary operators of the CFT Hilbert space. Each term in the sum represents a factorized wavefunction (between geometries 1 and 2), which, therefore, gives 22Ẑ sc = k 1 ,k 2 C k 1 ,k 2 φ k 1 (w 1 ,w 1 ) str φ k 2 (w 2 ,w 2 ) cyl , Z ss = k 1 ,k 2 C k 1 ,k 2 φ k 1 (w 1 ,w 1 ) str φ k 2 (w 2 ,w 2 ) str , Z cc = k 1 ,k 2 C k 1 ,k 2 φ k 1 (w 1 ,w 1 ) cyl φ k 2 (w 2 ,w 2 ) cyl (79) Here the subscripts str and cyl refer to "strip", and "cylinder" respectively. The one-point functions are evaluated on the respective geometries without any cut (see Section 2 for more details). The glued functional integral (76), (77) is recovered by summing over k 1 , k 2 with the coefficients C k 1 ,k 2 ; , as clear from (79) these are determined by the gluing condition and depend on the size of the cut [33] (see Section B.1 for more details).

B.1
The coefficients C k 1 ,k 2 As explained in [33] (see also Section B), the coefficients C k 1 ,k 2 are determined by the equation where C 2 represents two infinite planes glued along a cut A, Z 2 is the functional integral such a glued geometry and Z 1 is the functional integral over a single plane. This equation can be easily proved by inserting quasiprimary a operator at infinity in each plane in an equation like (76) or (77). The two point function in the glued geometry is to be determined by using the uniformizing map: The normalization constants n k are determined by the following orthogonality condition of the quasiprimary operators φ k 1 (z 1 ,z 1 )φ k 2 (z 2 ,z 2 ) C = n k 1 δ k 1 ,k 2 z h k 1 +h k 2 12zh k 1 +h k 2 12 (82) where n k 1 is a normalization constant. Note that C k 1 ,k 2 = C k 2 ,k 1 . Below we will use the notationĈ k 1 ,k 2 = C k 1 ,k 2 /C 0,0 Case (k 1 , k 2 ) = (0, 0): We will denote the identity operator as φ 0 = 1. It is obvious that All other C k,0 vanish as they are proportional to a one-point function of a primary operator on the Riemann surface (and hence to that on the complex plane). Case (k 1 , k 2 ) = (primary, primary): In case φ k 1 , φ k 2 are primary operators, (80) giveŝ C k 1 ,k 2 = 1 n k 1 δ k 1 ,k 2 le iπ/2 4 2(h k 1 +h k 1 ) Case (k 1 , k 2 ) = (descendent, descendent): In case φ k 1 is of the form L −n 1 L −n 2 ...L −m 1L −m 2 ...φ l 1 and φ k 2 is of the form L −r 1 L −r 2 ...L −s 1L −s 2 ...φ l 2 , we can show that C k 1 ,k 2 = δ l 1 ,l 2 δ n, r δ m, s A(n 1 , n 2 , ..., m 1 , m 2 , ...; r 1 , r 2 , ..., s 1 , s 2 , ...) l 2(h k 1 +h k 1 ) , where A(...) is a numerical coefficient. 23 Here Λ(z) = : T T :(z) − 3 10 ∂ 2 z T is the level 4 quasiprimary descendent of the identity.